Three accuracy bugs in reproject resampling kernels

[brendancol]

Describe the bug

Three accuracy bugs in the reproject resampling kernels (_interpolate.py):

1. Cubic resampling uses border replication instead of bilinear fallback for OOB neighbors

In _resample_cubic_jit and _resample_cubic_cuda, when a cubic stencil neighbor falls outside the source array, the index gets clamped to the edge pixel (if ric < 0: ric = 0). The clamped pixel usually has a valid value, so the NaN check passes and cubic interpolation proceeds with border-replicated data instead of falling back to bilinear.

GDAL's GWKCubicResample treats out-of-bounds neighbors as invalid, triggering bilinear fallback. The difference shows up as edge artifacts in a 1-2 pixel strip around source boundaries.

_interpolate.py:112-126 (numba kernel):

for di in range(4):
    ri = r0 - 1 + di
    ric = ri
    if ric < 0:
        ric = 0          # border replication
    elif ric >= sh:
        ric = sh - 1     # border replication
    ...
    sv = src[ric, cjc]
    if sv != sv:          # NaN check passes because edge pixel is valid
        has_nan = True
        break

Fix: check whether ri (the real index) is in bounds, not just whether ric (the clamped index) holds a NaN.

2. Nearest-neighbor rounds incorrectly for negative fractional pixel coordinates

int(r + 0.5) truncates toward zero in numba/C, not toward nearest integer. For coordinates in [-1.0, -0.5), this maps to pixel 0 instead of pixel -1 (out of bounds, should be nodata).

_interpolate.py:49:

ri = int(r + 0.5)    # r = -0.7: int(-0.2) = 0, should be -1

This extends the source extent by up to half a pixel at top/left edges. Fix: use int(math.floor(r + 0.5)).

3. Legacy CuPy map_coordinates path uses hard-coded NaN contamination threshold

In _resample_cupy(), NaN values are zeroed out before map_coordinates, then contamination is detected with nan_weight > 0.1. The threshold is arbitrary. Low contamination from one distant cubic neighbor slips through, producing values biased toward zero.

_interpolate.py:843-848:

nan_weight = map_coordinates(
    nan_mask.astype(cp.float64), coords,
    order=order, mode='constant', cval=1.0
).reshape(...)
oob = oob | (nan_weight > 0.1)   # arbitrary threshold

The native CUDA kernels handle NaN inline correctly. This only affects the fallback GPU path, but the two paths disagree near nodata boundaries.

Expected behavior

• Cubic resampling near source edges falls back to bilinear (matching GDAL)
• Nearest-neighbor does not duplicate edge pixels beyond the half-pixel boundary
• Both GPU resampling paths produce consistent results near nodata

Reproducer

Create a 4x4 raster with values 1-16, resample with cubic at coordinates just outside the source (r=-0.6). Current code returns the edge pixel value; correct behavior is nodata.