.odc.geobox coordinate precision is not the same as DataArray coordinates

[saschahofmann]

There seems to be an issue with floating precisions when assigning the GeoBox to a xarray.DataArray related to floating point precision.

import numpy as np
import xarray as xr
import odc.geo.xr

lats = np.arange(21.95, 72.55, 0.05)
lons = np.arange(-44.55, 64.95, 0.05)
da = xr.DataArray(
    np.random.rand(lats.size, lons.size),
    coords=dict(latitude=("latitude", lats), longitude=("longitude", lons)),
)
(da.odc.geobox.coords["x"].values == da.longitude.values).sum()/da.longitude.size*100

Leads to a 0.41% agreement.

In my case, I am using a reference dataset to reproject another dataset but because of the imprecision the resulting reprojected dataset is not actually on the same grid as the reference. You can also see this when looking at the coordinates

The resulting coordinate resolutions are slightly off.

[Kirill888]

Yeah, that's expected, that's how IEEE floats are.

[saschahofmann]

Not sure I understand that argument. Is it because it has to work out the resolution from the distances between the coordinates. Even when I generate a GeoBox myself providing the resolution, I am not getting the same coordinates as np.arange would return?

[Kirill888]

@saschahofmann for starters instead of computing proportion of matches using exact comparison do something like max(abs(a-b)), and you'll see that the difference is something around 1e-14, maybe 1e-13

(da.odc.geobox.coords["x"].values - da.longitude.values).abs().max()

or something like that.

Related issues in ODC:

#115
#127
opendatacube/odc-stac#94
opendatacube/odc-stac#95

Xarray folks seem to accept that exact comparison for floating point coordinates is not always the right way, but there is no standard "Index" class that can handle those tiny deviations as of yet:

https://discourse.pangeo.io/t/example-which-highlights-the-limitations-of-netcdf-style-coordinates-for-large-geospatial-rasters/4140

you might need to dig through those discussions to find proposed Index solutions in xarray, last time I checked that was not standardized in xarray ecosystem yet.

[saschahofmann]

For now, I simply overwrite the coordinates with the original ones.

I think the biggest issue is that when I want to e.g. sum up the reference and regridded datarrays, xarray will produce an output of the coordinates that are exactly equal and drop the others. So it took me a little bit of debugging to find this.

Will check out the links you've sent. Thanks