Resample polish: fix cubic depth comment, tighten NaN threshold, add edge tests

xrspatial/resample.py

@@ -39,8 +39,9 @@
 # Overlap depth (input pixels) each interpolation kernel needs from
 # neighbouring chunks when processing dask arrays.  Cubic requires extra
 # depth because the B-spline prefilter is a global IIR filter whose
-# boundary transient decays as ~0.268^n; depth 10 brings it to machine
-# epsilon.
+# boundary transient decays as ~0.268^n; depth 10 brings the boundary
+# transient to ~1e-7 (sufficient for float32 output and well under
+# typical interpolation error).
 _INTERP_DEPTH = {'nearest': 1, 'bilinear': 1, 'cubic': 10}
 
 # Approximate working-set size per output cell for the eager backends:
@@ -191,7 +192,11 @@ def _nan_aware_interp_np(data, out_h, out_w, order):
     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,
+    # Gate on majority weight: an output pixel is valid only when more
+    # than half of the resampling kernel weight came from valid input
+    # pixels.  This rejects pixels lit only by cubic-kernel sidelobes
+    # leaking small positive weight from a single neighbour.
+    result = np.where(z_wt > 0.5,
                       z_data / np.maximum(z_wt, 1e-10),
                       np.nan)
     return result.reshape(out_h, out_w)
@@ -223,7 +228,8 @@ def _nan_aware_interp_cupy(data, out_h, out_w, order):
     z_data = _cupy_map_coords(filled, coords, order=order, mode='nearest')
     z_wt = _cupy_map_coords(weights, coords, order=order, mode='nearest')
 
