From 7809154f7b3802500a4d69c50438f5d9f1d100e0 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Thu, 28 May 2026 08:27:20 -0700 Subject: [PATCH 1/2] resample: refresh transform to match actual coord orientation (#2571) The transform refresh block at the end of resample() was hard-coding a north-up rasterio tuple (positive x scale, negative y scale, origin at the top-left edge). For arrays with ascending y or descending x coords that produced an attrs['transform'] that did not match the coords on the same DataArray, so downstream code that georeferenced via attrs['transform'] would get a flipped raster. Use the signed px / py returned by _new_coords (positive when the coord ascends along the axis, negative when it descends) and anchor the origin at *_edge_start, which is already the leading edge of the first row / column on the side of vals[0]. That matches rasterio's `(col=0, row=0) -> first array pixel` convention for every coord orientation. Tests cover all four orientation combinations (x asc/desc x y asc/desc) and assert the refreshed transform round-trips against the output coords. --- xrspatial/resample.py | 21 +++++++------- xrspatial/tests/test_resample.py | 48 ++++++++++++++++++++++++++++++++ 2 files changed, 58 insertions(+), 11 deletions(-) diff --git a/xrspatial/resample.py b/xrspatial/resample.py index c316590a3..9e5d5eacf 100644 --- a/xrspatial/resample.py +++ b/xrspatial/resample.py @@ -1361,19 +1361,18 @@ def _new_coords(vals, n_out): if has_nodata: new_attrs['_FillValue'] = float('nan') - # Refresh `transform` if the input had one. The rasterio 6-tuple is - # (res_x, 0.0, left, 0.0, -res_y, top). `top` is the upper edge of - # the first row, which is `y_edge_start` when y is descending and - # `y_edge_end` when y is ascending. `left` is the lower edge of the - # first column, which is `x_edge_start` when x is ascending and - # `x_edge_end` when x is descending. + # Refresh `transform` if the input had one. The rasterio 6-tuple + # `(a, b, c, d, e, f)` maps `(col, row) -> (x, y)` for the first + # array pixel at `(col=0, row=0)`, so the scale signs and the + # origin corner have to follow the actual array layout rather + # than assuming a north-up grid. `px` / `py` from `_new_coords` + # are already signed (positive when the coord ascends along the + # axis, negative when it descends), and `*_edge_start` is the + # leading edge of the first row / column on the side of + # `vals[0]` -- exactly what rasterio wants for `c` and `f`. if 'transform' in agg.attrs: - out_res_x = abs(px) - out_res_y = abs(py) - top = y_edge_start if y_vals[0] > y_vals[-1] else y_edge_end - left = x_edge_start if x_vals[0] < x_vals[-1] else x_edge_end new_attrs['transform'] = ( - out_res_x, 0.0, left, 0.0, -out_res_y, top, + px, 0.0, x_edge_start, 0.0, py, y_edge_start, ) # Resample replaces sentinel pixels with NaN regardless of input diff --git a/xrspatial/tests/test_resample.py b/xrspatial/tests/test_resample.py index 09b434bd7..d357f6022 100644 --- a/xrspatial/tests/test_resample.py +++ b/xrspatial/tests/test_resample.py @@ -151,6 +151,54 @@ def test_transform_refreshed_on_downsample(self): 2.0, 0.0, 100.0, 0.0, -2.0, 200.0, ) + @staticmethod + def _raster_with_orientation(x_asc, y_asc): + """4x4 raster with the chosen coord orientation and a matching + rasterio transform for `(col=0, row=0) -> upper-left edge of + pixel [0, 0]`. Issue #2571. + """ + res = 1.0 + x_step = res if x_asc else -res + y_step = res if y_asc else -res + # Pick first-pixel center so the leading edge sits on a round + # number; the actual value doesn't matter for the round-trip. + x0_center = 0.5 if x_asc else 3.5 + y0_center = 0.5 if y_asc else 3.5 + x = np.array([x0_center + i * x_step for i in range(4)]) + y = np.array([y0_center + i * y_step for i in range(4)]) + left = x0_center - x_step / 2 + top = y0_center - y_step / 2 + transform = (x_step, 0.0, left, 0.0, y_step, top) + data = np.arange(16, dtype=np.float32).reshape(4, 4) + return xr.DataArray( + data, + dims=['y', 'x'], + coords={'y': y, 'x': x}, + attrs={'transform': transform, 'res': (res, res)}, + ) + + @pytest.mark.parametrize('x_asc', [True, False]) + @pytest.mark.parametrize('y_asc', [True, False]) + def test_transform_matches_coord_orientation(self, x_asc, y_asc): + """The refreshed transform must round-trip against the output + coords for every coord orientation (issue #2571).""" + raster = self._raster_with_orientation(x_asc, y_asc) + out = resample(raster, scale_factor=0.5) + + a, b, c, d, e, f = out.attrs['transform'] + # Off-diagonals are always zero for an axis-aligned grid. + assert b == 0.0 and d == 0.0 + # Scale signs follow the coord direction. + assert (a > 0) == x_asc + assert (e > 0) == y_asc + # transform * (col + 0.5, row + 0.5) -> pixel-center coord. + out_x = out[out.dims[-1]].values + out_y = out[out.dims[-2]].values + for col, expected_x in enumerate(out_x): + assert pytest.approx(c + a * (col + 0.5)) == expected_x + for row, expected_y in enumerate(out_y): + assert pytest.approx(f + e * (row + 0.5)) == expected_y + def test_transform_absent_stays_absent(self): # grid_4x4 fixture has no transform attr. data = np.arange(16, dtype=np.float32).reshape(4, 4) From 34dc11002a022f65712e6abb67a96beee7a73a78 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Thu, 28 May 2026 08:29:18 -0700 Subject: [PATCH 2/2] Address review nits: clarify test docstring and split assert (#2571) - `_raster_with_orientation` docstring now says "leading-edge corner" rather than "upper-left edge", which only matches when both coords are north-up. - Split `assert b == 0.0 and d == 0.0` into two asserts so a failure identifies which off-diagonal drifted. --- xrspatial/tests/test_resample.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/xrspatial/tests/test_resample.py b/xrspatial/tests/test_resample.py index d357f6022..f61c208b0 100644 --- a/xrspatial/tests/test_resample.py +++ b/xrspatial/tests/test_resample.py @@ -154,8 +154,10 @@ def test_transform_refreshed_on_downsample(self): @staticmethod def _raster_with_orientation(x_asc, y_asc): """4x4 raster with the chosen coord orientation and a matching - rasterio transform for `(col=0, row=0) -> upper-left edge of - pixel [0, 0]`. Issue #2571. + rasterio transform for `(col=0, row=0) -> leading-edge corner + of pixel [0, 0]` (the geographic "upper-left" only when both + coords are in the conventional north-up direction). Issue + #2571. """ res = 1.0 x_step = res if x_asc else -res @@ -187,7 +189,8 @@ def test_transform_matches_coord_orientation(self, x_asc, y_asc): a, b, c, d, e, f = out.attrs['transform'] # Off-diagonals are always zero for an axis-aligned grid. - assert b == 0.0 and d == 0.0 + assert b == 0.0 + assert d == 0.0 # Scale signs follow the coord direction. assert (a > 0) == x_asc assert (e > 0) == y_asc