0

I am trying to create a time series object from extracted climate data (NEX-GDDP) using the Google Earth Engine (GEE). The data is daily meteorological data, and in the attached file, the data for January, 2005, is collected over an area of interest. The images from GEE are stored in the GeoTIFF as bands (numbered 1-31), and now I am struggling to get these individual bands into a dataset, and add a time dimension to the file. GEE will not export for more than ten years, so my idea is to create yearly files, which, when saved locally, will be merged (concatenated) on the lat/lon and time dimensions.

I am using Python in a Windows environment, so I am a bit limited (for example, I can't use cdo as this is a Linux based library), and I think that what I would like to do is possible with xarray, but I am missing the (learning) resources to solve this problem with code.

The image shows the xarray view of the metadata: https://i.stack.imgur.com/odWZK.png

Showing one band with Test data Missing the time dimension: https://i.stack.imgur.com/MIp4I.png

Can you provide any link to a training module on NetCDF?

Update:

import rioxarray as rxr

test = rxr.open_rasterio(filename, masked=True)
test

xarray.DataArrayband: 31y: 15x: 13
array([[[296.5076 , 296.37006, ..., 295.75443, 295.6565 ],
        [296.3718 , 295.8929 , ..., 295.39032, 294.9773 ],
        ...,
        [295.55945, 295.2863 , ..., 294.49417, 293.64276],
        [295.90247, 295.60196, ..., 294.41473, 292.2122 ]],

       [[296.21115, 295.9582 , ..., 295.90155, 295.75537],
        [296.14667, 295.56586, ..., 295.54822, 295.10535],
        ...,
        [295.1618 , 294.9685 , ..., 294.12482, 293.33273],
        [295.4936 , 295.26575, ..., 293.999  , 291.8243 ]],

       ...,

To work in the next step of my process, e.g. the concatenating of the other months/years, and later to analyse the data, the structure should be (I think):

dimensions: lat, lon, and time, and
data variable: temperature (or precipitation)

time would start at 01/01/1950 and goes until 31/12/2005

Improve this question
2
  • I think my problem is that I don't know how to change dimensions to variables, and change the type from xarray.DataArrayband to xarray.DataArray, but I am not sure, that's why there is no more code.
    – Gijs van den Dool
    Jan 17 at 12:52
  • the time dimension can be made with: <br> start_day = 2005-01-01 (from file name) <br> band_count = len(test.band) <br> time = pd.date_range(start_day, freq="D", periods=band_count) <br>
    – Gijs van den Dool
    Jan 17 at 13:55

1 Answer 1

Reset to default
1

I will assume that your GEE exports have geolocation information on them. Basically you can read things in, set the dates axis and then concatenate together. I am assuming here annual monthly datasets (so a GeoTIFF for each year with 12 bands), but you can modify that to suit your data format:

from osgeo import gdal
import numpy as np
import rioxarray as rio
import xarray as xr
import pandas as pd

x = []
for year in [2020, 2021]: # for example
    # Create a sample geotiff file 
    a = np.random.rand(12, 360, 720)
    dst = gdal.GetDriverByName("GTiff").Create("sample_file.tif", xsize=720, ysize=360,
                        bands=12, eType=gdal.GDT_Float32)
    dst.WriteArray(a)
    # Assume geolocation etc are properly set
    dst = None

    ds = rio.open_rasterio("sample_file.tif")
    # This opens the GeoTIFF as with x, y, band coordinates
    # Rename "band" to "time"
    ds = ds.rename({"band":"time"})
    time_idx = pd.to_datetime([f"{year}-{i:02d}" for i in range(1, 13)])

    ds = ds.assign_coords({"time":("time", time_idx)})
    x.append(ds)
    
# Stack the individual years along "time" dimension
ds = xr.concat(x, dim="time")

In my case, the final xarray dataset looks like:

<xarray.DataArray (time: 24, y: 360, x: 720)>
array([[[0.42069325, 0.17924134, 0.7396588 , ..., 0.7511819 ,
[...]
Coordinates:
  * time         (time) datetime64[ns] 2020-01-01 2020-02-01 ... 2021-12-01
  * x            (x) float64 0.5 1.5 2.5 3.5 4.5 ... 716.5 717.5 718.5 719.5
  * y            (y) float64 0.5 1.5 2.5 3.5 4.5 ... 356.5 357.5 358.5 359.5
    spatial_ref  int64 0
Attributes:
    scale_factor:  1.0
    add_offset:    0.0
Improve this answer
5
  • Thanks Jose, ds = ds.assign_coords({"time":("time", time_idx)}) is something I can add and ds = ds.rename({"band":"time"}) is also something I can do but wouldn't it be better, when concatenating, that the bands are stored as a value array instead of a dimension?
    – Gijs van den Dool
    Jan 17 at 14:02
  • It wasn't clear what your original input data looked like. I assumed one variable per file with N timesteps as individual bands. Do you mean that each geotiff has N bands that refer to N different magnitudes (temp, precip, etc) for a given timestep? In that case, you may need to check the band_as_variable option and modify the code above slightly (this requires a fairly new rioxarray)
    – Jose
    Jan 17 at 14:48
  • No, your assumption was correct; the tiff is parameter-specific (temp, precip) and has daily values, but the next process is looking at netCDF files where the parameter is defined as a variable.
    – Gijs van den Dool
    Jan 17 at 14:53
  • I had a look, and I don't think band_as_variable will work, as it will load bands in a raster to separate variables, and (again perhaps due to my lack of knowledge) I was looking for a way to name the block as "temperature" or "precipitation" so that the other step knows what kind of data is in the block.
    – Gijs van den Dool
    Jan 17 at 15:01
  • import xarray as xr x = [] x.append(test) x.append(test2) test3 = xr.concat(x, dim="time") this works, so I can get one file for the whole period.
    – Gijs van den Dool
    Jan 17 at 16:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.