-    result = cupy.where(z_wt > 0.01,
+    # Majority-weight gate (see _nan_aware_interp_np for rationale).
+    result = cupy.where(z_wt > 0.5,
                         z_data / cupy.maximum(z_wt, 1e-10),
                         cupy.nan)
     return result.reshape(out_h, out_w)
@@ -625,7 +631,8 @@ def _interp_block_np(block, global_in_h, global_in_w,
         weights = (~mask).astype(np.float64)
         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,
+        # Majority-weight gate (see _nan_aware_interp_np for rationale).
+        result = np.where(z_wt > 0.5,
                           z_data / np.maximum(z_wt, 1e-10), np.nan)
 
     return result.reshape(target_h, target_w).astype(np.float32)
@@ -667,7 +674,8 @@ def _interp_block_cupy(block, global_in_h, global_in_w,
         weights = (~mask).astype(cupy.float64)
         z_data = _cupy_map_coords(filled, coords, order=order, mode='nearest')
         z_wt = _cupy_map_coords(weights, coords, order=order, mode='nearest')
-        result = cupy.where(z_wt > 0.01,
+        # Majority-weight gate (see _nan_aware_interp_np for rationale).
+        result = cupy.where(z_wt > 0.5,
                             z_data / cupy.maximum(z_wt, 1e-10), cupy.nan)
 
     return result.reshape(target_h, target_w).astype(cupy.float32)

xrspatial/tests/test_resample.py

@@ -682,6 +682,12 @@ def test_target_resolution_list_form(self, grid_8x8):
             out_tuple.values, out_list.values, atol=1e-6
         )
 
+    def test_target_resolution_matches_scale_factor(self, grid_8x8):
+        # target_resolution=(2.0, 2.0) should match scale_factor=0.5.
+        a = resample(grid_8x8, target_resolution=(2.0, 2.0), method='nearest')
+        b = resample(grid_8x8, scale_factor=0.5, method='nearest')
+        np.testing.assert_allclose(a.values, b.values)
+
 
 # ---------------------------------------------------------------------------
 # Memory guard (#1295)
@@ -735,6 +741,136 @@ def test_dask_path_skips_guard(self, grid_4x4):
         assert out.shape == (400, 400)
 
 
+# ---------------------------------------------------------------------------
+# Coordinate orientation (T-3)
+# ---------------------------------------------------------------------------
+
+class TestCoordinateOrientation:
+    """Resampling must preserve the y-axis orientation of the input."""
+
+    def test_descending_y_preserved(self):
+        # 8x8 grid with strong vertical asymmetry: row 0 = 100, row 7 = 0.
+        data = np.zeros((8, 8), dtype=np.float32)
+        for r in range(8):
+            data[r, :] = (7 - r) * (100.0 / 7.0)
+        # Descending y: 10, 9, ..., 3.  x ascending.
+        y = np.arange(10, 2, -1, dtype=np.float64)
+        x = np.arange(8, dtype=np.float64)
+        agg = xr.DataArray(
+            data, dims=('y', 'x'), coords={'y': y, 'x': x},
+            attrs={'res': (1.0, 1.0)},
+        )
+
+        out = resample(agg, scale_factor=0.5, method='average')
+
+        # Output should have 4 y-coords, still descending.
+        assert out.shape == (4, 4)
+        assert np.all(np.diff(out.y.values) < 0), (
+            "output y must remain descending"
+        )
+
+        # Cell-centered y-coords for a 4-pixel output spanning the same
+        # extent as input y in [2.5, 10.5]: centers at 9.5, 7.5, 5.5, 3.5.
+        np.testing.assert_allclose(
+            out.y.values, [9.5, 7.5, 5.5, 3.5], atol=1e-6
+        )
+
+        # Top row of output should average the (high-value) top of the
+        # input, bottom row should average the (low-value) bottom.  No
+        # vertical flip.
+        assert out.values[0, 0] > out.values[-1, 0], (
+            "top of output should be greater than bottom (no flip)"
+        )
+        # Block-average of rows 0 and 1 = mean of 100 and ~85.7 = ~92.86.
+        np.testing.assert_allclose(out.values[0, 0], 92.857143, atol=1e-4)
+        # Bottom block averages rows 6 and 7 = mean of ~14.29 and 0 = ~7.14.
+        np.testing.assert_allclose(out.values[-1, 0], 7.142857, atol=1e-4)
+
+
+# ---------------------------------------------------------------------------
+# Inf input behavior (T-6)
+# ---------------------------------------------------------------------------
+
+class TestInfHandling:
+    """Lock in the current behavior: infs propagate, they are not NaN-coerced.
+
+    The NaN-aware interpolation path treats infs the same as any other
+    finite value (they pass through ``np.isnan`` checks), so the
+    interpolation kernel includes them in its weighted sum and the
+    aggregate kernel includes them in min / max / mean.  Output values
+    that depend on an inf input cell come back as +/- inf.
+    """
+
+    @pytest.fixture
+    def grid_with_infs(self):
+        data = np.ones((6, 6), dtype=np.float32)
+        data[1, 1] = np.inf
+        data[4, 4] = -np.inf
+        return create_test_raster(data, attrs={'res': (1.0, 1.0)})
+
+    def test_bilinear_propagates_inf(self, grid_with_infs):
+        out = resample(grid_with_infs, scale_factor=0.5, method='bilinear')
+        # +inf and -inf both reach the output, no NaN coercion.
+        assert np.isposinf(out.values).any()
+        assert np.isneginf(out.values).any()
+        assert not np.isnan(out.values).any()
+
+    def test_average_propagates_inf(self, grid_with_infs):
+        out = resample(grid_with_infs, scale_factor=0.5, method='average')
+        assert np.isposinf(out.values).any()
+        assert np.isneginf(out.values).any()
+        assert not np.isnan(out.values).any()
+
+    def test_nearest_passes_inf_through(self, grid_with_infs):
+        out = resample(grid_with_infs, scale_factor=0.5, method='nearest')
+        # Nearest only picks one input cell per output pixel, so at
+        # least one infinite input survives in the output.
+        assert np.isinf(out.values).any()
+
+
+# ---------------------------------------------------------------------------
+# 4D input rejection message (T-7)
+# ---------------------------------------------------------------------------
+
+class TestDimensionRejection:
+    def test_4d_input_error_message(self):
+        # 3D (band, y, x) is now valid; 4D should still be rejected.
+        data = np.zeros((2, 2, 4, 4), dtype=np.float32)
+        agg = xr.DataArray(
+            data, dims=('time', 'band', 'y', 'x'),
+            coords={
+                'time': [0, 1],
+                'band': [1, 2],
+                'y': np.arange(4, dtype=np.float64),
+                'x': np.arange(4, dtype=np.float64),
+            },
+            attrs={'res': (1.0, 1.0)},
+        )
+        with pytest.raises(ValueError, match="must be 2D"):
+            resample(agg, scale_factor=0.5, method='nearest')
+
+
+# ---------------------------------------------------------------------------
+# Bit reproducibility (T-9)
+# ---------------------------------------------------------------------------
+
+class TestReproducibility:
+    """Repeat calls with identical inputs must produce identical outputs."""
+
+    @pytest.fixture
+    def random_grid(self):
+        data = np.random.RandomState(1471).rand(20, 20).astype(np.float32)
+        return create_test_raster(data, attrs={'res': (1.0, 1.0)})
+
+    @pytest.mark.parametrize('method', ['bilinear', 'cubic', 'average'])
+    def test_bit_identical_repeat(self, random_grid, method):
+        out1 = resample(random_grid, scale_factor=0.5, method=method)
+        out2 = resample(random_grid, scale_factor=0.5, method=method)
+        assert np.array_equal(out1.values, out2.values), (
+            f"{method} resample is not bit-identical across runs"
+        )
+
+
 # ---------------------------------------------------------------------------
 # Inlined dask aggregate kernel (#1463)
 # ---------------------------------------------------------------------------
@@ -943,20 +1079,3 @@ def test_explicit_nodata_overrides_attr(self):
         # Without override the attr would say -999 (no match) and -1 would
         # leak through; with override the top-left output pixel is NaN.
         assert np.isnan(out.values[0, 0])
-
-
-# ---------------------------------------------------------------------------
-# target_resolution as 2-tuple (issue #1466)
-# ---------------------------------------------------------------------------
-
-class TestTargetResolutionTuple:
-    def test_tuple_resolution_independent_axes(self, grid_8x8):
-        # 8x8 grid with res=(1, 1) -> target (2, 4) -> output (4, 2).
-        out = resample(grid_8x8, target_resolution=(2.0, 4.0))
-        assert out.shape == (4, 2)
-
-    def test_tuple_resolution_matches_scale_factor(self, grid_8x8):
-        # target_resolution=(2.0, 2.0) should match scale_factor=0.5.
-        a = resample(grid_8x8, target_resolution=(2.0, 2.0), method='nearest')
-        b = resample(grid_8x8, scale_factor=0.5, method='nearest')
-        np.testing.assert_allclose(a.values, b.values)