Add overview resampling options: nearest, min, max, median, mode, cubic

_make_overview() now accepts a method parameter instead of hardcoding
2x2 block averaging. Available methods:

- mean (default): nanmean of each 2x2 block, same as before
- nearest: top-left pixel of each block (no interpolation)
- min/max: nanmin/nanmax of each block
- median: nanmedian of each block
- mode: most frequent value per block (for classified rasters)
- cubic: scipy.ndimage.zoom with order=3 (requires scipy)

All methods work on both 2D and 3D (multi-band) arrays. Exposed via
overview_resampling= parameter on write() and write_geotiff().

12 new tests covering each method, NaN handling, multi-band, COG
round-trips with nearest and mode, the public API, and error on
invalid method names.

Author: brendancol

xrspatial/geotiff/__init__.py

@@ -186,7 +186,8 @@ def write_geotiff(data: xr.DataArray | np.ndarray, path: str, *,
                   tile_size: int = 256,
                   predictor: bool = False,
                   cog: bool = False,
-                  overview_levels: list[int] | None = None) -> None:
+                  overview_levels: list[int] | None = None,
+                  overview_resampling: str = 'mean') -> None:
     """Write data as a GeoTIFF or Cloud Optimized GeoTIFF.
 
     Parameters
@@ -211,6 +212,9 @@ def write_geotiff(data: xr.DataArray | np.ndarray, path: str, *,
         Write as Cloud Optimized GeoTIFF.
     overview_levels : list[int] or None
         Overview decimation factors. Only used when cog=True.
+    overview_resampling : str
+        Resampling method for overviews: 'mean' (default), 'nearest',
+        'min', 'max', 'median', 'mode', or 'cubic'.
     """
     geo_transform = None
     epsg = crs
@@ -243,6 +247,7 @@ def write_geotiff(data: xr.DataArray | np.ndarray, path: str, *,
         predictor=predictor,
         cog=cog,
         overview_levels=overview_levels,
+        overview_resampling=overview_resampling,
         raster_type=raster_type,
     )
 

xrspatial/geotiff/_writer.py

@@ -69,40 +69,85 @@ def _compression_tag(compression_name: str) -> int:
     return _map[name]
 
 
-def _make_overview(arr: np.ndarray) -> np.ndarray:
-    """Generate a 2x decimated overview using 2x2 block averaging.
+OVERVIEW_METHODS = ('mean', 'nearest', 'min', 'max', 'median', 'mode', 'cubic')
+
+
+def _block_reduce_2d(arr2d, method):
+    """2x block-reduce a single 2D plane using *method*."""
+    h, w = arr2d.shape
+    h2 = (h // 2) * 2
+    w2 = (w // 2) * 2
+    cropped = arr2d[:h2, :w2]
+    oh, ow = h2 // 2, w2 // 2
+
+    if method == 'nearest':
+        # Top-left pixel of each 2x2 block
+        return cropped[::2, ::2].copy()
+
+    if method == 'cubic':
+        try:
+            from scipy.ndimage import zoom
+        except ImportError:
+            raise ImportError(
+                "scipy is required for cubic overview resampling. "
+                "Install it with: pip install scipy")
+        return zoom(arr2d, 0.5, order=3).astype(arr2d.dtype)
+
+    if method == 'mode':
+        # Most-common value per 2x2 block (useful for classified rasters)
+        blocks = cropped.reshape(oh, 2, ow, 2).transpose(0, 2, 1, 3).reshape(oh, ow, 4)
+        out = np.empty((oh, ow), dtype=arr2d.dtype)
+        for r in range(oh):
+            for c in range(ow):
+                vals, counts = np.unique(blocks[r, c], return_counts=True)
+                out[r, c] = vals[counts.argmax()]
+        return out
+
+    # Block reshape for mean/min/max/median
+    if arr2d.dtype.kind == 'f':
+        blocks = cropped.reshape(oh, 2, ow, 2)
+    else:
+        blocks = cropped.astype(np.float64).reshape(oh, 2, ow, 2)
+
+    if method == 'mean':
+        result = np.nanmean(blocks, axis=(1, 3))
+    elif method == 'min':
+        result = np.nanmin(blocks, axis=(1, 3))
+    elif method == 'max':
+        result = np.nanmax(blocks, axis=(1, 3))
+    elif method == 'median':
+        flat = blocks.transpose(0, 2, 1, 3).reshape(oh, ow, 4)
+        result = np.nanmedian(flat, axis=2)
+    else:
+        raise ValueError(
+            f"Unknown overview resampling method: {method!r}. "
+            f"Use one of: {OVERVIEW_METHODS}")
+
+    if arr2d.dtype.kind != 'f':
+        return np.round(result).astype(arr2d.dtype)
+    return result.astype(arr2d.dtype)
+
+
+def _make_overview(arr: np.ndarray, method: str = 'mean') -> np.ndarray:
+    """Generate a 2x decimated overview.
 
     Parameters
     ----------
     arr : np.ndarray
         2D or 3D (height, width, bands) array.
