Reprojection produces a GeoBox that does not match when destination CRS is EPSG 4326

[fbunt]

I ran into this issue while reprojecting some data but I think the underlying cause is that GeoBox objects created by different means for the same grid location will have affine transforms that are slightly different. It may be specific to EPSG 4326 since I haven't run into this issue with other projections.

When reprojecting to EPSG 4326, the geobox for the result does not match what I input. A short example:

import numpy as np
import xarray as xr
import rioxarray as riox

src = xr.DataArray(
    np.ones((1, 10, 10)),
    dims=("band", "y", "x"),
    coords=([1], np.arange(10), np.arange(10)[::-1]),
).rio.write_crs("EPSG:3310")

dst_gb = src.odc.geobox.to_crs(4326)
dst = src.odc.reproject(dst_gb)
print("GeoBox comparison:", dst.odc.geobox == dst_gb)
print("Shape comparison:", dst.odc.geobox.shape == dst_gb.shape)
print("CRS comparison:", dst.odc.geobox.crs == dst_gb.crs)
print("Affine comparison:", dst.odc.geobox.affine == dst_gb.affine)
print()
print("dst_gb affine")
print(repr(dst_gb.affine))
print("result affine")
print(repr(dst.odc.geobox.affine))
print()
print("Affine comparison with tolerance:", np.allclose(list(dst_gb.affine), list(dst.odc.geobox.affine)))

GeoBox comparison: False
Shape comparison: True
CRS comparison: True
Affine comparison: False

dst_gb affine
Affine(1.0200436314046532e-05, 0.0, -120.00002388788778,
       0.0, -1.0200436314046532e-05, 38.01647531889047)
result affine
Affine(1.0200435637076365e-05, 0.0, -120.00001592863761,
       0.0, -1.0200435637131023e-05, 38.01647279736898)

Affine comparison with tolerance: True

I'm using version 0.4.2.

[Kirill888]

@fbunt that's a limitation of Pixel->World transform representation we use on Xarray objects, and is not .reproject specific. To report correct geobox after slicing we have to recompute it from coordinates, and not just read from some attribute, actual "affine" matrix is not stored on the xarray object, only pixel center coordinates are stored.

Example of the same issue without reproject:

from odc.geo.geobox import GeoBox
from odc.geo.math import snap_affine
from odc.geo.xr import xr_zeros

src_gb = GeoBox.from_bbox((0.5, 0.5, 10.5, 10.5), shape=(10, 10), crs=3310, tight=True)
dst_gb = src_gb.to_crs(4326)

dst = xr_zeros(dst_gb)

print("GeoBox comparison:", dst.odc.geobox == dst_gb)
print("Shape comparison:", dst.odc.geobox.shape == dst_gb.shape)
print("CRS comparison:", dst.odc.geobox.crs == dst_gb.crs)
print("Affine comparison:", dst.odc.geobox.affine == dst_gb.affine)
print("Area comparison:", (dst.odc.geobox.extent ^ dst_gb.extent).area == 0.0)
print("\ndst_gb affine")
print(repr(dst_gb.affine))
print("\nresult affine")
print(repr(dst.odc.geobox.affine))

# Difference in Affine
A = dst_gb.affine  # < desired
A_ = dst.odc.affine  # < actual, reconstructed from coords
dA = (~A) * A_  # pixel relation Actual -> Requested

print("\ndiff affine, pixel domain")
print(repr(dA))
print("\ndiff affine, rounded")
display(snap_affine(dA, tol=1e-6))

display(src_gb, dst_gb)

out:

GeoBox comparison: False
Shape comparison: True
CRS comparison: True
Affine comparison: False
Area comparison: False

dst_gb affine
Affine(1.020043632879809e-05, 0.0, -120.00000366055548,
       0.0, -1.020043632879809e-05, 38.01648557430506)

result affine
Affine(1.0200436329035676e-05, 0.0, -120.00000366055548,
       0.0, -1.0200436328489104e-05, 38.01648557430505)

diff affine, pixel domain
Affine(1.0000000000232918, 0.0, 0.0,
       0.0, 0.9999999999697087, 4.656612873077393e-10)

diff affine, rounded

Affine(1.0, 0.0, 0.0,
       0.0, 1.0, 0.0)

1-meter pixel resolution from your example converted to degrees is just too small for float representation to handle without any losses. This issue is more common when working with geographic coordinates because of floating point precision.

[Kirill888]

We need an FAQ we can point people to when confusion about pixel coords comes up,
related

• Task: add documentation for pixel grids odc-stac#95
• Issue with cropped xarray data causing xarray's .where() to behave differently #115

[fbunt]

Ah ok. That makes sense. Thanks! I'm happy to close this if that works for you.