OpenEO sometimes returns a datacube in a different projection than requested

openEO Platform
1

Hi all,

I am retrieving data with load_collection, using the following code:

s2_cube = connection.load_collection(
                    "SENTINEL2_L2A_SENTINELHUB",
                    spatial_extent={
                        "west": cfg.bounds[0],
                        "south": cfg.bounds[1],
                        "east": cfg.bounds[2],
                        "north": cfg.bounds[3],
                        "crs": cfg.crs,
                    },

I use UTM zone EPSG code as the crs, and utm coordinates. Most of the time s2_cube has the exact requested crs and bounds, but sometimes it gets the crs and coordinates of the neighboring UTM zone.

For instance, with the following UTM 50N request:

    "aoi": [
        750690.0,
        5372910.0,
        760590.0,
        5382810.0
    ],
    "crs": "EPSG:32650",

The retrieved cube actually is in UTM 51N, with coordinates ranging from 307125. to 317775. in x and from 5369865. to 5380505. in y:

>>> ds=xr.open_dataset("sentinel2.nc")
>>> ds.crs
<xarray.DataArray 'crs' ()> Size: 1B
[1 values with dtype=|S1]
Attributes:
    crs_wkt:      PROJCS["WGS 84 / UTM zone 51N", GEOGCS["WGS 84", DATUM["Wor...
    spatial_ref:  PROJCS["WGS 84 / UTM zone 51N", GEOGCS["WGS 84", DATUM["Wor...
>>> ds.x.min()
<xarray.DataArray 'x' ()> Size: 8B
array(307125.)
>>> ds.x.max()
<xarray.DataArray 'x' ()> Size: 8B
array(317775.)
>>> ds.y.min()
<xarray.DataArray 'y' ()> Size: 8B
array(5369865.)
>>> ds.y.max()
<xarray.DataArray 'y' ()> Size: 8B
array(5380505.)

The expected number of pixels is also different: the request has (760590.0-750690.0)/10.=990 pixels while we actually retrieve (317775.-307125.)/10.=1065 pixels in the datacube.

This is problematic for precise AOI retrieval. Is there a way to force openeo to produce datacube using the requested CRS and bounds ?

Regards,

Julien

2

in the openEO backend we try to automatically pick the best UTM for your area of interest

and these coordinates seems to fall in the UTM 51N zone even though you specified them as UTM 50N

I think it should be possible to force your desired projection with something like

s2_cube = s2_cube.resample_spatial(projection="EPSG:32650")
3

Actually, it falls in the overlap between both tiles, though maybe closer to the center of 51N.

image
image811Ă—477 111 KB

I can reproject with resample_spatial, but it will incur blur. I wanted raw Sentinel-2 pixels from the 50UQU tile, without any resampling.

4

Hi Julien,

OpenEO can use raw data that is already in the target CRS in this case, avoiding resampling artifacts. You can check the used products in the “derived_from” metadata. I ran a test with your extent, and it matched the raw data perfectly.

Additionally, you can specify the resolution in the resample to be sure all stays fine
s2_cube = s2_cube.resample_spatial(resolution=10.0, projection="EPSG:32650")

Emile

5

Great! Just tested it, works like a charm. I get exactly the requested crs and aoi.

Thanks a lot to you @emile.sonneveld and @stefaan.lippens

1 Like