+    method : str
+        Resampling method: 'mean' (default), 'nearest', 'min', 'max',
+        'median', 'mode', or 'cubic'.
 
     Returns
     -------
     np.ndarray
         Half-resolution array.
     """
-    h, w = arr.shape[:2]
-    h2 = (h // 2) * 2
-    w2 = (w // 2) * 2
-    cropped = arr[:h2, :w2]
-
     if arr.ndim == 3:
-        # Multi-band: average each band independently
-        bands = arr.shape[2]
-        if arr.dtype.kind == 'f':
-            blocks = cropped.reshape(h2 // 2, 2, w2 // 2, 2, bands)
-            return np.nanmean(blocks, axis=(1, 3)).astype(arr.dtype)
-        else:
-            blocks = cropped.astype(np.float64).reshape(h2 // 2, 2, w2 // 2, 2, bands)
-            return np.round(blocks.mean(axis=(1, 3))).astype(arr.dtype)
-    else:
-        if arr.dtype.kind == 'f':
-            blocks = cropped.reshape(h2 // 2, 2, w2 // 2, 2)
-            return np.nanmean(blocks, axis=(1, 3)).astype(arr.dtype)
-        else:
-            blocks = cropped.astype(np.float64).reshape(h2 // 2, 2, w2 // 2, 2)
-            return np.round(blocks.mean(axis=(1, 3))).astype(arr.dtype)
+        bands = [_block_reduce_2d(arr[:, :, b], method) for b in range(arr.shape[2])]
+        return np.stack(bands, axis=2)
+    return _block_reduce_2d(arr, method)
 
 
 # ---------------------------------------------------------------------------
@@ -619,6 +664,7 @@ def write(data: np.ndarray, path: str, *,
           predictor: bool = False,
           cog: bool = False,
           overview_levels: list[int] | None = None,
+          overview_resampling: str = 'mean',
           raster_type: int = 1) -> None:
     """Write a numpy array as a GeoTIFF or COG.
 
