Fix three accuracy bugs in reproject resampling kernels (#1087)

* Add design spec for geotiff performance and memory controls

Covers dtype on read, compression_level on write, and VRT tiled
output from to_geotiff for streaming dask writes.

* Fix three accuracy bugs in reproject resampling kernels (#1086)

1. Nearest-neighbor: use floor(r + 0.5) instead of int(r + 0.5) to
   round correctly for negative fractional pixel coordinates.
   int() truncates toward zero, which maps [-1.0, -0.5) to pixel 0
   instead of out-of-bounds.

2. Cubic: check real neighbor indices against source bounds instead of
   clamping OOB indices to the edge pixel. Border replication masked
   the OOB condition, preventing the bilinear fallback that GDAL uses
   when the 4x4 stencil extends outside the source.

3. CuPy map_coordinates fallback: tighten the NaN contamination
   threshold from 0.1 to 1e-6 so that small NaN weights from distant
   cubic neighbors are caught, matching the native CUDA kernel behavior.

* Address code review: edge-case tests and threshold comment (#1086)

- Test nearest-neighbor at exact boundary r=-0.5
- Test cubic fallback at far (bottom-right) edge
- Comment the 1e-6 NaN threshold rationale

Author: brendancol

xrspatial/reproject/_interpolate.py

@@ -46,16 +46,11 @@ def _resample_nearest_jit(src, row_coords, col_coords, nodata):
             if r < -1.0 or r > sh or c < -1.0 or c > sw:
                 out[i, j] = nodata
                 continue
-            ri = int(r + 0.5)
-            ci = int(c + 0.5)
-            if ri < 0:
-                ri = 0
-            if ri >= sh:
-                ri = sh - 1
-            if ci < 0:
-                ci = 0
-            if ci >= sw:
-                ci = sw - 1
+            ri = int(np.floor(r + 0.5))
+            ci = int(np.floor(c + 0.5))
+            if ri < 0 or ri >= sh or ci < 0 or ci >= sw:
+                out[i, j] = nodata
+                continue
             v = src[ri, ci]
             # NaN check: works for float64
             if v != v:
@@ -109,21 +104,17 @@ def _resample_cubic_jit(src, row_coords, col_coords, nodata):
             has_nan = False
             for di in range(4):
                 ri = r0 - 1 + di
-                ric = ri
-                if ric < 0:
-                    ric = 0
-                elif ric >= sh:
-                    ric = sh - 1
+                if ri < 0 or ri >= sh:
+                    has_nan = True
+                    break
                 # Interpolate along columns for this row
                 rv = 0.0
                 for dj in range(4):
                     cj = c0 - 1 + dj
-                    cjc = cj
-                    if cjc < 0:
-                        cjc = 0
-                    elif cjc >= sw:
-                        cjc = sw - 1
-                    sv = src[ric, cjc]
+                    if cj < 0 or cj >= sw:
+                        has_nan = True
+                        break
+                    sv = src[ri, cj]
                     if sv != sv:  # NaN check
                         has_nan = True
                         break
@@ -339,16 +330,11 @@ def _resample_nearest_cuda(src, row_coords, col_coords, out, nodata):
         if r < -1.0 or r > sh or c < -1.0 or c > sw:
             out[i, j] = nodata
             return
-        ri = int(r + 0.5)
-        ci = int(c + 0.5)
-        if ri < 0:
-            ri = 0
-        if ri >= sh:
-            ri = sh - 1
-        if ci < 0:
-            ci = 0
-        if ci >= sw:
-            ci = sw - 1
+        ri = int(math.floor(r + 0.5))
+        ci = int(math.floor(c + 0.5))
+        if ri < 0 or ri >= sh or ci < 0 or ci >= sw:
+            out[i, j] = nodata
+            return
         v = src[ri, ci]
         # NaN check
         if v != v:
@@ -453,6 +439,11 @@ def _resample_cubic_cuda(src, row_coords, col_coords, out, nodata):
         val = 0.0
         has_nan = False
 
+        # If any of the 4x4 stencil neighbors is outside source, fall
+        # back to bilinear (matches GDAL behavior for boundary pixels).
+        if r0 - 1 < 0 or r0 + 2 >= sh or c0 - 1 < 0 or c0 + 2 >= sw:
+            has_nan = True
+
         # Row 0
         ric = r0 - 1
         if ric < 0:
@@ -844,7 +835,10 @@ def _resample_cupy(source_window, src_row_coords, src_col_coords,
             nan_mask.astype(cp.float64), coords,
             order=order, mode='constant', cval=1.0
         ).reshape(src_row_coords.shape)
-        oob = oob | (nan_weight > 0.1)
+        # Any nonzero weight means NaN pixels contributed to the output.
+        # Use 1e-6 instead of 0.0 to absorb floating-point interpolation
+        # noise from map_coordinates on the binary mask.
+        oob = oob | (nan_weight > 1e-6)
 
     result[oob] = nodata
 

xrspatial/tests/test_reproject.py

@@ -195,6 +195,78 @@ def test_resample_oob_fills_nodata(self):
         result = _resample_numpy(src, rows, cols, nodata=-999)
         assert result[0, 0] == -999
 
+    def test_nearest_negative_rounding(self):
+        """int(r + 0.5) must round toward -inf, not toward zero (#1086)."""
+        from xrspatial.reproject._interpolate import _resample_numpy
+        src = np.arange(1, 17, dtype=np.float64).reshape(4, 4)
+        # r = -0.6 is beyond the half-pixel boundary of pixel 0 -> nodata
+        rows = np.array([[-0.6]])
+        cols = np.array([[1.0]])
+        result = _resample_numpy(src, rows, cols, resampling='nearest', nodata=-999)
+        assert result[0, 0] == -999, (
+            f"r=-0.6 should be nodata, got {result[0, 0]}"
+        )
+        # r = -0.4 is within pixel 0's domain -> pixel 0
+        rows2 = np.array([[-0.4]])
+        result2 = _resample_numpy(src, rows2, cols, resampling='nearest', nodata=-999)
+        assert result2[0, 0] == src[0, 1], (
+            f"r=-0.4 should map to pixel 0, got {result2[0, 0]}"
+        )
+        # r = -0.5 is exactly on the half-pixel boundary: floor(-0.5+0.5)=0 -> pixel 0
+        rows3 = np.array([[-0.5]])
+        result3 = _resample_numpy(src, rows3, cols, resampling='nearest', nodata=-999)
+        assert result3[0, 0] == src[0, 1], (
+            f"r=-0.5 should map to pixel 0, got {result3[0, 0]}"
+        )
+
+    def test_cubic_oob_fallback(self):
+        """Cubic must fall back to bilinear when stencil extends outside source (#1086)."""
+        from xrspatial.reproject._interpolate import _resample_numpy
+        # 6x6 source with a gradient
+        src = np.arange(36, dtype=np.float64).reshape(6, 6)
+        # Query at r=0.5, c=0.5: cubic stencil needs row -1, which is OOB.
+        # Should fall back to bilinear using pixels (0,0),(0,1),(1,0),(1,1).
+        rows = np.array([[0.5]])
+        cols = np.array([[0.5]])
+        cubic_result = _resample_numpy(src, rows, cols, resampling='cubic', nodata=-999)
+        bilinear_result = _resample_numpy(src, rows, cols, resampling='bilinear', nodata=-999)
+        # At the boundary, cubic should produce the same result as bilinear
+        np.testing.assert_allclose(
+            cubic_result, bilinear_result, atol=1e-10,
+            err_msg="Cubic near boundary should fall back to bilinear"
+        )
+        # Interior query at r=2.5, c=2.5: full stencil fits, cubic should differ from bilinear
+        rows_int = np.array([[2.5]])
+        cols_int = np.array([[2.5]])
+        cubic_int = _resample_numpy(src, rows_int, cols_int, resampling='cubic', nodata=-999)
+        bilinear_int = _resample_numpy(src, rows_int, cols_int, resampling='bilinear', nodata=-999)
+        # For a linear gradient, cubic and bilinear should agree closely
+        # but the point is the code path exercises the non-fallback branch
+        assert cubic_int[0, 0] != -999
+
+    def test_cubic_oob_fallback_far_edge(self):
+        """Cubic at bottom-right boundary: stencil needs row sh, same fallback (#1086)."""
+        from xrspatial.reproject._interpolate import _resample_numpy
+        src = np.arange(36, dtype=np.float64).reshape(6, 6)
+        # r=4.5: cubic stencil needs row 6 (= sh), which is OOB
+        rows = np.array([[4.5]])
+        cols = np.array([[4.5]])
+        cubic = _resample_numpy(src, rows, cols, resampling='cubic', nodata=-999)
+        bilinear = _resample_numpy(src, rows, cols, resampling='bilinear', nodata=-999)
+        np.testing.assert_allclose(cubic, bilinear, atol=1e-10)
+
+    def test_cubic_oob_bilinear_fallback_renormalizes(self):
+        """Cubic at (-0.8,-0.8): stencil OOB triggers bilinear, which
+        finds pixel (0,0) as the only valid neighbor and returns it (#1086)."""
+        from xrspatial.reproject._interpolate import _resample_numpy
+        src = np.arange(1, 17, dtype=np.float64).reshape(4, 4)
+        rows = np.array([[-0.8]])
+        cols = np.array([[-0.8]])
+        result = _resample_numpy(src, rows, cols, resampling='cubic', nodata=-999)
+        # bilinear fallback: r0=-1 (OOB), r1=0, c0=-1 (OOB), c1=0
+        # only (r1,c1)=(0,0) is valid -> returns src[0,0]=1.0
+        assert result[0, 0] == 1.0
+
     def test_invalid_resampling(self):
         from xrspatial.reproject._interpolate import _validate_resampling
         with pytest.raises(ValueError, match="resampling"):