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

Updating coordinates on a dataset doesn't update geobox for dataarrays #157

Closed alexgleith closed 4 months ago

alexgleith commented 4 months ago

Updating coordiantes works at the dataset level, but isn't reflected in the arrays.

Reproduce with the below code using odc-geo==0.4.5

import xarray as xr
import odc.geo.xr  # noqa
from odc.geo.xr import xr_coords, assign_crs
from affine import Affine
from odc.geo.geobox import GeoBox

# File from https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20231106090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc
test_file = "/tmp/20240412090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc"

data = xr.open_dataset(test_file)
data = assign_crs(data, crs="EPSG:4326")

# Set up a new Affine and GeoBox
new_affine = Affine(0.01, 0.0, -180.0, 0.0, -0.01, 89.995, 0.0, 0.0, 1.0)
new_geobox = GeoBox(data.odc.geobox.shape, new_affine, data.odc.geobox.crs)

# First flip the dataset
data = data.reindex(lat=data.lat[::-1])

# Update the coordinates of the xarray to be precise
data = data.assign_coords(xr_coords(new_geobox))

# These should be the same
print(data.odc.geobox.affine[0])
print(data.analysed_sst.odc.geobox.affine[0])

0.01 0.010009765625

SpacemanPaul commented 4 months ago

You are calculating coordinates from an affine transform, then calculating the affine transform back from the coordinates. Why would you expect this round trip conversion to be exact?

alexgleith commented 4 months ago

Why would you expect this round trip conversion to be exact?

It is exact, and in previous versions, it applied to the dataarrays as well as the dataset. I'm saying that it's no longer working on the individual dataarrays.

To be more complete, the inferred coordinates from the NetCDF are messy (i.e., numbers like 0.010009765625). Providing a specific transform and computing coordinates leads to them being exact, and my intent is to write a COG with this exact affine transform defining the coordinates. This did work prior to 0.4.5 and I think it's changed in that release... haven't looked at the code though to find out where it might be.

SpacemanPaul commented 4 months ago

Ah OK, yes, NetCDF's raster representation is a poor match for xarray, odc-geo, GeoTiff, Zarr, and generally the rest of the ecosystem.

Kirill888 commented 4 months ago

@alexgleith does this work when coords are assigned per variable? are the coords on patched data variables as you expect? what's the meaning of reindex and why it's needed?

alexgleith commented 4 months ago

Oh, it's my mistake actually... not naming the x/y dims right 🤦