Why are the coordinates in the ODC different from the input data?

Question

Why are the coordinates in the ODC different from the input data?

I prepared the following example for clarification and discussion. The used dataset is shown in the following two images:

• Simple plot of the data using matplotlib:
• The same data rendered by panoply:
• the data looks like the following:

  height = xarray.DataArray([[1., 2., 3.], [4., 5., 6.]], dims=(
      "latitude", "longitude"), coords={"latitude": [1., 2.], "longitude": [3., 4., 5.]})
  ds = xarray.Dataset({"height": height})
  ds.rio.write_crs(4326, inplace=True)
  

  resulting in the following output:

  Dimensions:      (latitude: 2, longitude: 3)
  Coordinates:
    * latitude     (latitude) float64 1.0 2.0
    * longitude    (longitude) float64 3.0 4.0 5.0
      spatial_ref  int64 ...
  Data variables:
      height       (latitude, longitude) float64 ...
  Attributes:
      grid_mapping:  spatial_ref
  
  Additional information
  height values:
  [[1. 2. 3.]
   [4. 5. 6.]]
  spatial_ref: EPSG:4326
  

This data is created using xarray and rioxarray and written to disk as netcdf. The netcdf is then registered in a datacube and loaded.

I expected to get the same output in terms of data (values and coordinates), but the datacube is using shifted coordinates:

<xarray.Dataset>
Dimensions:      (latitude: 2, longitude: 3, time: 1)
Coordinates:
  * time         (time) datetime64[ns] 1970-01-01
  * latitude     (latitude) float64 0.5 1.5
  * longitude    (longitude) float64 2.5 3.5 4.5
    spatial_ref  int32 4326
Data variables:
    height       (time, latitude, longitude) float64 1.0 2.0 3.0 4.0 5.0 6.0
Attributes:
    crs:           EPSG:4326
    grid_mapping:  spatial_ref

In addition, the coordinates change depending on the submitted query, but the data variables are the same:

Query 1 : {'latitude': (0, 2), 'longitude': (2, 5)}
Coords 1: Coordinates:
  * time         (time) datetime64[ns] 1970-01-01
  * latitude     (latitude) float64 0.5 1.5
  * longitude    (longitude) float64 2.5 3.5 4.5
    spatial_ref  int32 4326
Values 1: [[[1. 2. 3.]
  [4. 5. 6.]]]

Query 2 : {'latitude': (1, 3), 'longitude': (3, 6)}
Coords 2: Coordinates:
  * time         (time) datetime64[ns] 1970-01-01
  * latitude     (latitude) float64 1.5 2.5
  * longitude    (longitude) float64 3.5 4.5 5.5
    spatial_ref  int32 4326
Values 2: [[[1. 2. 3.]
  [4. 5. 6.]]]

I could not find anything regarding this topic in the documentation, the issues, and gis.stackexchange.

Do you need any additional information?

The script to create these outputs is published as gist. It requires a running datacube environment, psql installed, and the delete_odc_product.sql file. The script performs the following steps:

1. Create some files:
   1. netcdf file
   2. datacube product specification
   3. datacube dataset specification
   4. datacube.conf
2. Init datacube database if not done already
3. Add the product if not done already
4. Add the dataset if not done already
5. Loads the dataset
6. Print metadata and data
7. Remove product from datacube using psql and the above mentioned delete_odc_product.sql
8. Removes all created files

Versions

• os: Ubuntu 20.04.2 LTS
• python: 3.8.5
• datacube: Open Data Cube core, version 1.8.3

Answer 1

This issue can be solved by using the align parameter (...of datacube.Datacube.load or via the product definition .yaml using the load:align element).

If you want to have the coordinates beeing in the center of the pixel, align must be pixel-resolution/2.

Code snippet:

query = {"align": (0.5, 0.5)}
data = dc.load('dim_example', **query)
print("""Query   : {0}
Coords  : {1}
Values  : {2}
Min     : {3}
Mean    : {4}
Max     : {5}
""".format(query, data.coords, data.height.values,
           data.height.min(), data.height.mean(), data.height.max()))

Output:

Query   : {'align': (0.5, 0.5)}
Coords  : Coordinates:
  * time         (time) datetime64[ns] 1970-01-01
  * latitude     (latitude) float64 1.0 2.0
  * longitude    (longitude) float64 3.0 4.0 5.0
    spatial_ref  int32 4326
Values  : [[[1. 2. 3.]
  [4. 5. 6.]]]
Min     : <xarray.DataArray 'height' ()>
array(1.)
Coordinates:
    spatial_ref  int32 4326
Mean    : <xarray.DataArray 'height' ()>
array(3.5)
Coordinates:
    spatial_ref  int32 4326
Max     : <xarray.DataArray 'height' ()>
array(6.)
Coordinates:
    spatial_ref  int32 4326

See the documentation of datacube.Datacube.load or product definition for more details.