@@ -676,7 +722,7 @@ def write(data: np.ndarray, path: str, *,
 
         current = data
         for _ in overview_levels:
-            current = _make_overview(current)
+            current = _make_overview(current, method=overview_resampling)
             oh, ow = current.shape[:2]
             if tiled:
                 o_off, o_bc, o_data = _write_tiled(current, comp_tag, predictor, tile_size)

xrspatial/geotiff/tests/test_features.py

@@ -256,6 +256,157 @@ def test_zstd_public_api(self, tmp_path):
         np.testing.assert_array_equal(result.values, arr)
 
 
+# -----------------------------------------------------------------------
+# Overview resampling methods
+# -----------------------------------------------------------------------
+
+class TestOverviewResampling:
+
+    def test_mean_default(self, tmp_path):
+        """Default mean resampling produces correct 2x2 block averages."""
+        from xrspatial.geotiff._writer import _make_overview
+        arr = np.array([[1, 3, 5, 7],
+                        [2, 4, 6, 8],
+                        [10, 20, 30, 40],
+                        [10, 20, 30, 40]], dtype=np.float32)
+        ov = _make_overview(arr, 'mean')
+        assert ov.shape == (2, 2)
+        # (1+3+2+4)/4 = 2.5
+        assert ov[0, 0] == pytest.approx(2.5)
+
+    def test_nearest(self, tmp_path):
+        """Nearest resampling picks top-left pixel of each 2x2 block."""
+        from xrspatial.geotiff._writer import _make_overview
+        arr = np.array([[10, 20, 30, 40],
+                        [50, 60, 70, 80],
+                        [90, 100, 110, 120],
+                        [130, 140, 150, 160]], dtype=np.uint8)
+        ov = _make_overview(arr, 'nearest')
+        assert ov.shape == (2, 2)
+        assert ov[0, 0] == 10
+        assert ov[0, 1] == 30
+        assert ov[1, 0] == 90
+        assert ov[1, 1] == 110
+
+    def test_min(self, tmp_path):
+        from xrspatial.geotiff._writer import _make_overview
+        arr = np.array([[10, 1, 5, 3],
+                        [20, 2, 6, 4],
+                        [30, 3, 7, 5],
+                        [40, 4, 8, 6]], dtype=np.float32)
+        ov = _make_overview(arr, 'min')
+        assert ov[0, 0] == pytest.approx(1.0)
+        assert ov[0, 1] == pytest.approx(3.0)
+
+    def test_max(self, tmp_path):
+        from xrspatial.geotiff._writer import _make_overview
+        arr = np.array([[10, 1, 5, 3],
+                        [20, 2, 6, 4],
+                        [30, 3, 7, 5],
+                        [40, 4, 8, 6]], dtype=np.float32)
+        ov = _make_overview(arr, 'max')
+        assert ov[0, 0] == pytest.approx(20.0)
+        assert ov[1, 1] == pytest.approx(8.0)
+
+    def test_median(self, tmp_path):
+        from xrspatial.geotiff._writer import _make_overview
+        arr = np.array([[1, 2, 10, 20],
+                        [3, 100, 30, 40],
+                        [0, 0, 0, 0],
+                        [0, 0, 0, 0]], dtype=np.float32)
+        ov = _make_overview(arr, 'median')
+        assert ov.shape == (2, 2)
+        # median of [1, 2, 3, 100] = 2.5
+        assert ov[0, 0] == pytest.approx(2.5)
+
+    def test_mode(self, tmp_path):
+        """Mode picks the most common value in each 2x2 block."""
+        from xrspatial.geotiff._writer import _make_overview
+        arr = np.array([[1, 1, 2, 3],
+                        [1, 2, 2, 2],
+                        [5, 5, 5, 6],
+                        [5, 7, 6, 6]], dtype=np.uint8)
+        ov = _make_overview(arr, 'mode')
+        assert ov[0, 0] == 1   # 1 appears 3 times
+        assert ov[0, 1] == 2   # 2 appears 3 times
+        assert ov[1, 0] == 5   # 5 appears 3 times
+        assert ov[1, 1] == 6   # 6 appears 3 times
+
+    def test_mean_with_nan(self, tmp_path):
+        """Mean resampling ignores NaN values."""
+        from xrspatial.geotiff._writer import _make_overview
+        arr = np.array([[np.nan, 2, 4, 6],
+                        [1, 3, np.nan, 8],
+                        [10, 20, 30, 40],
+                        [10, 20, 30, 40]], dtype=np.float32)
+        ov = _make_overview(arr, 'mean')
+        # nanmean([nan, 2, 1, 3]) = 2.0
+        assert ov[0, 0] == pytest.approx(2.0)
+
+    def test_multiband(self, tmp_path):
+        """Resampling works on 3D (multi-band) arrays."""
+        from xrspatial.geotiff._writer import _make_overview
+        arr = np.zeros((4, 4, 3), dtype=np.uint8)
+        arr[:, :, 0] = 100
+        arr[:, :, 1] = 200
+        arr[:, :, 2] = 50
+        ov = _make_overview(arr, 'mean')
+        assert ov.shape == (2, 2, 3)
+        assert ov[0, 0, 0] == 100
+        assert ov[0, 0, 1] == 200
+        assert ov[0, 0, 2] == 50
+
+    def test_cog_round_trip_nearest(self, tmp_path):
+        """COG with nearest resampling writes and reads back."""
+        arr = np.arange(256, dtype=np.float32).reshape(16, 16)
+        path = str(tmp_path / 'cog_nearest.tif')
+        write(arr, path, compression='deflate', tiled=True, tile_size=8,
+              cog=True, overview_levels=[1], overview_resampling='nearest')
+
+        result, _ = read_to_array(path)
+        np.testing.assert_array_equal(result, arr)
+
+    def test_cog_round_trip_mode(self, tmp_path):
+        """COG with mode resampling for classified data."""
+        arr = np.array([[0, 0, 1, 1, 2, 2, 3, 3],
+                        [0, 0, 1, 1, 2, 2, 3, 3],
+                        [4, 4, 5, 5, 6, 6, 7, 7],
+                        [4, 4, 5, 5, 6, 6, 7, 7],
+                        [0, 0, 1, 1, 2, 2, 3, 3],
+                        [0, 0, 1, 1, 2, 2, 3, 3],
+                        [4, 4, 5, 5, 6, 6, 7, 7],
+                        [4, 4, 5, 5, 6, 6, 7, 7]], dtype=np.uint8)
+        path = str(tmp_path / 'cog_mode.tif')
+        write(arr, path, compression='deflate', tiled=True, tile_size=4,
+              cog=True, overview_levels=[1], overview_resampling='mode')
+
+        # Full res should be exact
+        result, _ = read_to_array(path)
+        np.testing.assert_array_equal(result, arr)
+
+        # Overview should have mode-reduced values
+        ov, _ = read_to_array(path, overview_level=1)
+        assert ov.shape == (4, 4)
+        assert ov[0, 0] == 0
+        assert ov[0, 1] == 1
+
+    def test_write_geotiff_api(self, tmp_path):
+        """overview_resampling kwarg works through the public API."""
+        arr = np.arange(64, dtype=np.float32).reshape(8, 8)
+        path = str(tmp_path / 'api_nearest.tif')
+        write_geotiff(arr, path, compression='deflate',
+                      cog=True, overview_resampling='nearest')
+
+        result = read_geotiff(path)
+        np.testing.assert_array_equal(result.values, arr)
+
+    def test_invalid_method(self):
+        from xrspatial.geotiff._writer import _make_overview
+        arr = np.ones((4, 4), dtype=np.float32)
+        with pytest.raises(ValueError, match="Unknown overview resampling"):
+            _make_overview(arr, 'bicubic_spline')
+
+
 # -----------------------------------------------------------------------
 # BigTIFF write
 # -----------------------------------------------------------------------