From 65d3f2958d320848dbd01d6730f036cad1962596 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Tue, 14 Apr 2026 08:34:00 -0700 Subject: [PATCH] Fix resample interpolation coordinate mapping (#1202) The interpolation path used scipy.ndimage.zoom's edge-aligned formula which maps output pixel o to input pixel o*(N_in-1)/(N_out-1). This pinned the first output pixel to the first input pixel and the last to the last, creating a systematic half-pixel spatial shift vs. the block-centered output coordinate metadata. Switch all four backend paths (numpy, cupy, dask+numpy, dask+cupy) to use map_coordinates with the block-centered formula: input_pixel = (out_pixel + 0.5) * (N_in / N_out) - 0.5 This brings the interpolation path into alignment with both the aggregation path and the output coordinate labels. --- xrspatial/resample.py | 108 +++++++++++++++++++------------ xrspatial/tests/test_resample.py | 55 ++++++++++++++++ 2 files changed, 120 insertions(+), 43 deletions(-) diff --git a/xrspatial/resample.py b/xrspatial/resample.py index 3840e45d..21b69ba6 100644 --- a/xrspatial/resample.py +++ b/xrspatial/resample.py @@ -10,7 +10,6 @@ import numpy as np import xarray as xr from scipy.ndimage import map_coordinates as _scipy_map_coords -from scipy.ndimage import zoom as _scipy_zoom try: import dask.array as da @@ -19,10 +18,8 @@ try: import cupy - import cupyx.scipy.ndimage as _cupy_ndimage except ImportError: cupy = None - _cupy_ndimage = None from xrspatial.utils import ( ArrayTypeFunctionMapping, @@ -64,52 +61,88 @@ def _output_chunks(in_chunks, scale): for i in range(len(in_chunks))) -# -- NaN-aware zoom (NumPy) -------------------------------------------------- +# -- Block-centered coordinate mapping --------------------------------------- -def _nan_aware_zoom_np(data, zoom_yx, order): - """``scipy.ndimage.zoom`` with NaN-aware weighting. +def _block_centered_coords(n_in, n_out): + """Return input coordinates for each output pixel using block-centered mapping. + + Maps output pixel ``o`` to input pixel ``(o + 0.5) * (n_in / n_out) - 0.5``. + This places each output pixel at the center of its spatial footprint, + matching the convention used by ``_new_coords`` for output coordinate + metadata. + """ + o = np.arange(n_out, dtype=np.float64) + return (o + 0.5) * (n_in / n_out) - 0.5 + + +# -- NaN-aware interpolation (NumPy) ---------------------------------------- + +def _nan_aware_interp_np(data, out_h, out_w, order): + """Interpolate *data* to *(out_h, out_w)* with NaN-aware weighting. + + Uses ``scipy.ndimage.map_coordinates`` with block-centered coordinate + mapping so that sample positions match the output coordinate metadata. For *order* 0 (nearest-neighbour) NaN propagates naturally. For higher orders the zero-fill / weight-mask trick is used so that NaN pixels do not corrupt their neighbours. """ + iy = _block_centered_coords(data.shape[0], out_h) + ix = _block_centered_coords(data.shape[1], out_w) + yy, xx = np.meshgrid(iy, ix, indexing='ij') + coords = np.array([yy.ravel(), xx.ravel()]) + if order == 0: - return _scipy_zoom(data, zoom_yx, order=0, mode='nearest') + result = _scipy_map_coords(data, coords, order=0, mode='nearest') + return result.reshape(out_h, out_w) mask = np.isnan(data) if not mask.any(): - return _scipy_zoom(data, zoom_yx, order=order, mode='nearest') + result = _scipy_map_coords(data, coords, order=order, mode='nearest') + return result.reshape(out_h, out_w) filled = np.where(mask, 0.0, data) weights = (~mask).astype(data.dtype) - z_data = _scipy_zoom(filled, zoom_yx, order=order, mode='nearest') - z_wt = _scipy_zoom(weights, zoom_yx, order=order, mode='nearest') + z_data = _scipy_map_coords(filled, coords, order=order, mode='nearest') + z_wt = _scipy_map_coords(weights, coords, order=order, mode='nearest') + + result = np.where(z_wt > 0.01, + z_data / np.maximum(z_wt, 1e-10), + np.nan) + return result.reshape(out_h, out_w) - return np.where(z_wt > 0.01, - z_data / np.maximum(z_wt, 1e-10), - np.nan) +# -- NaN-aware interpolation (CuPy) ----------------------------------------- -# -- NaN-aware zoom (CuPy) --------------------------------------------------- +def _nan_aware_interp_cupy(data, out_h, out_w, order): + """CuPy variant of :func:`_nan_aware_interp_np`.""" + from cupyx.scipy.ndimage import map_coordinates as _cupy_map_coords + + iy = cupy.asarray(_block_centered_coords(data.shape[0], out_h)) + ix = cupy.asarray(_block_centered_coords(data.shape[1], out_w)) + yy, xx = cupy.meshgrid(iy, ix, indexing='ij') + coords = cupy.array([yy.ravel(), xx.ravel()]) -def _nan_aware_zoom_cupy(data, zoom_yx, order): if order == 0: - return _cupy_ndimage.zoom(data, zoom_yx, order=0, mode='nearest') + result = _cupy_map_coords(data, coords, order=0, mode='nearest') + return result.reshape(out_h, out_w) mask = cupy.isnan(data) if not mask.any(): - return _cupy_ndimage.zoom(data, zoom_yx, order=order, mode='nearest') + result = _cupy_map_coords(data, coords, order=order, mode='nearest') + return result.reshape(out_h, out_w) filled = cupy.where(mask, 0.0, data) weights = (~mask).astype(data.dtype) - z_data = _cupy_ndimage.zoom(filled, zoom_yx, order=order, mode='nearest') - z_wt = _cupy_ndimage.zoom(weights, zoom_yx, order=order, mode='nearest') + z_data = _cupy_map_coords(filled, coords, order=order, mode='nearest') + z_wt = _cupy_map_coords(weights, coords, order=order, mode='nearest') - return cupy.where(z_wt > 0.01, - z_data / cupy.maximum(z_wt, 1e-10), - cupy.nan) + result = cupy.where(z_wt > 0.01, + z_data / cupy.maximum(z_wt, 1e-10), + cupy.nan) + return result.reshape(out_h, out_w) # -- Block-aggregation kernels (NumPy, numba) -------------------------------- @@ -283,13 +316,9 @@ def _interp_block_np(block, global_in_h, global_in_w, oy = np.arange(cum_out_y[yi], cum_out_y[yi + 1], dtype=np.float64) ox = np.arange(cum_out_x[xi], cum_out_x[xi + 1], dtype=np.float64) - # Map to global input coordinates (same formula scipy.ndimage.zoom uses) - iy = (oy * (global_in_h - 1) / (global_out_h - 1) - if global_out_h > 1 - else np.full(len(oy), (global_in_h - 1) / 2.0)) - ix = (ox * (global_in_w - 1) / (global_out_w - 1) - if global_out_w > 1 - else np.full(len(ox), (global_in_w - 1) / 2.0)) + # Map to global input coordinates using block-centered formula + iy = (oy + 0.5) * (global_in_h / global_out_h) - 0.5 + ix = (ox + 0.5) * (global_in_w / global_out_w) - 0.5 # Convert to local block coordinates (overlap shifts the origin) iy_local = iy - (cum_in_y[yi] - depth) @@ -331,12 +360,9 @@ def _interp_block_cupy(block, global_in_h, global_in_w, ox = cupy.arange(int(cum_out_x[xi]), int(cum_out_x[xi + 1]), dtype=cupy.float64) - iy = (oy * (global_in_h - 1) / (global_out_h - 1) - if global_out_h > 1 - else cupy.full(len(oy), (global_in_h - 1) / 2.0)) - ix = (ox * (global_in_w - 1) / (global_out_w - 1) - if global_out_w > 1 - else cupy.full(len(ox), (global_in_w - 1) / 2.0)) + # Map to global input coordinates using block-centered formula + iy = (oy + 0.5) * (global_in_h / global_out_h) - 0.5 + ix = (ox + 0.5) * (global_in_w / global_out_w) - 0.5 iy_local = iy - float(cum_in_y[yi] - depth) ix_local = ix - float(cum_in_x[xi] - depth) @@ -410,10 +436,8 @@ def _run_numpy(data, scale_y, scale_x, method): out_h, out_w = _output_shape(*data.shape, scale_y, scale_x) if method in INTERP_METHODS: - zy = out_h / data.shape[0] - zx = out_w / data.shape[1] - return _nan_aware_zoom_np(data, (zy, zx), - INTERP_METHODS[method]).astype(np.float32) + return _nan_aware_interp_np(data, out_h, out_w, + INTERP_METHODS[method]).astype(np.float32) return _AGG_FUNCS[method](data, out_h, out_w).astype(np.float32) @@ -423,10 +447,8 @@ def _run_cupy(data, scale_y, scale_x, method): out_h, out_w = _output_shape(*data.shape, scale_y, scale_x) if method in INTERP_METHODS: - zy = out_h / data.shape[0] - zx = out_w / data.shape[1] - return _nan_aware_zoom_cupy(data, (zy, zx), - INTERP_METHODS[method]).astype(cupy.float32) + return _nan_aware_interp_cupy(data, out_h, out_w, + INTERP_METHODS[method]).astype(cupy.float32) # Aggregate: GPU reshape+reduce for integer factors, CPU fallback otherwise fy, fx = data.shape[0] / out_h, data.shape[1] / out_w diff --git a/xrspatial/tests/test_resample.py b/xrspatial/tests/test_resample.py index 226cd0ff..a02600a7 100644 --- a/xrspatial/tests/test_resample.py +++ b/xrspatial/tests/test_resample.py @@ -182,6 +182,48 @@ def test_bilinear_upsample_smooth(self, grid_8x8): # Verify interior is within tolerance of the linear gradient assert np.all(np.isfinite(out.values)) + @pytest.mark.parametrize('method', ['nearest', 'bilinear', 'cubic']) + def test_interp_coordinate_alignment_downsample(self, method): + """Interpolated values should match output coordinate labels (#1202). + + On a linear gradient where value == x-coordinate, a correct + block-centered resample produces output values equal to the output + coordinate labels (within floating-point tolerance). + """ + data = np.tile(np.arange(8, dtype=np.float32), (8, 1)) + agg = create_test_raster(data, attrs={'res': (1.0, 1.0)}, + dims=['y', 'x']) + out = resample(agg, scale_factor=0.5, method=method) + + # Output x-coords are block-centered: 0.5, 2.5, 4.5, 6.5 + # Values should match because the input is a linear gradient + np.testing.assert_allclose( + out.values[0], out.x.values, atol=0.6, + err_msg=f"{method}: values should be close to x-coordinates" + ) + # For bilinear on a linear gradient, the match should be exact + if method == 'bilinear': + np.testing.assert_allclose( + out.values[0], out.x.values, atol=1e-5, + err_msg="bilinear on linear gradient must be exact" + ) + + def test_bilinear_coordinate_alignment_upsample(self): + """Upsampled interior pixels should match coordinates on a gradient.""" + data = np.tile(np.arange(8, dtype=np.float32), (8, 1)) + agg = create_test_raster(data, attrs={'res': (1.0, 1.0)}, + dims=['y', 'x']) + out = resample(agg, scale_factor=2.0, method='bilinear') + + # Interior pixels (away from boundary clamping) should be exact + # for bilinear on a linear gradient. Skip first and last pixel + # which may be clamped by mode='nearest' boundary handling. + interior = slice(1, -1) + np.testing.assert_allclose( + out.values[0, interior], out.x.values[interior], atol=1e-4, + err_msg="bilinear: interior values should match x-coordinates" + ) + # --------------------------------------------------------------------------- # NaN handling @@ -281,6 +323,19 @@ def test_interp_parity(self, numpy_and_dask_rasters, method, sf): np.testing.assert_allclose(dk_out.values, np_out.values, atol=1e-5, equal_nan=True) + @pytest.mark.parametrize('method', ['bilinear']) + def test_dask_coordinate_alignment(self, method): + """Dask bilinear on a linear gradient should match coordinates (#1202).""" + data = np.tile(np.arange(20, dtype=np.float32), (20, 1)) + dk_agg = create_test_raster(data, backend='dask+numpy', + attrs={'res': (1.0, 1.0)}, + chunks=(8, 8)) + out = resample(dk_agg, scale_factor=0.5, method=method) + np.testing.assert_allclose( + out.values[0], out.x.values, atol=1e-4, + err_msg="dask bilinear values should match x-coordinates" + ) + @pytest.mark.parametrize('method', ['average', 'min', 'max', 'median', 'mode']) def test_aggregate_parity(self, numpy_and_dask_rasters, method): np_agg, dk_agg = numpy_and_dask_rasters