corteva / rioxarray

geospatial xarray extension powered by rasterio
https://corteva.github.io/rioxarray
Other
502 stars 80 forks source link

`rio.transform()` does not retrieve exact transform stored in `rio.write_transform()` #739

Closed DahnJ closed 4 months ago

DahnJ commented 5 months ago

Code Sample, a copy-pastable example if possible

import xarray as xr
import rioxarray
from affine import Affine
lat = [2.5, 4.5]
lon = [2.5, 4.5]

da = xr.DataArray(
    [[1., 0], [-1, 2]],
    coords = [lat, lon],
    dims = ["latitude", "longitude"]
)

da.rio.write_crs("epsg:4326", inplace=True)
da.rio.set_spatial_dims(x_dim="longitude", y_dim="latitude", inplace=True)
da = da.rio.write_transform(Affine(5, 0, 100.5, 0, -10, 1.5), inplace=True)

# Stored as expected
da[da.rio.grid_mapping].attrs['GeoTransform']
# '100.5 5.0 0.0 1.5 0.0 -10.0'

# Transform returned with different translation part
da.rio.transform()

# Affine(5.0, 0.0, 0.0,
#       0.0, -10.0, 7.5)

Problem description

I would expect rio.transform(recalc=False) to return the exact transform stored using write_transform.

However, this is only the case of the linear part of the transformation – the offset/translation part is recalculated using rio._unordered_bounds() (source) which takes them from rio._internal_bounds() (source)

This may not be a bug, but it's a surprising behaviour I decided to note.

Expected Output

# Affine(5.0, 0.0, 100.5,
#       0.0, -10.0, 1.5)

Environment Information

rioxarray (0.15.1) deps:
  rasterio: 1.3.9
    xarray: 2024.1.1
      GDAL: 3.8.3
      GEOS: 3.12.1
      PROJ: 9.3.1

Other python deps:
     scipy: 1.12.0
    pyproj: 3.6.1

System:
    python: 3.11.7 | packaged by conda-forge | (main, Dec 23 2023, 14:38:07) [Clang 16.0.6 ]

Installation method

Conda

snowman2 commented 4 months ago

However, this is only the case of the linear part of the transformation – the offset/translation part is recalculated using rio._unordered_bounds() (source) which takes them from rio._internal_bounds() (source)

With xarray you can slice/select subsets of data and have the same GeoTransform stored in the metadata. The GeoTransform would be incorrect in this scenario. Due to this, the resolution component of the GeoTransform is the only part that can be trusted.

DahnJ commented 4 months ago

Thanks for the explanation @snowman2, this makes complete sense.