mlubej
commented
2 years ago I've been playing around with this and trying to make it work in openEO via the VITO backend.
My starting point was a geodataframe of 10 parcels (to start simple) which are relatively close together:
This is the code used to get to the final datacube:
# set extent from the geodataframe
minx, miny, maxx, maxy = gdf.unary_union.envelope.bounds
spat_ext = dict(west=minx, east=maxx, north=maxy, south=miny, crs=gdf.crs.to_epsg())
temp_ext = ["2019-01-01", "2019-12-31"]
# link to collection, specify bands
s2 = connection.load_collection(
'SENTINEL2_L2A_SENTINELHUB',
spatial_extent=spat_ext,
temporal_extent=temp_ext,
bands=["B04","B08","SCL"]
)
# filter to specific geometries
s2 = s2.filter_spatial(geometries = json.loads(gdf.geometry.to_json()))
# SCL cleanup/filtering
s2 = s2.process("mask_scl_dilation", data=s2, scl_band_name="SCL")
# calculate NDVI
s2 = s2.ndvi(red="B04", nir="B08", target_band='NDVI')
# download the data and load it
job = s2.execute_batch("cube.nc", out_format="netCDF")
ds = xarray.load_dataset('./cube.nc')Observations:
- If I understand correctly, the above process filters the data after it downloads it, which is not optimal, since it would be better to filter already only the area of interest while downloading.
- for 2/10/100 parcels it takes 5/6/19 min to execute the datacube pipeline for a single year (2019)
- there seems to be some initialization period regardless of the amount of data you're downloading
- it seems to me this is not so optimal and makes me question a bit if I'm doing things right, because this might prove to be a bottleneck when downloading for a large number of parcels (order of 1M) and spanning a larger spatial extent
- should any parallelization be implemented on the user side?
- didn't work for any CRS other than 4326 (WGS 84)
The steps from this point (having a .nc dataset on disk) start to differ now w.r.t. to the approach we would be taking:
- pixel based
- parcel based
Pixel-based approach
The goal here would be to obtain a pandas dataframe, where each row represents a timeseries for a single pixel, and where only the pixels from the area of interest are taken, but one should still be able to subset the dataset to a specific geometry - spatial context should be available.
I'm not sure what the proper way of doing this is, but my initial approach took me down the rioxarray path, where the idea was to first spatially clip the xarray to the extent of the geometry of interest and then convert the xarray to a dataframe.
Another way would be to rasterize the geometry mask into the xarray and then convert it, but I guess you would need to add a timeless feature to the xarray, haven't tried that yet.
Parcel-based approach
This one is a bit more manageable, because you can just add the spatial aggregation function to the datacube pipeline.
However, this requires re-downloading the data, where the spatial aggregation is done on the fly. It would make much more sense to just run the aggregation on the pixel-level dataframe if it is already available.
If the pixel-level dataset is not available, then this is fine and can be achieved via:
from openeo.rest.conversions import timeseries_json_to_pandas
aggregated_data_json = s2.aggregate_spatial(geometries=json.loads(gdf.geometry.to_json()), reducer='mean').execute()
timeseries_json_to_pandas(aggregated_data_json)This returns an output in the following format:
Which means that we would need to reshape the dataframe a bit to get the row-based format which is easier to work with.
thoughts/comments/suggestions?
jdries
Inputs: