opendatacube / odc-geo

GeoBox and geometry utilities extracted from datacube-core
https://odc-geo.readthedocs.io/en/latest/
Apache License 2.0
80 stars 12 forks source link

Strange resulting grid after xr_reproject #135

Closed Riccardo7-DS closed 6 months ago

Riccardo7-DS commented 6 months ago

After reprojecting my xarray dataset to a target geobox, the obtained dataset has the geobox (ds_repr.odc.geobox) of the target dataset and the same spatial resolution. However, some lat lon coordinates do not match by a very small value (e.g. when I do: ds_repr["lat"][1] - ds_target["lat"][1] I obtain: -1.7763568394002505e-15). This is a bit weird, do you have any ideas?

from odc.geo.xr import xr_reproject
import xarray as xr 

chunks={"time":-1, "lat":90, "lon":90}

ds = xr.open_zarr("/path/to/my/zarr"), 
                        chunks=chunks)
target_ds = xr.open_zarr("/path/to/my/zarr2"), 
                        chunks=chunks)

def odc_reproject(ds, target_ds, resampling:str):
    if target_ds.odc.geobox.crs.geographic:
        ds = ds.rename({"lon": "longitude", "lat": "latitude"})
    return xr_reproject(ds, target_ds.odc.geobox, 
                        resampling=resampling)

ds_repr = odc_reproject(ds, target_ds, resampling="bilinear")\
    .rename({"longitude": "lon", "latitude": `"lat"})
SpacemanPaul commented 6 months ago

Sounds like typical floating point errors to me.

E.g.

>>> (0.1 * 0.2) - 0.02
3.469446951953614e-18
Kirill888 commented 6 months ago

@Riccardo7-DS the error you are reporting is due to numeric precision. Your reference dataset was loaded from zarr, so it has whatever coordinates that were put in there by the process that created it. odc-geo then computes linear mapping from pixel coordinates to world coordinates (scale*(pix_idx + 0.5) + offset x2 for X/Y), and then reproject function recomputes coordinates given that linear mapping.

This problem is common when working with lon/lat coordinates #127, and especially when pixels snapping is not used. In this specific case there is nothing we can do at the library level, other than adding better docs and an FAQ.

Best I can offer is a work-around: after reproject you can copy coordinate values from your reference dataset into the warped one, use something like ds_repr.lon.data[:] = target_ds.lon.data.

Riccardo7-DS commented 6 months ago

Thank you guys for the clarification!