Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Strange resulting grid after xr_reproject #135

Closed
Riccardo7-DS opened this issue Mar 15, 2024 · 3 comments
Closed

Strange resulting grid after xr_reproject #135

Riccardo7-DS opened this issue Mar 15, 2024 · 3 comments
Labels
documentation Improvements or additions to documentation

Comments

@Riccardo7-DS
Copy link

Riccardo7-DS commented Mar 15, 2024

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
Copy link
Contributor

SpacemanPaul commented Mar 18, 2024

Sounds like typical floating point errors to me.

E.g.

>>> (0.1 * 0.2) - 0.02
3.469446951953614e-18

@Kirill888
Copy link
Member

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

@Kirill888 Kirill888 added the documentation Improvements or additions to documentation label Mar 18, 2024
@Riccardo7-DS
Copy link
Author

Thank you guys for the clarification!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
documentation Improvements or additions to documentation
Projects
None yet
Development

No branches or pull requests

3 participants