2

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 enter image description here

I can also confirm that the data is indexed correctly on datacube-explorer enter image description here

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

2

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:

enter image description here

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),
}
1
  • I have an update for you in my main post, thanks for your response.
    – prpp2021
    Aug 10 at 10:24
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

1

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'.

3
  • EPSG:6933 is a global equal area grid (WGS 84 / NSIDC EASE-Grid 2.0 Global), so I think that part should be all good: nsidc.org/data/ease Aug 10 at 0:24
  • @RobbiBishop-Taylor If the data is natively stored in albers (as the spatial_ref coordinate suggests), then is ODC smart enough to convert the lat/lon values to the Albers CRS to clip the datasets? Which is different from having them stored as lat/lon then outputting in a projected CRS...
    – cbur
    Aug 10 at 4:12
  • I really appreciate your response. The data is stored in EPSG:32635, not 6933. Clipping should be possible with any projection, as far I understood! The measurements are all okay, and it matches 's2a_sen2cor_granule'; this is why loading the whole dataset works perfectly; the only thing that doesn't work is the partial loading.
    – prpp2021
    Aug 10 at 10:40
-1

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

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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