Fix cubic NaN, add merge tests, validate grids, improve docs (#1045)

Cubic NaN handling:
- When any of the 16 Catmull-Rom neighbors is NaN, falls back to
  bilinear with weight renormalization instead of returning nodata.
  Eliminates the one-pixel nodata halo around NaN regions that
  cubic resampling previously produced.

Merge strategy tests:
- Added end-to-end tests for last, max, min strategies (were only
  tested at the internal _merge_arrays_numpy level).

Datum grid validation:
- load_grid() now validates band shapes match, grid is >= 2x2,
  and bounds are sensible. Invalid grids return None (Helmert
  fallback) instead of producing garbage.

Documentation:
- reproject() and merge() docstrings now note output CRS is WKT
  format in attrs['crs'], and merge() documents CRS selection
  when target_crs=None.

Author: brendancol

xrspatial/reproject/__init__.py

@@ -429,6 +429,7 @@ def reproject(
     Returns
     -------
     xr.DataArray
+        The output's ``attrs['crs']`` is always stored in WKT format.
     """
     from ._crs_utils import _require_pyproj
 
@@ -805,7 +806,8 @@ def merge(
     rasters : list of xr.DataArray
         Input rasters to merge.
     target_crs : optional
-        Target CRS. Defaults to the CRS of the first raster.
+        Target CRS. When None, the CRS of the first raster in the list
+        is used.
     resolution : float or (float, float) or None
         Output resolution in target CRS units.
     bounds : (left, bottom, right, top) or None
@@ -824,6 +826,7 @@ def merge(
     Returns
     -------
     xr.DataArray
+        The output's ``attrs['crs']`` is always stored in WKT format.
     """
     from ._crs_utils import _require_pyproj
 

xrspatial/reproject/_datum_grids.py

@@ -169,6 +169,13 @@ def load_grid(grid_key):
             b = ds.bounds
             bounds = (b.left, b.bottom, b.right, b.top)
             h, w = ds.height, ds.width
+            # Validate grid shape and bounds
+            if dlat.shape != dlon.shape:
+                return None
+            if h < 2 or w < 2:
+                return None
+            if b.left >= b.right or b.bottom >= b.top:
+                return None
             # Compute resolution from bounds and shape (avoids ds.res ordering ambiguity)
             res_x = (b.right - b.left) / w if w > 1 else 0.25
             res_y = (b.top - b.bottom) / h if h > 1 else 0.25
@@ -191,10 +198,19 @@ def load_grid(grid_key):
     else:
         return None
 
+    # Validate grid shape and bounds
+    if dlat.shape != dlon.shape:
+        return None
+    if dlat.shape[0] < 2 or dlat.shape[1] < 2:
+        return None
+
     y_coords = da.coords['y'].values
     x_coords = da.coords['x'].values
     bounds = (float(x_coords[0]), float(y_coords[-1]),
               float(x_coords[-1]), float(y_coords[0]))
+    left, bottom, right, top = bounds
+    if left >= right or bottom >= top:
+        return None
     res_x = abs(float(x_coords[1] - x_coords[0])) if len(x_coords) > 1 else 0.25
     res_y = abs(float(y_coords[1] - y_coords[0])) if len(y_coords) > 1 else 0.25
     return dlat, dlon, bounds, (res_x, res_y)

xrspatial/reproject/_interpolate.py

@@ -67,10 +67,11 @@ def _resample_nearest_jit(src, row_coords, col_coords, nodata):
 
 @njit(nogil=True, cache=True)
 def _resample_cubic_jit(src, row_coords, col_coords, nodata):
-    """Catmull-Rom cubic resampling with NaN propagation.
+    """Catmull-Rom cubic resampling with NaN-aware fallback to bilinear.
 
     Separable: interpolate 4 row-slices along columns, then combine
-    along rows.  Handles NaN inline (no second pass needed).
+    along rows.  When any of the 16 neighbors is NaN, falls back to
+    bilinear with weight renormalization (matching GDAL behavior).
     """
     h_out, w_out = row_coords.shape
     sh, sw = src.shape
@@ -145,7 +146,51 @@ def _resample_cubic_jit(src, row_coords, col_coords, nodata):
                 else:
                     val += rv * wr3
 
-            out[i, j] = nodata if has_nan else val
+            if not has_nan:
+                out[i, j] = val
+            else:
+                # Fall back to bilinear with weight renormalization
+                r1 = r0 + 1
+                c1 = c0 + 1
+                dr = r - r0
+                dc = c - c0
+
+                w00 = (1.0 - dr) * (1.0 - dc)
+                w01 = (1.0 - dr) * dc
+                w10 = dr * (1.0 - dc)
+                w11 = dr * dc
+
+                accum = 0.0
+                wsum = 0.0
+
+                if 0 <= r0 < sh and 0 <= c0 < sw:
+                    v = src[r0, c0]
+                    if v == v:
+                        accum += w00 * v
+                        wsum += w00
+
+                if 0 <= r0 < sh and 0 <= c1 < sw:
+                    v = src[r0, c1]
+                    if v == v:
+                        accum += w01 * v
+                        wsum += w01
+
+                if 0 <= r1 < sh and 0 <= c0 < sw:
+                    v = src[r1, c0]
+                    if v == v:
+                        accum += w10 * v
+                        wsum += w10
+
+                if 0 <= r1 < sh and 0 <= c1 < sw:
+                    v = src[r1, c1]
+                    if v == v:
+                        accum += w11 * v
+                        wsum += w11
+
+                if wsum > 1e-10:
+                    out[i, j] = accum / wsum
+                else:
+                    out[i, j] = nodata
     return out
 
 

xrspatial/tests/test_reproject.py

@@ -577,6 +577,54 @@ def test_merge_invalid_strategy(self):
         with pytest.raises(ValueError, match="strategy"):
             merge([raster], strategy='median')
 
+    def test_merge_strategy_last(self):
+        """merge() with strategy='last' uses the last valid value."""
+        from xrspatial.reproject import merge
+        a = _make_raster(
+            np.full((16, 16), 10.0), x_range=(-5, 5), y_range=(-5, 5)
+        )
+        b = _make_raster(
+            np.full((16, 16), 20.0), x_range=(-5, 5), y_range=(-5, 5)
+        )
+        result = merge([a, b], strategy='last', resolution=1.0)
+        vals = result.values
+        interior = vals[2:-2, 2:-2]
+        valid = ~np.isnan(interior) & (interior != 0)
+        if valid.any():
+            np.testing.assert_allclose(interior[valid], 20.0, atol=1.0)
+
+    def test_merge_strategy_max(self):
+        """merge() with strategy='max' takes the maximum."""
+        from xrspatial.reproject import merge
+        a = _make_raster(
+            np.full((16, 16), 10.0), x_range=(-5, 5), y_range=(-5, 5)
+        )
+        b = _make_raster(
+            np.full((16, 16), 20.0), x_range=(-5, 5), y_range=(-5, 5)
+        )
+        result = merge([a, b], strategy='max', resolution=1.0)
+        vals = result.values
+        interior = vals[2:-2, 2:-2]
+        valid = ~np.isnan(interior) & (interior != 0)
+        if valid.any():
+            np.testing.assert_allclose(interior[valid], 20.0, atol=1.0)
+
+    def test_merge_strategy_min(self):
+        """merge() with strategy='min' takes the minimum."""
+        from xrspatial.reproject import merge
+        a = _make_raster(
+            np.full((16, 16), 10.0), x_range=(-5, 5), y_range=(-5, 5)
+        )
+        b = _make_raster(
+            np.full((16, 16), 20.0), x_range=(-5, 5), y_range=(-5, 5)
+        )
+        result = merge([a, b], strategy='min', resolution=1.0)
+        vals = result.values
+        interior = vals[2:-2, 2:-2]
+        valid = ~np.isnan(interior) & (interior != 0)
+        if valid.any():
+            np.testing.assert_allclose(interior[valid], 10.0, atol=1.0)
+
     @pytest.mark.skipif(not HAS_DASK, reason="dask required")
     def test_merge_dask(self):
         from xrspatial.reproject import merge