Partial data loading from the Open Data Cube is not working

Question

I have deployed the cube in a box docker in our local environment and managed to index a couple of Sentinel-2 level 2 data. However, when I tried to load my AOI from those datasets I receive an empty list. The code used is copied below:

# Create a query object
lat, lon = 22.821, 28.518
buffer = 0.05

query = {
    'time': ('2020-12-11', '2020-12-13'),
    'x': (lon - buffer, lon + buffer),
    'y': (lat + buffer, lat - buffer),
    'output_crs': 'epsg:3857',
    'resolution':(-10,10),
}

# Load an image, I am loading a Senitnel2 Level1 data from the datacube
bands = ['B02_10m', 'B03_10m', 'B04_10m']

ds = dc.load(product='s2a_sen2cor_granule',
             measurements=bands,
             dask_chunks={},
             **query)

print(ds)

If I ignore the longitude and latitude then it works okay

# Create a query object
lat, lon = 22.821, 28.518
buffer = 0.05

query = {
    'time': ('2020-12-11', '2020-12-13'),
    'output_crs': 'epsg:3857',
    'resolution':(-10,10),
}

#load data
ds = dc.load(product='s2_l2a',
             dask_chunks={},
             **query)

print(ds)

And the output I receive is as follows:

<xarray.Dataset>
Dimensions:      (time: 2, x: 484, y: 591)
Coordinates:
  * time         (time) datetime64[ns] 2021-07-10T08:53:55 2021-07-10T08:54:10
  * y            (y) float64 2.814e+06 2.814e+06 ... 2.803e+06 2.803e+06
  * x            (x) float64 2.765e+06 2.765e+06 ... 2.775e+06 2.775e+06
    spatial_ref  int32 6933
Data variables:
    B01          (time, y, x) uint16 dask.array<chunksize=(1, 591, 484), meta=np.ndarray>
    B02          (time, y, x) uint16 dask.array<chunksize=(1, 591, 484), meta=np.ndarray>
    B03          (time, y, x) uint16 dask.array<chunksize=(1, 591, 484), meta=np.ndarray>
    B04          (time, y, x) uint16 dask.array<chunksize=(1, 591, 484), meta=np.ndarray>
    B05          (time, y, x) uint16 dask.array<chunksize=(1, 591, 484), meta=np.ndarray>
    B06          (time, y, x) uint16 dask.array<chunksize=(1, 591, 484), meta=np.ndarray>
    B07          (time, y, x) uint16 dask.array<chunksize=(1, 591, 484), meta=np.ndarray>
    B08          (time, y, x) uint16 dask.array<chunksize=(1, 591, 484), meta=np.ndarray>
    B8A          (time, y, x) uint16 dask.array<chunksize=(1, 591, 484), meta=np.ndarray>
    B09          (time, y, x) uint16 dask.array<chunksize=(1, 591, 484), meta=np.ndarray>
    B11          (time, y, x) uint16 dask.array<chunksize=(1, 591, 484), meta=np.ndarray>
    B12          (time, y, x) uint16 dask.array<chunksize=(1, 591, 484), meta=np.ndarray>
    SCL          (time, y, x) uint8 dask.array<chunksize=(1, 591, 484), meta=np.ndarray>
    AOT          (time, y, x) uint16 dask.array<chunksize=(1, 591, 484), meta=np.ndarray>
    WVP          (time, y, x) uint16 dask.array<chunksize=(1, 591, 484), meta=np.ndarray>
Attributes:
    crs:           epsg:6933
    grid_mapping:  spatial_ref

What could be the possible reason?

UPDATE for @Robbi Bishop-Taylor: Hi, I tried the query you posted below, but no joy. So, I cleared the ODC DB and started from the beginning with a new dataset, please see the image below on google map and the red mark is the coordinate I used in the query

I can also confirm that the data is indexed correctly on datacube-explorer

The new query parameters are as follows:

lat, lon = 22.821, 28.518
buffer = 0.05

query = {
    'time': ('2021-07-24', '2021-07-27'),
    'x': (lon - buffer, lon + buffer),
    'y': (lat + buffer, lat - buffer),
    'output_crs': 'epsg:6933',
    'resolution':(-10,10),
}

Still didn't manage to load AOI/partial data.

NOTE: Sentinel-2 data I have used is S2A_MSIL2A_20210725T083601_N0301_R064_T35QPF_20210725T115425.SAFE

Answer 1

As far as I can tell, your query is set up correctly - it's entirely possible to query by geographic (lat-lon) coords and return data in a different coordinate reference system (e.g. EPSG:3857 with metre units like you do above, although @whatnick is correct that this CRS isn't usually used for analysis). Without specifying any spatial extents you simply get back all data regardless of location - in this case it's loaded in the data's native CRS which is EPSG:6933 - this is a global grid (WGS 84 / NSIDC EASE-Grid 2.0 Global) so I think that should be OK too.

I've plotted out your input geographic coordinates with their 0.05 buffer, and the bounds of the data that you loaded into EPSG:6933 without any spatial query. This looks like your indexed data only covers a small area that does not overlap with your query bounds:

To rule this out, I'd recommend re-running your code using the following lat-lon query input instead:

# Create a query object
lat, lon = 22.58318610, 28.70269922
buffer = 0.5

query = {
    'time': ('2020-12-11', '2020-12-13'),
    'x': (lon - buffer, lon + buffer),
    'y': (lat + buffer, lat - buffer),
    'output_crs': 'epsg:3857',
    'resolution':(-10, 10),
}

Answer 2

I managed to fix this issue. The query was absolutely fine; however, the preparation script I have used from the data cube Github repository needed a tweak. FYI, I used the sen2cor_prepare.py to index the L2A data downloaded from the Copernicus hub, this script generated the YAML file with the wrong coordinates; I just added the following line in the get_coords method to make it work. Thanks to my friend who helped me to sort this out.

spatial_ref.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)

NOTE: It is due to the version of GDAL I am using, may only work with version >3

Answer 3

It looks like your data is stored in 'epsg:6933' (see the attributes on the dataset you loaded without partial data loading), which is an Albers equal area projection for Africa (coordinates are in metres), but you're trying to use lat/lon coordinates to clip and load the data. You best bet is to work out what's going wrong with your coordinate systems. Also, the measurements you are passing to dc.load (B04_10m etc) do not match those in the data (B04.), and your product names don't match: 's2_l2a' vs 's2a_sen2cor_granule'.

Answer 4

There are all sorts of things going on with this particular set of variables. EPSG:3857 uses a co-ordinate system in meters rather than degrees. Also it is only really used in slippy maps, for analysis an equal area projection would be better.

Such as Albers Conic tuned to your area : https://en.m.wikipedia.org/wiki/Albers_projection