From c047e761908315fa3b25e4cfdf8c61bc5369b974 Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 20:53:52 +0200
Subject: [PATCH 01/22] perf(vectorize): auto-tile to_feature_collection by
 array size instead of always reading whole

to_feature_collection read the whole raster into memory by default (tile=False).
Make tile default to None and auto-select: read the whole array when it is at most
256 MiB, else read it tile by tile, so a huge or /vsicurl source is never
materialised whole (#969). The tiled path already restores row-major order, so the
result is byte-identical; an explicit tile=True/False still overrides the choice.

Refs #969
---
 src/pyramids/dataset/engines/vectorize.py     | 32 ++++++++++---
 .../analysis/test_dataset_vectorize.py        | 48 +++++++++++++++++++
 2 files changed, 73 insertions(+), 7 deletions(-)

diff --git a/src/pyramids/dataset/engines/vectorize.py b/src/pyramids/dataset/engines/vectorize.py
index 5a82f301cb..b3a4491ba0 100644
--- a/src/pyramids/dataset/engines/vectorize.py
+++ b/src/pyramids/dataset/engines/vectorize.py
@@ -52,6 +52,12 @@
 # out of order can be restored to row-major before the frame is returned.
 _CELL_ORDER = "_pyramids_cell_order"
 
+# When `to_feature_collection(tile=None)` auto-selects, read the whole array only if
+# it is at most this many bytes; a larger raster is read tile by tile so a huge or
+# /vsicurl source is never materialised whole (#969). 256 MiB keeps the fast
+# whole-array path for everyday rasters while bounding the pathological case.
+_AUTO_TILE_BYTES = 256 * 1024 * 1024
+
 
 class Vectorize(_Engine["Dataset"]):
     """Mixin providing vectorization, clustering, and translate methods for Dataset."""
@@ -206,7 +212,7 @@ def to_feature_collection(
         self,
         mask: GeoDataFrame | None = None,
         add_geometry: str | None = None,
-        tile: bool = False,
+        tile: bool | None = None,
         tile_size: int = 256,
         touch: bool = True,
     ) -> DataFrame | GeoDataFrame:
@@ -229,12 +235,15 @@ def to_feature_collection(
             add_geometry (str):
                 "Polygon" or "Point" if you want to add a polygon geometry of the cells as column in dataframe.
                 Default is None.
-            tile (bool):
-                True to read the raster tile by tile rather than in one pass, which bounds
-                the peak allocation on a large raster. Default is False. The rows are the
-                same cells in the same row-major order either way -- `mask` included -- so
-                the choice is a memory/throughput trade, not a change of result, and
-                `add_geometry` is safe with either.
+            tile (bool | None):
+                Whether to read the raster tile by tile rather than in one pass, which
+                bounds the peak allocation on a large raster. `None` (default)
+                auto-selects: the whole array is read when it is at most ~256 MiB, else
+                the raster is tiled, so a huge or `/vsicurl` source is never
+                materialised whole. Pass `True`/`False` to force the choice. The rows
+                are the same cells in the same row-major order either way -- `mask`
+                included -- so it is a memory/throughput trade, not a change of result,
+                and `add_geometry` is safe with either.
             tile_size (int):
                 Tile size in cells, applied to both axes. Default is 256.
             touch (bool):
@@ -367,6 +376,15 @@ def to_feature_collection(
         else:
             src_ds = self._ds
 
+        # None auto-selects on the array's byte size: keep the fast whole-array read
+        # for everyday rasters, tile a large one so it is never materialised whole
+        # (#969). The tiled path is byte-identical (same row-major rows), so this only
+        # trades memory for throughput. An explicit True/False overrides.
+        if tile is None:
+            itemsize = np.dtype(src_ds.dtype[0]).itemsize
+            full_bytes = src_ds.rows * src_ds.columns * src_ds.band_count * itemsize
+            tile = full_bytes > _AUTO_TILE_BYTES
+
         # Both branches must read `src_ds` -- the cropped dataset when a mask was
         # given. Reading `self` here silently discarded the mask, so a tiled call
         # returned values for the whole raster while the geometry attached below
diff --git a/tests/dataset/analysis/test_dataset_vectorize.py b/tests/dataset/analysis/test_dataset_vectorize.py
index 1cb4ae003c..c35207b770 100644
--- a/tests/dataset/analysis/test_dataset_vectorize.py
+++ b/tests/dataset/analysis/test_dataset_vectorize.py
@@ -647,3 +647,51 @@ def test_no_data_cells_are_dropped_on_both_paths(self, uneven):
         assert not (tiled.to_numpy() == -9999.0).any(), (
             "no sentinel may survive into the frame"
         )
+
+    def test_auto_reads_the_whole_array_below_the_threshold(self, uneven, mocker):
+        """`tile=None` reads the whole array when it fits under the byte threshold.
+
+        Test scenario:
+            The 37x53 fixture is far below 256 MiB, so the default call takes the fast
+            whole-array path, not the tiled one.
+        """
+        full_spy = mocker.spy(Vectorize, "_extract_values_full")
+        tiled_spy = mocker.spy(Vectorize, "_extract_values_tiled")
+        uneven.to_feature_collection()
+        assert full_spy.call_count == 1, "small raster should read the whole array"
+        assert tiled_spy.call_count == 0, "small raster should not tile"
+
+    def test_auto_tiles_above_the_threshold_byte_identical(
+        self, uneven, mocker, monkeypatch
+    ):
+        """`tile=None` tiles a raster above the threshold, byte-identical to the full read.
+
+        Test scenario:
+            Drop the auto-tile threshold to 1 byte so the fixture exceeds it; the default
+            call then takes the tiled path and returns exactly the rows the untiled path
+            returns, in the same order.
+        """
+        monkeypatch.setattr("pyramids.dataset.engines.vectorize._AUTO_TILE_BYTES", 1)
+        tiled_spy = mocker.spy(Vectorize, "_extract_values_tiled")
+        auto = uneven.to_feature_collection(tile_size=8)
+        assert tiled_spy.call_count == 1, "large raster should tile"
+        untiled = uneven.to_feature_collection(tile=False)
+        pd.testing.assert_frame_equal(
+            auto.reset_index(drop=True), untiled.reset_index(drop=True)
+        )
+
+    def test_explicit_tile_false_overrides_the_auto_threshold(
+        self, uneven, mocker, monkeypatch
+    ):
+        """An explicit `tile=False` reads the whole array even above the threshold.
+
+        Test scenario:
+            With the threshold dropped to 1 byte (which would auto-tile), passing
+            `tile=False` still takes the whole-array path.
+        """
+        monkeypatch.setattr("pyramids.dataset.engines.vectorize._AUTO_TILE_BYTES", 1)
+        full_spy = mocker.spy(Vectorize, "_extract_values_full")
+        tiled_spy = mocker.spy(Vectorize, "_extract_values_tiled")
+        uneven.to_feature_collection(tile=False)
+        assert full_spy.call_count == 1, "explicit tile=False must read the whole array"
+        assert tiled_spy.call_count == 0, "explicit tile=False must not tile"

From f3c07a3256a0bceba36343c08a21ee7165f17bc4 Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 21:15:22 +0200
Subject: [PATCH 02/22] perf(analysis): seed sieve target via gdal.Translate,
 not a NumPy band copy

sieve() seeded its in-memory target with dst_band.WriteArray(src_band.ReadAsArray()),
materialising the whole source band as a NumPy array on top of the MEM buffer (peak
memory ~2x the band). Replace the manual MEM Create + NumPy round trip with a
gdal.Translate to a single-band MEM dataset: GDAL copies the band block by block in
the C layer, carrying the geotransform, CRS, dtype, and no-data across, so no NumPy
array is ever allocated. The gdal.SieveFilter call is unchanged; output is byte-identical.

Refs #969
---
 src/pyramids/dataset/engines/analysis.py | 15 ++++++---------
 1 file changed, 6 insertions(+), 9 deletions(-)

diff --git a/src/pyramids/dataset/engines/analysis.py b/src/pyramids/dataset/engines/analysis.py
index 1b8a26deca..8cfa02346e 100644
--- a/src/pyramids/dataset/engines/analysis.py
+++ b/src/pyramids/dataset/engines/analysis.py
@@ -1045,17 +1045,14 @@ def sieve(
             raise ValueError(f"connectedness must be 4 or 8, got {connectedness}.")
         validate_band_index(band, self._ds.band_count)
 
-        src_band = self._ds.raster.GetRasterBand(band + 1)
-        out_ds = gdal.GetDriverByName("MEM").Create(
-            "", self._ds.columns, self._ds.rows, 1, src_band.DataType
+        # Seed the sieve target with GDAL's block-based copy of the one band
+        # (geotransform, CRS, dtype, and no-data carried across in the C layer)
+        # instead of a full-band ``ReadAsArray`` -> ``WriteArray`` NumPy round
+        # trip, so the whole band is never materialised as a NumPy array (#969).
+        out_ds = gdal.Translate(
+            "", self._ds.raster, format="MEM", bandList=[band + 1]
         )
-        out_ds.SetGeoTransform(self._ds.geotransform)
-        out_ds.SetProjection(self._ds.crs)
         dst_band = out_ds.GetRasterBand(1)
-        dst_band.WriteArray(src_band.ReadAsArray())
-        no_data_value = src_band.GetNoDataValue()
-        if no_data_value is not None:
-            dst_band.SetNoDataValue(no_data_value)
 
         mask_band = mask.raster.GetRasterBand(1) if mask is not None else None
         gdal.SieveFilter(dst_band, mask_band, dst_band, threshold, connectedness)

From 9b23f4c21b8d523bff0d8f11f92e24fe38c56a92 Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 21:21:18 +0200
Subject: [PATCH 03/22] perf(spatial): tile the raster-mask crop instead of
 reading source+mask whole

_crop_aligned() read the full source (all bands) AND the full mask band into memory
before stamping the mask's no-data layout onto the source (peak ~ source + mask). For a
raster mask with no gap-filling, apply the mask tile by tile via a new
_crop_aligned_tiled helper: each 256x256 window reads its own source and mask block,
applies the no-data value, and writes straight into the destination, so the full arrays
are never held together. The mask-apply is purely per-pixel, so the result is
byte-identical to the whole-array path. The numpy-array mask (already in memory) and the
fill_gaps path (interpolates across the whole array) keep the eager path.

Refs #969
---
 src/pyramids/dataset/engines/spatial.py       | 77 +++++++++++++++----
 tests/dataset/spatial/test_dataset_spatial.py | 77 +++++++++++++++++++
 2 files changed, 141 insertions(+), 13 deletions(-)

diff --git a/src/pyramids/dataset/engines/spatial.py b/src/pyramids/dataset/engines/spatial.py
index 1955d7f158..be44ac701d 100644
--- a/src/pyramids/dataset/engines/spatial.py
+++ b/src/pyramids/dataset/engines/spatial.py
@@ -1327,15 +1327,11 @@ def _crop_aligned(
             row = mask.rows
             col = mask.columns
             mask_noval = mask.no_data_value[0]
-            # read_array() is called with no chunks=, so it always returns a plain
-            # ndarray here (the dask.Array arm of ArrayLike is unreachable).
-            mask_array = cast(np.typing.NDArray, mask.read_array(band=0))
         elif isinstance(mask, np.ndarray):
             if mask_noval is None:
                 raise ValueError(
                     "You have to enter the value of the no_val parameter when the mask is a numpy array"
                 )
-            mask_array = mask.copy()
             row, col = mask.shape
         else:
             raise TypeError(
@@ -1345,18 +1341,9 @@ def _crop_aligned(
 
         band_count = self._ds.band_count
         src_sref = sr_from_wkt(self._ds.crs)
-        # read_array() is called with no chunks=, so it always returns a plain
-        # ndarray here (the dask.Array arm of ArrayLike is unreachable).
-        src_array = cast(np.typing.NDArray, self._ds.read_array())
 
         self._assert_crop_aligned(mask, row, col)
 
-        mask_no_data = is_no_data(mask_array, mask_noval)
-        self._apply_mask_nodata(src_array, mask_no_data, band_count)
-
-        if fill_gaps:
-            src_array = self.fill_gaps(mask, src_array)
-
         dst = self._ds.__class__._create_dataset(
             col, row, band_count, self._ds.gdal_dtype[0], driver="MEM"
         )
@@ -1373,9 +1360,73 @@ def _crop_aligned(
         dst_obj = self._ds.__class__(dst)
         # set the no data value
         dst_obj._set_no_data_value(self._ds.no_data_value)
+
+        # Apply the mask's no-data layout tile by tile for the raster-mask,
+        # no-gap-fill case, so neither the full source (all bands) nor the full
+        # mask is ever materialised at once (#969). The mask-apply is purely
+        # per-pixel, so the tiled result is byte-identical to the whole-array
+        # path. The numpy-array mask (already in memory) and the gap-filling
+        # path (which interpolates across the whole array) stay eager below.
+        if isinstance(mask, RasterBase) and not fill_gaps:
+            self._crop_aligned_tiled(mask, mask_noval, dst_obj, band_count)
+            return dst_obj
+
+        # read_array() is called with no chunks=, so it always returns a plain
+        # ndarray here (the dask.Array arm of ArrayLike is unreachable).
+        if isinstance(mask, RasterBase):
+            mask_array = cast(np.typing.NDArray, mask.read_array(band=0))
+        else:
+            mask_array = mask.copy()
+        src_array = cast(np.typing.NDArray, self._ds.read_array())
+
+        mask_no_data = is_no_data(mask_array, mask_noval)
+        self._apply_mask_nodata(src_array, mask_no_data, band_count)
+
+        if fill_gaps:
+            src_array = self.fill_gaps(mask, src_array)
+
         self._write_bands(dst_obj, src_array, band_count)
         return dst_obj
 
+    def _crop_aligned_tiled(
+        self,
+        mask: RasterBase,
+        mask_noval: int | float | None,
+        dst_obj: Any,
+        band_count: int,
+    ) -> None:
+        """Stamp the mask's no-data layout onto the source one window at a time.
+
+        Reads the source (all bands) and the mask band a square tile at a time,
+        applies the source no-data value where the mask is no-data, and writes
+        each masked block straight into `dst_obj`, so the full arrays are never
+        held together. The operation is purely per-pixel, hence byte-identical
+        to the whole-array path (#969).
+
+        Args:
+            mask: The already-aligned raster mask supplying the no-data layout.
+            mask_noval: The mask's no-data value used to locate masked cells.
+            dst_obj: The destination Dataset the masked blocks are written into.
+            band_count: Number of bands in the source raster.
+        """
+        for xoff, yoff, xsize, ysize in self._ds.io._tile_offsets():
+            window = [xoff, yoff, xsize, ysize]
+            # read_array() is called with no chunks=, so it always returns a
+            # plain ndarray here (the dask.Array arm of ArrayLike is unreachable).
+            mask_tile = cast(
+                np.typing.NDArray, mask.read_array(band=0, window=window)
+            )
+            src_tile = cast(np.typing.NDArray, self._ds.read_array(window=window))
+            mask_no_data = is_no_data(mask_tile, mask_noval)
+            self._apply_mask_nodata(src_tile, mask_no_data, band_count)
+            if band_count > 1:
+                for band in range(band_count):
+                    dst_obj.raster.GetRasterBand(band + 1).WriteArray(
+                        src_tile[band, :, :], xoff, yoff
+                    )
+            else:
+                dst_obj.raster.GetRasterBand(1).WriteArray(src_tile, xoff, yoff)
+
     def _check_alignment(self, mask) -> bool:
         """Check if raster is aligned with a given mask raster."""
         if not isinstance(mask, RasterBase):
diff --git a/tests/dataset/spatial/test_dataset_spatial.py b/tests/dataset/spatial/test_dataset_spatial.py
index f2ac1a783d..2f272156d1 100644
--- a/tests/dataset/spatial/test_dataset_spatial.py
+++ b/tests/dataset/spatial/test_dataset_spatial.py
@@ -552,6 +552,83 @@ def test_crop_un_aligned(
         aligned_raster.spatial._crop_with_raster(mask_obj)
 
 
+class TestCropAlignedTiling:
+    """The tiled raster-mask path (`_crop_aligned_tiled`) is byte-identical across tile seams."""
+
+    @pytest.mark.parametrize("bands", [1, 2])
+    def test_tiled_matches_direct_apply_across_tile_boundaries(self, bands):
+        """A raster larger than one 256-px tile is masked exactly like a direct NumPy apply.
+
+        Args:
+            bands: Number of source bands (1 exercises the 2-D path, 2 the per-band path).
+
+        Test scenario:
+            No-data cells placed in three different 256-px tiles (including one on the column
+            seam) are stamped with the source no-data value identically to a whole-array apply.
+        """
+        rng = np.random.default_rng(0)
+        shape = (bands, 300, 300) if bands > 1 else (300, 300)
+        arr = (rng.random(shape) * 100).astype("float32")
+        src_no_data = -9999.0
+        geotransform = (0.0, 0.05, 0.0, 15.0, 0.0, -0.05)
+        src = Dataset.create_from_array(
+            arr, geo=geotransform, epsg=4326, no_data_value=src_no_data
+        )
+
+        mask_no_data = -1.0
+        mask_arr = np.ones((300, 300), dtype="float32")
+        mask_arr[10, 20] = mask_no_data
+        mask_arr[275, 290] = mask_no_data
+        mask_arr[128, 260] = mask_no_data
+        mask = Dataset.create_from_array(
+            mask_arr, geo=geotransform, epsg=4326, no_data_value=mask_no_data
+        )
+
+        cropped = src.spatial._crop_aligned(mask).read_array()
+
+        expected = arr.copy()
+        holes = mask_arr == mask_no_data
+        if bands > 1:
+            for band in range(bands):
+                expected[band][holes] = src_no_data
+        else:
+            expected[holes] = src_no_data
+        np.testing.assert_array_equal(
+            cropped,
+            expected,
+            err_msg="Tiled mask-apply must match the direct whole-array apply",
+        )
+
+    def test_tiled_and_fill_gaps_paths_agree(self):
+        """The tiled default and the eager fill_gaps=False result are identical on one raster.
+
+        Test scenario:
+            Running `_crop_aligned` (tiled) and forcing the eager whole-array branch via a
+            NumPy mask on the same source yield the same masked array.
+        """
+        rng = np.random.default_rng(1)
+        arr = (rng.random((260, 260)) * 10).astype("float64")
+        src_no_data = -32768.0
+        geotransform = (0.0, 0.05, 0.0, 13.0, 0.0, -0.05)
+        src = Dataset.create_from_array(
+            arr, geo=geotransform, epsg=4326, no_data_value=src_no_data
+        )
+        mask_arr = np.ones((260, 260), dtype="float64")
+        mask_arr[5, 5] = -1.0
+        mask_arr[259, 259] = -1.0
+        mask = Dataset.create_from_array(
+            mask_arr, geo=geotransform, epsg=4326, no_data_value=-1.0
+        )
+
+        tiled = src.spatial._crop_aligned(mask).read_array()
+        eager = src.spatial._crop_aligned(mask_arr, mask_noval=-1.0).read_array()
+        np.testing.assert_array_equal(
+            tiled,
+            eager,
+            err_msg="Tiled raster-mask path must match the eager numpy-mask path",
+        )
+
+
 class TestCropWithPolygon:
     def test_inplace(
         self,

From 25a4c1b881232014a95fbd57fdf222f5b500ef8a Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 21:25:12 +0200
Subject: [PATCH 04/22] perf(analysis): opt-in elementwise streaming for apply

apply() read the whole band and applied func to every domain value at once. Add an
opt-in elementwise=True that applies func one 256x256 tile at a time via a new
_apply_elementwise_tiled helper, so a very large or /vsicurl source is never
materialised whole. For a genuine per-pixel func the tiled result is byte-identical to
the whole-array pass; the shared domain-apply logic is factored into
_apply_func_to_domain so both paths behave identically. func that depend on the whole
array (normalisation, ranking, any global reduction) must keep the default
elementwise=False, whose behaviour is unchanged.

Refs #969
---
 src/pyramids/dataset/engines/analysis.py      |  80 ++++++++++---
 .../analysis/test_apply_elementwise.py        | 106 ++++++++++++++++++
 2 files changed, 173 insertions(+), 13 deletions(-)
 create mode 100644 tests/dataset/analysis/test_apply_elementwise.py

diff --git a/src/pyramids/dataset/engines/analysis.py b/src/pyramids/dataset/engines/analysis.py
index 8cfa02346e..8026d48edc 100644
--- a/src/pyramids/dataset/engines/analysis.py
+++ b/src/pyramids/dataset/engines/analysis.py
@@ -301,7 +301,14 @@ def _count(acc: int, strip: np.ndarray, _window: list[int]) -> int:
         domain_count = self._ds.rows * self._ds.columns - no_data_count
         return int(domain_count)
 
-    def apply(self, func, band: int = 0, inplace: bool = False) -> Dataset | None:
+    def apply(
+        self,
+        func,
+        band: int = 0,
+        inplace: bool = False,
+        *,
+        elementwise: bool = False,
+    ) -> Dataset | None:
         """Apply a function to all domain cells.
 
         - apply method executes a mathematical operation on the raster array.
@@ -315,6 +322,16 @@ def apply(self, func, band: int = 0, inplace: bool = False) -> Dataset | None:
             inplace (bool):
                 If True, the original dataset will be modified. If False, a new dataset will be created.
                 Default is False.
+            elementwise (bool):
+                Opt-in streaming mode. When `True`, `func` is applied one tile at
+                a time instead of to the whole band at once, so a very large or
+                `/vsicurl` source is never materialised whole. Only pass `True`
+                when `func` is a genuine **per-pixel** map (e.g. `np.abs`,
+                `lambda v: v * 2 + 1`): the tiled result is then byte-identical to
+                the default whole-array pass. A `func` that depends on the whole
+                array -- a min/max normalisation, a rank, any global reduction --
+                would give a different result tiled, so it must keep the default
+                `False`. Default `False` (whole-array, unchanged behaviour).
 
         Returns:
             Dataset | None:
@@ -362,19 +379,8 @@ def apply(self, func, band: int = 0, inplace: bool = False) -> Dataset | None:
             raise TypeError("The second argument should be a function")
 
         no_data_value = self._ds.no_data_value[band]
-        src_array = self._ds.read_array(band)
         dtype = self._ds.gdal_dtype[band]
 
-        new_array = np.full(
-            (self._ds.rows, self._ds.columns), no_data_value, dtype=src_array.dtype
-        )
-        domain_mask = inside_domain(src_array, no_data_value)
-        domain_values = src_array[domain_mask]
-        try:
-            new_array[domain_mask] = func(domain_values)
-        except (ValueError, TypeError):
-            new_array[domain_mask] = np.vectorize(func)(domain_values)
-
         dst_obj = self._ds.__class__._build_dataset(
             self._ds.columns,
             self._ds.rows,
@@ -384,13 +390,61 @@ def apply(self, func, band: int = 0, inplace: bool = False) -> Dataset | None:
             self._ds.crs,
             no_data_value,
         )
-        dst_obj.raster.GetRasterBand(1).WriteArray(new_array)
+        if elementwise:
+            self._apply_elementwise_tiled(func, band, no_data_value, dst_obj)
+        else:
+            src_array = self._ds.read_array(band)
+            new_array = np.full(
+                (self._ds.rows, self._ds.columns),
+                no_data_value,
+                dtype=src_array.dtype,
+            )
+            self._apply_func_to_domain(func, src_array, new_array, no_data_value)
+            dst_obj.raster.GetRasterBand(1).WriteArray(new_array)
 
         if inplace:
             self._ds._update_inplace(dst_obj.raster)
             return None
         return dst_obj
 
+    @staticmethod
+    def _apply_func_to_domain(func, src_array, out_array, no_data_value) -> None:
+        """Apply `func` to the domain (non-no-data) cells of `src_array` into `out_array`.
+
+        Args:
+            func: The per-domain-values callable to apply.
+            src_array: The source array supplying the domain values.
+            out_array: The pre-filled output array written in place.
+            no_data_value: The value marking cells to exclude from the domain.
+        """
+        domain_mask = inside_domain(src_array, no_data_value)
+        domain_values = src_array[domain_mask]
+        try:
+            out_array[domain_mask] = func(domain_values)
+        except (ValueError, TypeError):
+            out_array[domain_mask] = np.vectorize(func)(domain_values)
+
+    def _apply_elementwise_tiled(self, func, band, no_data_value, dst_obj) -> None:
+        """Apply an elementwise `func` over one band tile by tile, out of core.
+
+        Reads the band a square window at a time, applies `func` to that tile's
+        domain values, and writes the block straight into `dst_obj`, so the full
+        band is never materialised. For a per-pixel `func` the result is
+        byte-identical to the whole-array path (#969).
+
+        Args:
+            func: The per-pixel callable to apply to each tile's domain values.
+            band: Zero-based index of the source band to transform.
+            no_data_value: The source no-data value, preserved in excluded cells.
+            dst_obj: The single-band destination Dataset written in place.
+        """
+        dst_band = dst_obj.raster.GetRasterBand(1)
+        for xoff, yoff, xsize, ysize in self._ds.io._tile_offsets():
+            tile = self._ds.read_array(band=band, window=[xoff, yoff, xsize, ysize])
+            new_tile = np.full(tile.shape, no_data_value, dtype=tile.dtype)
+            self._apply_func_to_domain(func, tile, new_tile, no_data_value)
+            dst_band.WriteArray(new_tile, xoff, yoff)
+
     def fill(
         self, value: float | int, inplace: bool = False, path: str | Path | None = None
     ) -> Dataset | None:
diff --git a/tests/dataset/analysis/test_apply_elementwise.py b/tests/dataset/analysis/test_apply_elementwise.py
new file mode 100644
index 0000000000..a03229178f
--- /dev/null
+++ b/tests/dataset/analysis/test_apply_elementwise.py
@@ -0,0 +1,106 @@
+"""Tests for the opt-in streaming path of ``Analysis.apply`` (``elementwise=True``).
+
+The tiled path must be byte-identical to the default whole-array pass for a per-pixel
+``func``, including across 256-px tile seams, for both a plain band and a band carrying a
+no-data value, and it must honour band selection and the ``inplace`` flag.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from pyramids.dataset import Dataset
+
+pytestmark = pytest.mark.core
+
+
+class TestApplyElementwise:
+    """The ``elementwise=True`` tiled path matches the default whole-array apply."""
+
+    @pytest.mark.parametrize(
+        "func",
+        [np.abs, lambda v: v * 2 + 1, np.sqrt],
+        ids=["abs", "affine", "sqrt"],
+    )
+    def test_matches_whole_array_across_tile_boundaries(self, func):
+        """A raster larger than one 256-px tile transforms identically streamed and whole.
+
+        Args:
+            func: A per-pixel callable applied to the domain values.
+
+        Test scenario:
+            ``apply(elementwise=True)`` on a 300x300 raster equals ``apply()`` cell for cell.
+        """
+        arr = (np.random.default_rng(0).random((300, 300)) * 50).astype("float64")
+        ds = Dataset.create_from_array(
+            arr, top_left_corner=(0.0, 15.0), cell_size=0.05, epsg=4326
+        )
+        streamed = ds.apply(func, elementwise=True).read_array()
+        whole = ds.apply(func, elementwise=False).read_array()
+        np.testing.assert_array_equal(
+            streamed, whole, err_msg="Streamed apply must match the whole-array apply"
+        )
+
+    def test_preserves_no_data_cells(self):
+        """No-data cells stay no-data and only domain cells are transformed.
+
+        Test scenario:
+            A raster with scattered no-data in different tiles keeps those cells and doubles
+            the rest, matching the whole-array pass.
+        """
+        arr = np.ones((300, 300), dtype="float64") * 3.0
+        arr[10, 10] = -9999.0
+        arr[290, 295] = -9999.0
+        ds = Dataset.create_from_array(
+            arr,
+            top_left_corner=(0.0, 15.0),
+            cell_size=0.05,
+            epsg=4326,
+            no_data_value=-9999.0,
+        )
+        streamed = ds.apply(lambda v: v * 2, elementwise=True).read_array()
+        whole = ds.apply(lambda v: v * 2, elementwise=False).read_array()
+        np.testing.assert_array_equal(
+            streamed, whole, err_msg="No-data layout and domain map must match"
+        )
+        assert streamed[10, 10] == -9999.0, "No-data cell must be preserved"
+
+    def test_band_selection(self):
+        """A non-default band is streamed and returned as a single-band result.
+
+        Test scenario:
+            ``apply(elementwise=True, band=1)`` on a 2-band raster transforms band 1 and the
+            result matches the whole-array apply on the same band.
+        """
+        arr = np.stack(
+            [
+                np.zeros((300, 300), dtype="float64"),
+                (np.random.default_rng(1).random((300, 300)) * 10).astype("float64"),
+            ]
+        )
+        ds = Dataset.create_from_array(
+            arr, top_left_corner=(0.0, 15.0), cell_size=0.05, epsg=4326
+        )
+        streamed = ds.apply(np.abs, band=1, elementwise=True).read_array()
+        whole = ds.apply(np.abs, band=1, elementwise=False).read_array()
+        assert streamed.shape == (300, 300), f"Expected single band, got {streamed.shape}"
+        np.testing.assert_array_equal(
+            streamed, whole, err_msg="Band-1 streamed apply must match whole-array"
+        )
+
+    def test_inplace_streaming(self):
+        """``elementwise=True`` with ``inplace=True`` updates the source in place.
+
+        Test scenario:
+            The call returns ``self`` and the source array is the doubled result.
+        """
+        arr = (np.random.default_rng(2).random((300, 300)) * 4).astype("float64")
+        ds = Dataset.create_from_array(
+            arr, top_left_corner=(0.0, 15.0), cell_size=0.05, epsg=4326
+        )
+        result = ds.apply(lambda v: v * 2, inplace=True, elementwise=True)
+        assert result is ds, "inplace streaming apply should return self"
+        np.testing.assert_array_equal(
+            ds.read_array(), arr * 2, err_msg="In-place streamed result must be doubled"
+        )

From 420c2cdb02d8d45bf73b4f180b6924a52fcc054d Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 21:33:14 +0200
Subject: [PATCH 05/22] perf(analysis): opt-in decimated reads for
 plot_histogram and footprint

plot_histogram() flattened the whole band and footprint() read the whole band to build
the coverage mask, materialising a very large raster just to reduce it to a plot or a
polygon. Add an opt-in max_samples cap to both: when set and the band has more cells
than that, GDAL reads a nearest-neighbour decimated grid (~max_samples cells) in the C
layer instead of the full band. A shared _read_decimated helper does the read; footprint
scales the mask geotransform to the coarser grid via _scaled_geotransform so the polygon
coordinates stay geographically correct. The result is then approximate (a subsample /
coarser trace), documented as such. max_samples=None (default) reads every pixel, so both
methods stay exact and byte-identical to before.

Refs #969
---
 src/pyramids/dataset/engines/analysis.py      | 89 ++++++++++++++++++-
 .../analysis/test_dataset_vectorize.py        | 39 ++++++++
 tests/dataset/plot/test_plot_dataset.py       | 37 ++++++++
 3 files changed, 162 insertions(+), 3 deletions(-)

diff --git a/src/pyramids/dataset/engines/analysis.py b/src/pyramids/dataset/engines/analysis.py
index 8026d48edc..4382c0f014 100644
--- a/src/pyramids/dataset/engines/analysis.py
+++ b/src/pyramids/dataset/engines/analysis.py
@@ -1502,6 +1502,8 @@ def footprint(
         self,
         band: int = 0,
         exclude_values: list[Any] | None = None,
+        *,
+        max_samples: int | None = None,
     ) -> GeoDataFrame | None:
         """Extract the real coverage of the values in a certain band.
 
@@ -1522,6 +1524,16 @@ def footprint(
                 - This parameter is introduced particularly in the case of rasters that has the no_data_value stored in
                   the `no_data_value` property does not match the value stored in the band, so this option can correct
                   this behavior.
+            max_samples (int, optional):
+                Opt-in cap on how many pixels of the band are read to build the
+                coverage mask. When set and the band has more than
+                ``max_samples`` cells, GDAL reads a nearest-neighbour
+                **decimated** grid (~``max_samples`` cells) instead of the full
+                band, so a very large raster is footprinted without materialising
+                it whole. The extracted polygon is then **approximate** -- traced
+                on the coarser grid, so its edges and area are coarser than the
+                exact footprint. ``None`` (default) reads every pixel, so the
+                footprint is exact.
 
         Returns:
             GeoDataFrame:
@@ -1563,8 +1575,13 @@ def footprint(
 
               ```
         """
-        arr = self._ds.read_array(band=band)
+        arr = self._read_decimated(band, max_samples)
         no_data_val = self._ds.no_data_value[band]
+        # A decimated read spans the same extent with fewer, larger cells, so the
+        # mask's geotransform must scale its pixel size (and rotation terms) to
+        # the coarser grid; the origin is unchanged. Full-resolution reads leave
+        # the geotransform untouched.
+        geotransform = self._scaled_geotransform(arr.shape)
 
         self._warn_if_nodata_absent(arr, no_data_val)
         if exclude_values:
@@ -1587,7 +1604,7 @@ def footprint(
 
         new_dataset = Dataset.create_from_array(
             arr,
-            geo=self._ds.geotransform,
+            geo=geotransform,
             epsg=crs_spec(self._ds.epsg, self._ds.crs),
             no_data_value=0,
         )
@@ -1762,12 +1779,69 @@ def get_histogram(
         )
         return hist, ranges
 
+    def _read_decimated(self, band: int, max_samples: int | None) -> np.ndarray:
+        """Read a band whole, or a nearest-neighbour decimated version of it.
+
+        When `max_samples` is set and the band has more cells than that, GDAL
+        reads a coarser grid of roughly `max_samples` cells (decimated in the C
+        layer) so the whole band is never materialised; otherwise the full band
+        is read. Nearest-neighbour keeps the samples real pixel values.
+
+        Args:
+            band: Zero-based band index to read.
+            max_samples: Approximate pixel budget, or `None` for an exact read.
+
+        Returns:
+            np.ndarray: The band array, full-resolution or decimated.
+        """
+        rows = self._ds.rows
+        cols = self._ds.columns
+        total = rows * cols
+        if max_samples is None or total <= max_samples:
+            return cast(np.ndarray, self._ds.read_array(band=band))
+        factor = (total / max_samples) ** 0.5
+        out_rows = max(1, round(rows / factor))
+        out_cols = max(1, round(cols / factor))
+        return cast(
+            np.ndarray,
+            self._ds.read_array(
+                band=band, out_shape=(out_rows, out_cols), resampling="nearest"
+            ),
+        )
+
+    def _scaled_geotransform(
+        self, shape: tuple[int, ...]
+    ) -> tuple[float, float, float, float, float, float]:
+        """Geotransform for an array covering the source extent at `shape` cells.
+
+        A full-resolution `shape` returns the source geotransform unchanged; a
+        decimated `shape` (fewer/larger cells over the same extent) scales the
+        pixel-size and rotation terms by the row/column decimation factors while
+        keeping the origin fixed.
+
+        Args:
+            shape: The `(rows, cols)` of the (possibly decimated) array.
+
+        Returns:
+            tuple[float, float, float, float, float, float]: The six-element
+            geotransform for that grid.
+        """
+        d_rows, d_cols = shape
+        gt = self._ds.geotransform
+        if (d_rows, d_cols) == (self._ds.rows, self._ds.columns):
+            return (gt[0], gt[1], gt[2], gt[3], gt[4], gt[5])
+        sx = self._ds.columns / d_cols
+        sy = self._ds.rows / d_rows
+        return (gt[0], gt[1] * sx, gt[2] * sy, gt[3], gt[4] * sx, gt[5] * sy)
+
     def plot_histogram(
         self,
         band: int = 0,
         bins: int = 15,
         exclude_value: Any | None = None,
         ax: Any | None = None,
+        *,
+        max_samples: int | None = None,
         **kwargs: Any,
     ):
         """Plot the value distribution of a band as a histogram.
@@ -1788,6 +1862,15 @@ def plot_histogram(
                 band's no-data value and ``NaN``. Default is ``None``.
             ax (matplotlib.axes.Axes, optional):
                 Axes to draw on. A new figure/axes is created when ``None``.
+            max_samples (int, optional):
+                Opt-in cap on how many pixels are read. When set and the band
+                has more than ``max_samples`` cells, GDAL reads a
+                nearest-neighbour **decimated** version (~``max_samples`` cells)
+                instead of the full band, so a very large raster is histogrammed
+                without materialising it whole. The distribution is then
+                **approximate** -- a subsample of the pixels, the usual
+                expectation for a large raster. ``None`` (default) reads every
+                pixel, so the histogram is exact.
             **kwargs:
                 Style options forwarded to the ``HistogramGlyph``
                 constructor, filtered via
@@ -1830,7 +1913,7 @@ def plot_histogram(
         require_cleopatra()
         from cleopatra.glyphs.stats.histogram_glyph import HistogramGlyph
 
-        arr = self._ds.read_array(band=band).flatten()
+        arr = self._read_decimated(band, max_samples).flatten()
         no_data_value = self._ds.no_data_value[band]
         mask = np.ones(arr.shape, dtype=bool)
         if np.issubdtype(arr.dtype, np.floating):
diff --git a/tests/dataset/analysis/test_dataset_vectorize.py b/tests/dataset/analysis/test_dataset_vectorize.py
index c35207b770..74bc684579 100644
--- a/tests/dataset/analysis/test_dataset_vectorize.py
+++ b/tests/dataset/analysis/test_dataset_vectorize.py
@@ -419,6 +419,45 @@ def test_era5_one_band_no_no_data_value_in_raster(
         # the class should be 2
         assert next(iter(set(extent[dataset.band_names[0]]))) == 2
 
+    def test_max_samples_none_is_exact_full_read(self):
+        """The default `max_samples=None` traces the exact full-resolution footprint.
+
+        Test scenario:
+            A 200x200 raster with a solid covered block footprints to a polygon whose area
+            equals the block's geographic area.
+        """
+        arr = np.zeros((200, 200), dtype="float32")
+        arr[50:150, 60:160] = 5.0
+        ds = Dataset.create_from_array(
+            arr, top_left_corner=(0.0, 200.0), cell_size=1.0, epsg=3857
+        )
+        extent = ds.footprint(exclude_values=[0])
+        assert extent is not None, "A covered block must yield a footprint"
+        assert float(extent.geometry.area.sum()) == pytest.approx(100 * 100), (
+            "Exact footprint area must equal the covered block area"
+        )
+
+    def test_max_samples_decimates_and_keeps_extent(self):
+        """`max_samples` traces an approximate footprint on a coarser grid over the same extent.
+
+        Test scenario:
+            The decimated footprint of a solid block stays inside the raster bounds and
+            covers roughly the block's area (coarser, not exact).
+        """
+        arr = np.zeros((200, 200), dtype="float32")
+        arr[50:150, 60:160] = 5.0
+        ds = Dataset.create_from_array(
+            arr, top_left_corner=(0.0, 200.0), cell_size=1.0, epsg=3857
+        )
+        approx = ds.footprint(exclude_values=[0], max_samples=400)
+        assert approx is not None, "Decimated footprint must still be produced"
+        minx, miny, maxx, maxy = approx.total_bounds
+        assert minx >= -1e-6 and miny >= -1e-6, "Footprint must stay within raster bounds"
+        assert maxx <= 200 + 1e-6 and maxy <= 200 + 1e-6, "Footprint must stay within bounds"
+        assert float(approx.geometry.area.sum()) == pytest.approx(10000, rel=0.5), (
+            "Decimated footprint area should be roughly the block area"
+        )
+
 
 class TestToFeatureCollectionMaskTiling:
     """The mask must be honoured on both the tiled and the non-tiled path."""
diff --git a/tests/dataset/plot/test_plot_dataset.py b/tests/dataset/plot/test_plot_dataset.py
index 485838d740..088635a762 100644
--- a/tests/dataset/plot/test_plot_dataset.py
+++ b/tests/dataset/plot/test_plot_dataset.py
@@ -139,6 +139,43 @@ def histogram(self, bins=15):
             3.0,
         ], f"nodata (-9999) and exclude_value (7.0) must be dropped; got {vals}"
 
+    @pytest.mark.plot
+    def test_plot_histogram_max_samples_decimates(self):
+        """``max_samples`` histograms a decimated subsample instead of every pixel.
+
+        Test scenario:
+            A 100x100 raster (10,000 cells) plotted with ``max_samples=400`` hands the glyph
+            far fewer samples than the full band, while the default reads all 10,000.
+        """
+        arr = np.arange(100 * 100, dtype="float32").reshape(100, 100)
+        dataset = Dataset.create_from_array(
+            arr, top_left_corner=(0, 0), cell_size=1.0, epsg=4326
+        )
+        captured: dict = {}
+
+        class _FakeSG:
+            @staticmethod
+            def filter_kwargs(kw):
+                return {}
+
+            def __init__(self, values, ax=None, **kwargs):
+                captured["n"] = np.asarray(values).size
+
+            def histogram(self, bins=15):
+                return ("fig", "ax", {})
+
+        with patch(
+            "cleopatra.glyphs.stats.histogram_glyph.HistogramGlyph", new=_FakeSG
+        ):
+            dataset.plot_histogram(band=0, max_samples=400)
+            decimated = captured["n"]
+            dataset.plot_histogram(band=0)
+            full = captured["n"]
+        assert full == 10000, f"Exact read must see every pixel, got {full}"
+        assert decimated <= 500, (
+            f"max_samples=400 must decimate to ~400 samples, got {decimated}"
+        )
+
     @pytest.mark.plot
     def test_invalid_color_scale_raises(self, src: Dataset):
         """A loose ``color_scale`` kwarg is rejected with a clear ``ValueError``.

From 3459f2fc4027507413a33e0003e4c77c6dbb187e Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 21:49:49 +0200
Subject: [PATCH 06/22] perf(bands): stream the change_no_data_value swap and
 add a disk output

change_no_data_value read each band whole (a full NumPy array per band) to swap the old
no-data value for the new one. Keep GDAL's block-based CreateCopy (so every band's colour
table, description, scale/offset, RAT and metadata survive exactly) but stream the swap
one 256x256 tile at a time via a new _swap_no_data_tiled helper, so a full band is never
held as a NumPy array. Add a path= option: when given, the clone is a disk-backed GeoTIFF
and the tiled swap makes the whole operation out-of-core; FlushCache persists it before
reopen. Output is byte-identical to before.

Because NetCDF.read_array takes variable as its first positional, the old positional
read_array(band) call incidentally raised a "pinned" ValueError on a variable-view
NetCDF; the windowed read (band= keyword) removes that accident, so change_no_data_value
now works on a variable view like fill() and the other band ops already do. The test that
asserted the accidental guard is updated to assert it works.

Refs #969
---
 src/pyramids/dataset/engines/bands.py         | 82 ++++++++++++++---
 tests/dataset/unit/test_unit_nodata.py        | 88 ++++++++++++++++++-
 .../netcdf/samples/test_inherited_arg_ops.py  | 16 +++-
 3 files changed, 167 insertions(+), 19 deletions(-)

diff --git a/src/pyramids/dataset/engines/bands.py b/src/pyramids/dataset/engines/bands.py
index 52f6448151..4756b5b863 100644
--- a/src/pyramids/dataset/engines/bands.py
+++ b/src/pyramids/dataset/engines/bands.py
@@ -9,6 +9,7 @@
 
 import warnings
 from collections.abc import Iterable
+from pathlib import Path
 from typing import TYPE_CHECKING, Any, cast
 
 import geopandas as gpd
@@ -1349,7 +1350,12 @@ def _change_no_data_value_attr(self, band: int, no_data_value) -> None:
         self._ds._no_data_value[band] = no_data_value
 
     def change_no_data_value(
-        self, new_value: Any, old_value: Any | None = None, inplace: bool = False
+        self,
+        new_value: Any,
+        old_value: Any | None = None,
+        inplace: bool = False,
+        *,
+        path: str | Path | None = None,
     ) -> Dataset | None:
         """Change No Data Value.
             - Set the no data value in all raster bands.
@@ -1366,6 +1372,11 @@ def change_no_data_value(
             inplace (bool):
                 If True, the original dataset will be modified. If False, a new dataset will be created.
                 Default is False.
+            path (str | Path | None):
+                Output `.tif` path for a disk-backed result. When given, the raster
+                is cloned to that GeoTIFF and the no-data swap is streamed tile by
+                tile, so the whole raster is never held in RAM (genuinely
+                out-of-core). `None` (default) keeps the result in memory.
 
         Returns:
             Dataset | None:
@@ -1381,7 +1392,8 @@ def change_no_data_value(
                 match `band_count`.
 
         Warning:
-            The `change_no_data_value` method creates a new dataset in memory in order to change the `no_data_value` in the raster bands.
+            With `path=None` the method clones the raster in memory to change the
+            `no_data_value`; pass `path` for a disk-backed, out-of-core result.
         Examples:
             - Create a Dataset (4 bands, 10 rows, 10 columns) at lon/lat (0, 0):
               ```python
@@ -1426,8 +1438,16 @@ def change_no_data_value(
                 f"old_value must be a scalar or a list of length band_count "
                 f"({self._ds.band_count}); got a list of length {len(old_value)}."
             )
-        dst = gdal.GetDriverByName("MEM").CreateCopy("", self._ds.raster, 0)
-        # create a new dataset
+        # Clone the full header + pixels with GDAL's block-based CreateCopy so
+        # every band's colour table, description, scale/offset, RAT and metadata
+        # survive exactly (no explicit per-property copy to drift out of sync).
+        # A `path` makes the clone a disk-backed GeoTIFF, so with the tiled swap
+        # below the whole raster is never held in RAM (out-of-core); `None` keeps
+        # it in memory. The old<->new no-data swap is then streamed one tile at a
+        # time, so a full band is never materialised as a NumPy array (#969).
+        driver = "GTiff" if path is not None else "MEM"
+        target = str(path) if path is not None else ""
+        dst = gdal.GetDriverByName(driver).CreateCopy(target, self._ds.raster, 0)
         new_dataset = self._ds.__class__(dst, "write")
         # the new_value could change inside the _set_no_data_value method before it is used to set the no_data_value
         # attribute in the gdal object/pyramids object and to fill the band.
@@ -1436,14 +1456,55 @@ def change_no_data_value(
         # updated.
         new_value = new_dataset.no_data_value
         for band in range(self._ds.band_count):
-            arr = self._ds.read_array(band)
             # old_value is normalized to a per-band list above (matching
             # new_value); index it per-band here too instead of comparing
             # against the whole list.
             band_old_value = old_value[band] if old_value is not None else None
+            self._swap_no_data_tiled(
+                new_dataset, band, band_old_value, new_value[band]
+            )
+        # Flush the block cache so a disk-backed GeoTIFF has the swapped pixels on
+        # disk before it is reopened (a no-op for the in-memory driver).
+        new_dataset.raster.FlushCache()
+
+        if inplace:
+            self._ds._update_inplace(new_dataset.raster)
+            return None
+        return new_dataset
+
+    def _swap_no_data_tiled(
+        self,
+        new_dataset: Dataset,
+        band: int,
+        band_old_value: Any,
+        new_band_value: Any,
+    ) -> None:
+        """Replace a band's old-no-data cells with the new value, one tile at a time.
+
+        The caller's `_set_no_data_value` fills the destination band with the new
+        no-data value, so every tile is read from the *source*, has its old
+        no-data cells swapped to the new value, and is written back at its offset
+        -- reconstructing the band exactly as the previous whole-band read/swap/
+        write did, but never holding more than one tile as a NumPy array. The swap
+        is attempted on every tile so an invalid `new_band_value` dtype raises just
+        as the whole-band assignment did, even for a band with no matching cells.
+
+        Args:
+            new_dataset: The destination Dataset written in place.
+            band: Zero-based index of the band to swap.
+            band_old_value: The old no-data value to locate (``None`` matches NaN).
+            new_band_value: The new no-data value written into the matched cells.
+
+        Raises:
+            NoDataValueError: `new_band_value` cannot be stored in the band dtype.
+        """
+        dst_band = new_dataset.raster.GetRasterBand(band + 1)
+        for xoff, yoff, xsize, ysize in self._ds.io._tile_offsets():
+            tile = self._ds.read_array(band=band, window=[xoff, yoff, xsize, ysize])
+            mask = is_no_data(tile, band_old_value)
             try:
                 with np.errstate(invalid="raise"):
-                    arr[is_no_data(arr, band_old_value)] = new_value[band]
+                    tile[mask] = new_band_value
             # A dtype mismatch surfaces differently across numpy paths: a None value
             # is not subscriptable (TypeError), a NaN cast into an integer band raises
             # ValueError ("cannot convert float NaN to integer"), and an invalid
@@ -1451,12 +1512,7 @@ def change_no_data_value(
             # to the package-level NoDataValueError.
             except (TypeError, ValueError, FloatingPointError):
                 raise NoDataValueError(
-                    f"The dtype of the given no_data_value: {new_value[band]} differs from the dtype of the "
+                    f"The dtype of the given no_data_value: {new_band_value} differs from the dtype of the "
                     f"band: {gdal_to_numpy_dtype(self._ds.gdal_dtype[band])}"
                 )
-            new_dataset.raster.GetRasterBand(band + 1).WriteArray(arr)
-
-        if inplace:
-            self._ds._update_inplace(new_dataset.raster)
-            return None
-        return new_dataset
+            dst_band.WriteArray(tile, xoff, yoff)
diff --git a/tests/dataset/unit/test_unit_nodata.py b/tests/dataset/unit/test_unit_nodata.py
index 6a63fd64ec..cfac648aab 100644
--- a/tests/dataset/unit/test_unit_nodata.py
+++ b/tests/dataset/unit/test_unit_nodata.py
@@ -626,9 +626,9 @@ def test_change_nodata_type_error_raises(self):
         )
         original_read = ds.read_array
 
-        def mock_read(band=None):
+        def mock_read(band=None, window=None, **kwargs):
             """Return array that raises TypeError on assignment."""
-            result = original_read(band=band)
+            result = original_read(band=band, window=window, **kwargs)
             mock_arr = MagicMock(wraps=result)
 
             def raise_type_error(key, value):
@@ -643,6 +643,90 @@ def raise_type_error(key, value):
                 ds.change_no_data_value(-1.0, old_value=-9999.0)
 
 
+class TestChangeNoDataValueStreaming:
+    """The tiled/CreateCopy redesign of change_no_data_value stays correct and lossless."""
+
+    def test_multi_tile_matches_direct_swap(self):
+        """A raster larger than one 256-px tile swaps no-data identically to a direct swap.
+
+        Test scenario:
+            Old no-data cells placed in three different tiles are all rewritten to the new
+            value and every other cell is preserved, matching a whole-array NumPy swap.
+        """
+        arr = (np.random.default_rng(0).random((300, 300)) * 10).astype("float32")
+        arr[10, 10] = -9999.0
+        arr[290, 295] = -9999.0
+        arr[130, 260] = -9999.0
+        ds = Dataset.create_from_array(
+            arr,
+            top_left_corner=(0.0, 300.0),
+            cell_size=0.05,
+            epsg=4326,
+            no_data_value=-9999.0,
+        )
+        result = ds.change_no_data_value(-1.0, old_value=-9999.0).read_array()
+        expected = arr.copy()
+        expected[expected == -9999.0] = -1.0
+        np.testing.assert_array_equal(
+            result, expected, err_msg="Tiled no-data swap must match the direct swap"
+        )
+        assert result.dtype == np.float32, "Band dtype must be preserved"
+
+    def test_preserves_color_table_and_scale_offset(self):
+        """CreateCopy carries the colour table, band description, and scale/offset across.
+
+        Test scenario:
+            A byte raster with a colour table, a band description, and a scale/offset keeps
+            all three after change_no_data_value.
+        """
+        from osgeo import gdal
+
+        mem = gdal.GetDriverByName("MEM").Create("", 4, 4, 1, gdal.GDT_Byte)
+        mem.SetGeoTransform((0.0, 1.0, 0.0, 4.0, 0.0, -1.0))
+        band = mem.GetRasterBand(1)
+        band.WriteArray(np.array([[0, 1, 2, 3]] * 4, dtype="uint8"))
+        band.SetNoDataValue(0)
+        band.SetDescription("classes")
+        band.SetScale(0.5)
+        band.SetOffset(2.0)
+        ct = gdal.ColorTable()
+        ct.SetColorEntry(1, (255, 0, 0, 255))
+        band.SetColorTable(ct)
+        ds = Dataset(mem)
+
+        result = ds.change_no_data_value(255, old_value=0)
+        out_band = result.raster.GetRasterBand(1)
+        assert out_band.GetDescription() == "classes", "Band description must survive"
+        assert out_band.GetScale() == 0.5, "Band scale must survive"
+        assert out_band.GetOffset() == 2.0, "Band offset must survive"
+        assert out_band.GetColorTable() is not None, "Colour table must survive"
+        assert out_band.GetColorTable().GetColorEntry(1) == (255, 0, 0, 255), (
+            "Colour table entries must survive"
+        )
+
+    def test_path_writes_disk_backed_result(self, tmp_path):
+        """`path=` writes a disk-backed GeoTIFF with the new no-data value.
+
+        Test scenario:
+            Passing a `.tif` path produces a file on disk whose band carries the new
+            no-data value and the swapped cells.
+        """
+        arr = np.array([[1.0, 2.0], [-9999.0, 4.0]], dtype="float32")
+        ds = Dataset.create_from_array(
+            arr,
+            top_left_corner=(0.0, 2.0),
+            cell_size=1.0,
+            epsg=4326,
+            no_data_value=-9999.0,
+        )
+        out = tmp_path / "changed.tif"
+        result = ds.change_no_data_value(-1.0, old_value=-9999.0, path=out)
+        assert out.exists(), "path= must write a file to disk"
+        assert result.no_data_value[0] == -1.0, "Disk result must carry the new no-data"
+        reopened = Dataset.read_file(str(out)).read_array()
+        assert reopened[1, 0] == -1.0, "Swapped cell must be persisted on disk"
+
+
 class TestChangeNoDataAttrConversion:
     """Tests for _change_no_data_value_attr type conversion."""
 
diff --git a/tests/netcdf/samples/test_inherited_arg_ops.py b/tests/netcdf/samples/test_inherited_arg_ops.py
index 81c4d9f121..f8c390eb65 100644
--- a/tests/netcdf/samples/test_inherited_arg_ops.py
+++ b/tests/netcdf/samples/test_inherited_arg_ops.py
@@ -121,10 +121,18 @@ def test_get_band_by_color_absent_returns_none(tos):
     assert tos.get_band_by_color("gray_index") is None
 
 
-def test_change_no_data_value_guarded_on_variable_view(tos):
-    """change_no_data_value is guarded on a variable-pinned view (clear error, not a crash)."""
-    with pytest.raises(ValueError, match="pinned"):
-        tos.change_no_data_value(-999.0, (tos.no_data_value or [None])[0])
+def test_change_no_data_value_on_variable_view(tos):
+    """change_no_data_value works on a variable-pinned view, like the other band ops.
+
+    It clones the view, swaps every band's old no-data cells to the new value, and returns
+    a fresh dataset with the updated no-data value -- the source view is left untouched.
+    """
+    old = (tos.no_data_value or [None])[0]
+    result = tos.change_no_data_value(-999.0, old)
+    assert result is not None, "change_no_data_value must return a dataset for a view"
+    assert all(nd == -999.0 for nd in result.no_data_value), (
+        f"every band's no-data must be -999.0, got {result.no_data_value}"
+    )
 
 
 def _classes(v):

From 54b018ee01d678b183adb165dd48cc72adb612d6 Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 22:24:52 +0200
Subject: [PATCH 07/22] style: apply black line-wrapping to the #969 streaming
 changes

---
 src/pyramids/dataset/engines/analysis.py         | 4 +---
 src/pyramids/dataset/engines/bands.py            | 4 +---
 src/pyramids/dataset/engines/spatial.py          | 4 +---
 tests/dataset/analysis/test_apply_elementwise.py | 4 +++-
 tests/dataset/analysis/test_dataset_vectorize.py | 8 ++++++--
 5 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/pyramids/dataset/engines/analysis.py b/src/pyramids/dataset/engines/analysis.py
index 4382c0f014..118c68f9fd 100644
--- a/src/pyramids/dataset/engines/analysis.py
+++ b/src/pyramids/dataset/engines/analysis.py
@@ -1103,9 +1103,7 @@ def sieve(
         # (geotransform, CRS, dtype, and no-data carried across in the C layer)
         # instead of a full-band ``ReadAsArray`` -> ``WriteArray`` NumPy round
         # trip, so the whole band is never materialised as a NumPy array (#969).
-        out_ds = gdal.Translate(
-            "", self._ds.raster, format="MEM", bandList=[band + 1]
-        )
+        out_ds = gdal.Translate("", self._ds.raster, format="MEM", bandList=[band + 1])
         dst_band = out_ds.GetRasterBand(1)
 
         mask_band = mask.raster.GetRasterBand(1) if mask is not None else None
diff --git a/src/pyramids/dataset/engines/bands.py b/src/pyramids/dataset/engines/bands.py
index 4756b5b863..d9a3d7b46c 100644
--- a/src/pyramids/dataset/engines/bands.py
+++ b/src/pyramids/dataset/engines/bands.py
@@ -1460,9 +1460,7 @@ def change_no_data_value(
             # new_value); index it per-band here too instead of comparing
             # against the whole list.
             band_old_value = old_value[band] if old_value is not None else None
-            self._swap_no_data_tiled(
-                new_dataset, band, band_old_value, new_value[band]
-            )
+            self._swap_no_data_tiled(new_dataset, band, band_old_value, new_value[band])
         # Flush the block cache so a disk-backed GeoTIFF has the swapped pixels on
         # disk before it is reopened (a no-op for the in-memory driver).
         new_dataset.raster.FlushCache()
diff --git a/src/pyramids/dataset/engines/spatial.py b/src/pyramids/dataset/engines/spatial.py
index be44ac701d..b1f2964379 100644
--- a/src/pyramids/dataset/engines/spatial.py
+++ b/src/pyramids/dataset/engines/spatial.py
@@ -1413,9 +1413,7 @@ def _crop_aligned_tiled(
             window = [xoff, yoff, xsize, ysize]
             # read_array() is called with no chunks=, so it always returns a
             # plain ndarray here (the dask.Array arm of ArrayLike is unreachable).
-            mask_tile = cast(
-                np.typing.NDArray, mask.read_array(band=0, window=window)
-            )
+            mask_tile = cast(np.typing.NDArray, mask.read_array(band=0, window=window))
             src_tile = cast(np.typing.NDArray, self._ds.read_array(window=window))
             mask_no_data = is_no_data(mask_tile, mask_noval)
             self._apply_mask_nodata(src_tile, mask_no_data, band_count)
diff --git a/tests/dataset/analysis/test_apply_elementwise.py b/tests/dataset/analysis/test_apply_elementwise.py
index a03229178f..8ca41233a5 100644
--- a/tests/dataset/analysis/test_apply_elementwise.py
+++ b/tests/dataset/analysis/test_apply_elementwise.py
@@ -84,7 +84,9 @@ def test_band_selection(self):
         )
         streamed = ds.apply(np.abs, band=1, elementwise=True).read_array()
         whole = ds.apply(np.abs, band=1, elementwise=False).read_array()
-        assert streamed.shape == (300, 300), f"Expected single band, got {streamed.shape}"
+        assert streamed.shape == (300, 300), (
+            f"Expected single band, got {streamed.shape}"
+        )
         np.testing.assert_array_equal(
             streamed, whole, err_msg="Band-1 streamed apply must match whole-array"
         )
diff --git a/tests/dataset/analysis/test_dataset_vectorize.py b/tests/dataset/analysis/test_dataset_vectorize.py
index 74bc684579..fc81a887db 100644
--- a/tests/dataset/analysis/test_dataset_vectorize.py
+++ b/tests/dataset/analysis/test_dataset_vectorize.py
@@ -452,8 +452,12 @@ def test_max_samples_decimates_and_keeps_extent(self):
         approx = ds.footprint(exclude_values=[0], max_samples=400)
         assert approx is not None, "Decimated footprint must still be produced"
         minx, miny, maxx, maxy = approx.total_bounds
-        assert minx >= -1e-6 and miny >= -1e-6, "Footprint must stay within raster bounds"
-        assert maxx <= 200 + 1e-6 and maxy <= 200 + 1e-6, "Footprint must stay within bounds"
+        assert minx >= -1e-6 and miny >= -1e-6, (
+            "Footprint must stay within raster bounds"
+        )
+        assert maxx <= 200 + 1e-6 and maxy <= 200 + 1e-6, (
+            "Footprint must stay within bounds"
+        )
         assert float(approx.geometry.area.sum()) == pytest.approx(10000, rel=0.5), (
             "Decimated footprint area should be roughly the block area"
         )

From dca6e63a821eb1ad17e6e3e6a7a90f09de4aa5d4 Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 22:46:23 +0200
Subject: [PATCH 08/22] fix(analysis): read band by keyword in apply's default
 path so it works on a NetCDF view

apply()'s default path called read_array(band) positionally; because NetCDF.read_array
puts variable first, that mis-bound band->variable and raised a spurious "pinned" error
on a variable view, while apply(elementwise=True) (which passes band=) worked. Both paths
now read with band= as a keyword, so apply behaves consistently on a variable view like
the other band ops. The test that asserted the accidental guard now asserts both paths
work and agree.

Refs #969
---
 src/pyramids/dataset/engines/analysis.py      |  4 +++-
 .../netcdf/samples/test_inherited_arg_ops.py  | 19 ++++++++++++++++---
 2 files changed, 19 insertions(+), 4 deletions(-)

diff --git a/src/pyramids/dataset/engines/analysis.py b/src/pyramids/dataset/engines/analysis.py
index 118c68f9fd..d633fcf748 100644
--- a/src/pyramids/dataset/engines/analysis.py
+++ b/src/pyramids/dataset/engines/analysis.py
@@ -393,7 +393,9 @@ def apply(
         if elementwise:
             self._apply_elementwise_tiled(func, band, no_data_value, dst_obj)
         else:
-            src_array = self._ds.read_array(band)
+            # `band=` as a keyword, never positional: NetCDF.read_array puts
+            # `variable` first, so read_array(band) mis-binds on a variable view.
+            src_array = self._ds.read_array(band=band)
             new_array = np.full(
                 (self._ds.rows, self._ds.columns),
                 no_data_value,
diff --git a/tests/netcdf/samples/test_inherited_arg_ops.py b/tests/netcdf/samples/test_inherited_arg_ops.py
index f8c390eb65..d3a6136313 100644
--- a/tests/netcdf/samples/test_inherited_arg_ops.py
+++ b/tests/netcdf/samples/test_inherited_arg_ops.py
@@ -161,9 +161,22 @@ def test_set_attribute_table(tos):
     assert df_out is not None and len(df_out) == len(df_in)
 
 
-def test_apply_guarded_on_variable_view(tos):
-    with pytest.raises(ValueError, match="pinned"):
-        tos.apply(lambda a: a + 1)
+def test_apply_on_variable_view(tos):
+    """apply works on a variable-pinned view (both the default and elementwise paths).
+
+    The old "pinned" error was an accidental positional read_array(band) mis-bind, not a
+    real guard; both paths now read with band= as a keyword and return a result.
+    """
+    default = tos.apply(lambda a: a + 1)
+    streamed = tos.apply(lambda a: a + 1, elementwise=True)
+    assert default is not None and streamed is not None, (
+        "apply must return a dataset on a variable view for both paths"
+    )
+    np.testing.assert_array_equal(
+        np.asarray(default.read_array(band=0)),
+        np.asarray(streamed.read_array(band=0)),
+        err_msg="default and elementwise apply must agree on a variable view",
+    )
 
 
 def test_set_rpcs_guarded_read_only(tos):

From 3bba72a514820429eed1920a33d93112957e7d16 Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 22:52:26 +0200
Subject: [PATCH 09/22] fix(bands): remove the partial GeoTIFF if
 change_no_data_value(path=) fails mid-swap

If the tiled no-data swap raised mid-stream (e.g. a dtype-mismatch NoDataValueError) on the
disk-backed path, a half-written GeoTIFF was left at path with the write handle still open.
Wrap the swap in try/except that, on the disk path, closes the wrapper handle, drops the
local dst reference (both are needed or GDAL keeps the file locked on Windows), and deletes
the partial file plus its .aux.xml sidecar before re-raising.

Refs #969
---
 src/pyramids/dataset/engines/bands.py | 39 ++++++++++++++++++---------
 1 file changed, 27 insertions(+), 12 deletions(-)

diff --git a/src/pyramids/dataset/engines/bands.py b/src/pyramids/dataset/engines/bands.py
index d9a3d7b46c..897c76daf4 100644
--- a/src/pyramids/dataset/engines/bands.py
+++ b/src/pyramids/dataset/engines/bands.py
@@ -1449,18 +1449,33 @@ def change_no_data_value(
         target = str(path) if path is not None else ""
         dst = gdal.GetDriverByName(driver).CreateCopy(target, self._ds.raster, 0)
         new_dataset = self._ds.__class__(dst, "write")
-        # the new_value could change inside the _set_no_data_value method before it is used to set the no_data_value
-        # attribute in the gdal object/pyramids object and to fill the band.
-        new_dataset._set_no_data_value(new_value)
-        # now we have to use the no_data_value value in the no_data_value attribute in the Dataset object as it is
-        # updated.
-        new_value = new_dataset.no_data_value
-        for band in range(self._ds.band_count):
-            # old_value is normalized to a per-band list above (matching
-            # new_value); index it per-band here too instead of comparing
-            # against the whole list.
-            band_old_value = old_value[band] if old_value is not None else None
-            self._swap_no_data_tiled(new_dataset, band, band_old_value, new_value[band])
+        try:
+            # the new_value could change inside the _set_no_data_value method before it is used to set the no_data_value
+            # attribute in the gdal object/pyramids object and to fill the band.
+            new_dataset._set_no_data_value(new_value)
+            # now we have to use the no_data_value value in the no_data_value attribute in the Dataset object as it is
+            # updated.
+            new_value = new_dataset.no_data_value
+            for band in range(self._ds.band_count):
+                # old_value is normalized to a per-band list above (matching
+                # new_value); index it per-band here too instead of comparing
+                # against the whole list.
+                band_old_value = old_value[band] if old_value is not None else None
+                self._swap_no_data_tiled(
+                    new_dataset, band, band_old_value, new_value[band]
+                )
+        except Exception:
+            # A mid-stream failure (e.g. a dtype-mismatch NoDataValueError) must not
+            # leave a half-written GeoTIFF behind on the disk path: release every
+            # handle to the file (both the wrapper and the local `dst` reference,
+            # or GDAL keeps the file locked on Windows) and delete the partial file
+            # and its sidecar before re-raising.
+            if path is not None:
+                new_dataset.close()
+                dst = None
+                Path(target).unlink(missing_ok=True)
+                Path(f"{target}.aux.xml").unlink(missing_ok=True)
+            raise
         # Flush the block cache so a disk-backed GeoTIFF has the swapped pixels on
         # disk before it is reopened (a no-op for the in-memory driver).
         new_dataset.raster.FlushCache()

From 124d79491140ba4fbced3c17247fa4bf244ef24f Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 22:52:27 +0200
Subject: [PATCH 10/22] test(bands): pin change_no_data_value streaming
 guarantees

Add regression tests for the tiled/CreateCopy redesign: a per-band old_value list on a
multi-band raster, the dtype guard firing on a band with zero matching cells (proving the
swap is attempted per tile), RAT preservation through CreateCopy, inplace=True + path=
updating the source and persisting to disk, and a mid-swap failure removing the partial
disk file.

Refs #969
---
 tests/dataset/unit/test_unit_nodata.py | 138 +++++++++++++++++++++++++
 1 file changed, 138 insertions(+)

diff --git a/tests/dataset/unit/test_unit_nodata.py b/tests/dataset/unit/test_unit_nodata.py
index cfac648aab..6d2d8014c7 100644
--- a/tests/dataset/unit/test_unit_nodata.py
+++ b/tests/dataset/unit/test_unit_nodata.py
@@ -726,6 +726,144 @@ def test_path_writes_disk_backed_result(self, tmp_path):
         reopened = Dataset.read_file(str(out)).read_array()
         assert reopened[1, 0] == -1.0, "Swapped cell must be persisted on disk"
 
+    def test_multi_band_per_band_old_value_list(self):
+        """A per-band old_value list swaps each band's own no-data cells.
+
+        Test scenario:
+            A 2-band raster whose band 0 uses 7 and band 1 uses 8 as no-data has each swapped
+            to its own new value, non-matching cells preserved, and reports the new per-band
+            no-data.
+        """
+        band0 = np.array([[1.0, 2.0], [7.0, 4.0]], dtype="float32")
+        band1 = np.array([[5.0, 8.0], [5.0, 5.0]], dtype="float32")
+        ds = Dataset.create_from_array(
+            np.stack([band0, band1]),
+            top_left_corner=(0.0, 2.0),
+            cell_size=1.0,
+            epsg=4326,
+            no_data_value=[7.0, 8.0],
+        )
+        out = ds.change_no_data_value([-1.0, -2.0], old_value=[7.0, 8.0]).read_array()
+        assert out[0][1, 0] == -1.0, "band 0's old no-data (7) must become -1"
+        assert out[1][0, 1] == -2.0, "band 1's old no-data (8) must become -2"
+        assert out[0][0, 0] == 1.0 and out[1][0, 0] == 5.0, "non-matching cells preserved"
+
+    def test_dtype_error_fires_on_band_with_no_matching_cells(self):
+        """The dtype guard fires even when the old value matches zero cells in the band.
+
+        Test scenario:
+            The swap is attempted on every tile (via a setitem that raises), so a band with no
+            matching cells still raises NoDataValueError rather than silently skipping the guard.
+        """
+        arr = np.array([[1.0, 2.0], [3.0, 4.0]], dtype="float32")
+        ds = Dataset.create_from_array(
+            arr,
+            top_left_corner=(0.0, 2.0),
+            cell_size=1.0,
+            epsg=4326,
+            no_data_value=-9999.0,
+        )
+        original_read = ds.read_array
+
+        def mock_read(band=None, window=None, **kwargs):
+            """Wrap the real read but raise on assignment (mimic a bad dtype)."""
+            result = original_read(band=band, window=window, **kwargs)
+            mock_arr = MagicMock(wraps=result)
+            mock_arr.__setitem__ = lambda key, value: (_ for _ in ()).throw(
+                TypeError("incompatible type")
+            )
+            mock_arr.__getitem__ = result.__getitem__
+            return mock_arr
+
+        with patch.object(ds, "read_array", mock_read):
+            with pytest.raises(NoDataValueError):
+                ds.change_no_data_value(-1.0, old_value=-12345.0)
+
+    def test_preserves_raster_attribute_table(self):
+        """CreateCopy carries the raster attribute table (RAT) across the swap.
+
+        Test scenario:
+            An int band with a 2-row RAT keeps the RAT (row count and a value) after
+            change_no_data_value.
+        """
+        from osgeo import gdal
+
+        mem = gdal.GetDriverByName("MEM").Create("", 4, 4, 1, gdal.GDT_Int32)
+        mem.SetGeoTransform((0.0, 1.0, 0.0, 4.0, 0.0, -1.0))
+        band = mem.GetRasterBand(1)
+        band.WriteArray(np.array([[0, 1, 2, 3]] * 4, dtype="int32"))
+        band.SetNoDataValue(0)
+        rat = gdal.RasterAttributeTable()
+        rat.CreateColumn("class_name", gdal.GFT_String, gdal.GFU_Name)
+        rat.SetRowCount(2)
+        rat.SetValueAsString(0, 0, "water")
+        rat.SetValueAsString(1, 0, "land")
+        band.SetDefaultRAT(rat)
+        ds = Dataset(mem)
+
+        result = ds.change_no_data_value(255, old_value=0)
+        out_rat = result.raster.GetRasterBand(1).GetDefaultRAT()
+        assert out_rat is not None, "RAT must survive change_no_data_value"
+        assert out_rat.GetRowCount() == 2, "RAT row count must survive"
+        assert out_rat.GetValueAsString(0, 0) == "water", "RAT values must survive"
+
+    def test_inplace_with_path_updates_source_and_persists(self, tmp_path):
+        """inplace=True combined with path= updates the source and persists to disk.
+
+        Test scenario:
+            The source dataset is re-pointed at the disk-backed GeoTIFF, reports the new
+            no-data, reads the swapped cell, and the file reopens with the swap on disk.
+        """
+        arr = np.array([[1.0, 2.0], [-9999.0, 4.0]], dtype="float32")
+        ds = Dataset.create_from_array(
+            arr,
+            top_left_corner=(0.0, 2.0),
+            cell_size=1.0,
+            epsg=4326,
+            no_data_value=-9999.0,
+        )
+        out = tmp_path / "inplace.tif"
+        result = ds.change_no_data_value(-1.0, old_value=-9999.0, inplace=True, path=out)
+        assert result is ds, "inplace should return the source dataset"
+        assert out.exists(), "inplace + path must still write the file"
+        assert ds.no_data_value[0] == -1.0, "source must report the new no-data"
+        assert np.asarray(ds.read_array(band=0))[1, 0] == -1.0, "source cell swapped"
+        reopened = np.asarray(Dataset.read_file(str(out)).read_array(band=0))
+        assert reopened[1, 0] == -1.0, "disk file must hold the swap"
+
+    def test_partial_file_removed_when_swap_fails(self, tmp_path):
+        """A mid-stream failure on the disk path leaves no partial GeoTIFF behind.
+
+        Test scenario:
+            An assignment that raises during the tiled swap makes change_no_data_value raise
+            NoDataValueError and delete the partially-written file at path.
+        """
+        arr = np.array([[1.0, 2.0], [3.0, 4.0]], dtype="float32")
+        ds = Dataset.create_from_array(
+            arr,
+            top_left_corner=(0.0, 2.0),
+            cell_size=1.0,
+            epsg=4326,
+            no_data_value=-9999.0,
+        )
+        original_read = ds.read_array
+
+        def mock_read(band=None, window=None, **kwargs):
+            """Wrap the real read but raise on assignment to force a mid-stream failure."""
+            result = original_read(band=band, window=window, **kwargs)
+            mock_arr = MagicMock(wraps=result)
+            mock_arr.__setitem__ = lambda key, value: (_ for _ in ()).throw(
+                TypeError("incompatible type")
+            )
+            mock_arr.__getitem__ = result.__getitem__
+            return mock_arr
+
+        out = tmp_path / "partial.tif"
+        with patch.object(ds, "read_array", mock_read):
+            with pytest.raises(NoDataValueError):
+                ds.change_no_data_value(-1.0, old_value=-9999.0, path=out)
+        assert not out.exists(), "a failed disk-backed call must remove the partial file"
+
 
 class TestChangeNoDataAttrConversion:
     """Tests for _change_no_data_value_attr type conversion."""

From 88c1b26c3fef35079487d7275320974317f5c8e7 Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 22:54:21 +0200
Subject: [PATCH 11/22] refactor(vectorize): correct the auto-tile byte
 estimate and soften its memory claim

Sum the per-band itemsizes instead of assuming band 0's dtype, so a mixed-dtype stack is
estimated correctly for the tile=None auto-select. Reword the comments/docstring: the
tiled path avoids the full-band ndarray, but the resulting DataFrame still scales with the
surviving non-no-data cells, so 'never materialised whole' overstated the bound.

Refs #969
---
 src/pyramids/dataset/engines/vectorize.py | 30 +++++++++++++----------
 1 file changed, 17 insertions(+), 13 deletions(-)

diff --git a/src/pyramids/dataset/engines/vectorize.py b/src/pyramids/dataset/engines/vectorize.py
index b3a4491ba0..aa147468b5 100644
--- a/src/pyramids/dataset/engines/vectorize.py
+++ b/src/pyramids/dataset/engines/vectorize.py
@@ -53,9 +53,10 @@
 _CELL_ORDER = "_pyramids_cell_order"
 
 # When `to_feature_collection(tile=None)` auto-selects, read the whole array only if
-# it is at most this many bytes; a larger raster is read tile by tile so a huge or
-# /vsicurl source is never materialised whole (#969). 256 MiB keeps the fast
-# whole-array path for everyday rasters while bounding the pathological case.
+# it is at most this many bytes; a larger raster is read tile by tile so the full-band
+# ndarray is never allocated (#969) -- the resulting DataFrame still scales with the
+# surviving non-no-data cells. 256 MiB keeps the fast whole-array path for everyday
+# rasters while bounding the pathological case.
 _AUTO_TILE_BYTES = 256 * 1024 * 1024
 
 
@@ -239,11 +240,12 @@ def to_feature_collection(
                 Whether to read the raster tile by tile rather than in one pass, which
                 bounds the peak allocation on a large raster. `None` (default)
                 auto-selects: the whole array is read when it is at most ~256 MiB, else
-                the raster is tiled, so a huge or `/vsicurl` source is never
-                materialised whole. Pass `True`/`False` to force the choice. The rows
-                are the same cells in the same row-major order either way -- `mask`
-                included -- so it is a memory/throughput trade, not a change of result,
-                and `add_geometry` is safe with either.
+                the raster is tiled so the full-band ndarray is never allocated (the
+                resulting DataFrame still scales with the surviving non-no-data cells).
+                Pass `True`/`False` to force the choice. The rows are the same cells in
+                the same row-major order either way -- `mask` included -- so it is a
+                memory/throughput trade, not a change of result, and `add_geometry` is
+                safe with either.
             tile_size (int):
                 Tile size in cells, applied to both axes. Default is 256.
             touch (bool):
@@ -377,12 +379,14 @@ def to_feature_collection(
             src_ds = self._ds
 
         # None auto-selects on the array's byte size: keep the fast whole-array read
-        # for everyday rasters, tile a large one so it is never materialised whole
-        # (#969). The tiled path is byte-identical (same row-major rows), so this only
-        # trades memory for throughput. An explicit True/False overrides.
+        # for everyday rasters, tile a large one so the full-band ndarray is never
+        # allocated (#969). The tiled path is byte-identical (same row-major rows), so
+        # this only trades memory for throughput. An explicit True/False overrides.
+        # Sum the per-band itemsizes rather than assuming band 0's dtype, so a
+        # mixed-dtype stack is estimated correctly.
         if tile is None:
-            itemsize = np.dtype(src_ds.dtype[0]).itemsize
-            full_bytes = src_ds.rows * src_ds.columns * src_ds.band_count * itemsize
+            bytes_per_cell = sum(np.dtype(dt).itemsize for dt in src_ds.dtype)
+            full_bytes = src_ds.rows * src_ds.columns * bytes_per_cell
             tile = full_bytes > _AUTO_TILE_BYTES
 
         # Both branches must read `src_ds` -- the cropped dataset when a mask was

From 42d5b575ef310d1432fd60e16a47332202144803 Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 22:55:38 +0200
Subject: [PATCH 12/22] test(analysis): cover decimated footprint on a rotated,
 non-square-cell raster

The existing max_samples footprint tests used axis-aligned square-cell rasters, so the
geotransform rotation terms (gt[2]/gt[4]) were always 0 and a regression dropping their
scaling in _scaled_geotransform would pass. Add a fully-covered raster with non-zero
rotation and non-square cells, asserting exact and decimated footprints share total_bounds.

Refs #969
---
 .../analysis/test_dataset_vectorize.py        | 23 +++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/tests/dataset/analysis/test_dataset_vectorize.py b/tests/dataset/analysis/test_dataset_vectorize.py
index fc81a887db..da545bc4f6 100644
--- a/tests/dataset/analysis/test_dataset_vectorize.py
+++ b/tests/dataset/analysis/test_dataset_vectorize.py
@@ -462,6 +462,29 @@ def test_max_samples_decimates_and_keeps_extent(self):
             "Decimated footprint area should be roughly the block area"
         )
 
+    def test_max_samples_preserves_extent_on_rotated_non_square_raster(self):
+        """Decimated footprint keeps the exact extent on a rotated, non-square-cell raster.
+
+        Test scenario:
+            A fully-covered raster with non-zero rotation terms and non-square cells
+            footprints to the same total_bounds exact vs decimated, so a regression dropping
+            the geotransform rotation/scale scaling in `_scaled_geotransform` would be caught.
+        """
+        arr = np.ones((100, 100), dtype="float32")
+        geotransform = (0.0, 1.0, 0.2, 50.0, 0.15, -1.0)
+        ds = Dataset.create_from_array(
+            arr, geo=geotransform, epsg=3857, no_data_value=-9999.0
+        )
+        exact = ds.footprint()
+        approx = ds.footprint(max_samples=400)
+        assert exact is not None and approx is not None, "Both footprints must exist"
+        np.testing.assert_allclose(
+            approx.total_bounds,
+            exact.total_bounds,
+            atol=1e-6,
+            err_msg="decimated footprint must preserve the exact extent (rotated/non-square)",
+        )
+
 
 class TestToFeatureCollectionMaskTiling:
     """The mask must be honoured on both the tiled and the non-tiled path."""

From a1db517788d30eaa3a72e2c19b41bfc8fdfd8b17 Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 22:56:51 +0200
Subject: [PATCH 13/22] refactor(spatial): coerce the crop no-data value once,
 not per tile

_crop_aligned_tiled called _apply_mask_nodata per tile, which re-ran the per-band
_check_no_data_value dtype coercion on every window. Precompute the coerced list once
before the loop and pass it in; _apply_mask_nodata still validates it itself when called
without one (the non-tiled path), so behaviour is unchanged.

Refs #969
---
 src/pyramids/dataset/engines/spatial.py | 24 ++++++++++++++++++++----
 1 file changed, 20 insertions(+), 4 deletions(-)

diff --git a/src/pyramids/dataset/engines/spatial.py b/src/pyramids/dataset/engines/spatial.py
index b1f2964379..2cc238d8c4 100644
--- a/src/pyramids/dataset/engines/spatial.py
+++ b/src/pyramids/dataset/engines/spatial.py
@@ -1277,13 +1277,23 @@ def _assert_crop_aligned(
                 )
 
     def _apply_mask_nodata(
-        self, src_array: np.ndarray, mask_no_data: np.ndarray, band_count: int
+        self,
+        src_array: np.ndarray,
+        mask_no_data: np.ndarray,
+        band_count: int,
+        no_data_value: list | None = None,
     ) -> None:
-        """Write the source no-data value into the masked cells (per band)."""
+        """Write the source no-data value into the masked cells (per band).
+
+        `no_data_value` may be a caller-precomputed, dtype-checked per-band list; the
+        tiled crop passes it so the coercion runs once instead of per tile. When
+        `None` the multi-band path validates it here as before.
+        """
         if band_count > 1:
             # check the no_data_value complies with the src dtype before writing it
             # into cells (a band full of values may never use its no_data_value).
-            no_data_value = self._ds._check_no_data_value(self._ds.no_data_value)
+            if no_data_value is None:
+                no_data_value = self._ds._check_no_data_value(self._ds.no_data_value)
             for band in range(self._ds.band_count):
                 src_array[band, mask_no_data] = no_data_value[band]
         else:
@@ -1409,6 +1419,12 @@ def _crop_aligned_tiled(
             dst_obj: The destination Dataset the masked blocks are written into.
             band_count: Number of bands in the source raster.
         """
+        # Coerce the per-band no-data value once here rather than on every tile.
+        no_data_value = (
+            self._ds._check_no_data_value(self._ds.no_data_value)
+            if band_count > 1
+            else None
+        )
         for xoff, yoff, xsize, ysize in self._ds.io._tile_offsets():
             window = [xoff, yoff, xsize, ysize]
             # read_array() is called with no chunks=, so it always returns a
@@ -1416,7 +1432,7 @@ def _crop_aligned_tiled(
             mask_tile = cast(np.typing.NDArray, mask.read_array(band=0, window=window))
             src_tile = cast(np.typing.NDArray, self._ds.read_array(window=window))
             mask_no_data = is_no_data(mask_tile, mask_noval)
-            self._apply_mask_nodata(src_tile, mask_no_data, band_count)
+            self._apply_mask_nodata(src_tile, mask_no_data, band_count, no_data_value)
             if band_count > 1:
                 for band in range(band_count):
                     dst_obj.raster.GetRasterBand(band + 1).WriteArray(

From dd1c99fbeb418ec7bf6b5237f5435c2d2ca7354e Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 23:14:52 +0200
Subject: [PATCH 14/22] fix(analysis): short-circuit an empty domain in apply's
 func helper

When apply(elementwise=True) streamed a fully-no-data tile, domain_values was size 0; a
func that raises on array input (a legitimate per-pixel map like lambda v: 1 if v > 5
else 0) fell back to np.vectorize(func), which raises "cannot call 'vectorize' on size 0
inputs unless otypes is set". The whole-array path never hit this (the band domain is
non-empty), so the tiled path diverged. Return early on an empty domain -- out_array is
already the no-data fill -- restoring byte-identity and also hardening the whole-array
path for an all-no-data band.

Refs #969
---
 src/pyramids/dataset/engines/analysis.py      |  7 +++++
 .../analysis/test_apply_elementwise.py        | 27 +++++++++++++++++++
 2 files changed, 34 insertions(+)

diff --git a/src/pyramids/dataset/engines/analysis.py b/src/pyramids/dataset/engines/analysis.py
index d633fcf748..8e64ff55cf 100644
--- a/src/pyramids/dataset/engines/analysis.py
+++ b/src/pyramids/dataset/engines/analysis.py
@@ -421,6 +421,13 @@ def _apply_func_to_domain(func, src_array, out_array, no_data_value) -> None:
         """
         domain_mask = inside_domain(src_array, no_data_value)
         domain_values = src_array[domain_mask]
+        # An empty domain (an all-no-data tile, common when streaming) needs no
+        # write -- out_array is already the no-data fill -- and short-circuiting
+        # here avoids `np.vectorize(func)` raising "cannot call 'vectorize' on
+        # size 0 inputs" on a fully-masked tile, keeping the tiled path
+        # byte-identical to the whole-array pass (#969).
+        if domain_values.size == 0:
+            return
         try:
             out_array[domain_mask] = func(domain_values)
         except (ValueError, TypeError):
diff --git a/tests/dataset/analysis/test_apply_elementwise.py b/tests/dataset/analysis/test_apply_elementwise.py
index 8ca41233a5..595166ef42 100644
--- a/tests/dataset/analysis/test_apply_elementwise.py
+++ b/tests/dataset/analysis/test_apply_elementwise.py
@@ -106,3 +106,30 @@ def test_inplace_streaming(self):
         np.testing.assert_array_equal(
             ds.read_array(), arr * 2, err_msg="In-place streamed result must be doubled"
         )
+
+    def test_all_nodata_tile_with_vectorize_fallback_func(self):
+        """A fully-no-data tile does not crash the np.vectorize fallback path.
+
+        Test scenario:
+            A raster whose last 256-px tile is entirely no-data, transformed with a scalar
+            conditional func (which forces the np.vectorize fallback on array input), streams
+            without raising and matches the whole-array pass.
+        """
+        arr = (np.random.default_rng(3).random((300, 300)) * 10).astype("float64")
+        arr[256:300, 256:300] = -9999.0
+        ds = Dataset.create_from_array(
+            arr,
+            top_left_corner=(0.0, 15.0),
+            cell_size=0.05,
+            epsg=4326,
+            no_data_value=-9999.0,
+        )
+        func = lambda v: 1.0 if v > 5 else 0.0  # noqa: E731 - forces vectorize fallback
+        streamed = ds.apply(func, elementwise=True).read_array()
+        whole = ds.apply(func, elementwise=False).read_array()
+        np.testing.assert_array_equal(
+            streamed,
+            whole,
+            err_msg="empty-domain tile must not crash and must match the whole-array pass",
+        )
+        assert streamed[270, 270] == -9999.0, "no-data cell in the empty tile preserved"

From fc78362bd6c40161ec3a0fa14a339e488a36cf95 Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 23:15:47 +0200
Subject: [PATCH 15/22] fix(analysis): validate max_samples in
 footprint/plot_histogram decimation

max_samples is a public keyword on footprint and plot_histogram; max_samples=0 reached the
decimation factor and raised ZeroDivisionError, and a negative value produced a complex
factor and a 'complex has no __round__' TypeError. Reject max_samples < 1 at the read
boundary with a clear ValueError.

Refs #969
---
 src/pyramids/dataset/engines/analysis.py       |  7 +++++++
 .../dataset/analysis/test_dataset_vectorize.py | 18 ++++++++++++++++++
 2 files changed, 25 insertions(+)

diff --git a/src/pyramids/dataset/engines/analysis.py b/src/pyramids/dataset/engines/analysis.py
index 8e64ff55cf..9ac84df15c 100644
--- a/src/pyramids/dataset/engines/analysis.py
+++ b/src/pyramids/dataset/engines/analysis.py
@@ -1800,7 +1800,14 @@ def _read_decimated(self, band: int, max_samples: int | None) -> np.ndarray:
 
         Returns:
             np.ndarray: The band array, full-resolution or decimated.
+
+        Raises:
+            ValueError: `max_samples` is not `None` and is less than 1.
         """
+        if max_samples is not None and max_samples < 1:
+            raise ValueError(
+                f"max_samples must be a positive integer or None, got {max_samples}."
+            )
         rows = self._ds.rows
         cols = self._ds.columns
         total = rows * cols
diff --git a/tests/dataset/analysis/test_dataset_vectorize.py b/tests/dataset/analysis/test_dataset_vectorize.py
index da545bc4f6..262b351bf5 100644
--- a/tests/dataset/analysis/test_dataset_vectorize.py
+++ b/tests/dataset/analysis/test_dataset_vectorize.py
@@ -462,6 +462,24 @@ def test_max_samples_decimates_and_keeps_extent(self):
             "Decimated footprint area should be roughly the block area"
         )
 
+    @pytest.mark.parametrize("bad", [0, -5])
+    def test_max_samples_below_one_raises(self, bad):
+        """A non-positive max_samples is rejected with a clear ValueError.
+
+        Args:
+            bad: An invalid max_samples value (zero or negative).
+
+        Test scenario:
+            footprint(max_samples=0) and (max_samples=-5) raise ValueError instead of a
+            cryptic ZeroDivisionError / complex-round TypeError.
+        """
+        arr = np.ones((8, 8), dtype="float32")
+        ds = Dataset.create_from_array(
+            arr, top_left_corner=(0.0, 8.0), cell_size=1.0, epsg=3857, no_data_value=-9.0
+        )
+        with pytest.raises(ValueError, match="max_samples must be a positive integer"):
+            ds.footprint(max_samples=bad)
+
     def test_max_samples_preserves_extent_on_rotated_non_square_raster(self):
         """Decimated footprint keeps the exact extent on a rotated, non-square-cell raster.
 

From 058337147066ceb65512ab7481ef7939b6eb227d Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 23:16:25 +0200
Subject: [PATCH 16/22] fix(bands): make change_no_data_value(path=) cleanup
 best-effort so it never masks the real error

The disk-path failure cleanup unlinked the partial file with a bare unlink; if a
GDAL/OS build still held the GeoTIFF locked, that unlink would raise PermissionError and
replace the original NoDataValueError. Wrap each unlink in try/except OSError so a residual
lock only lets the file linger and the original exception always propagates.

Refs #969
---
 src/pyramids/dataset/engines/bands.py | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/pyramids/dataset/engines/bands.py b/src/pyramids/dataset/engines/bands.py
index 897c76daf4..7028241918 100644
--- a/src/pyramids/dataset/engines/bands.py
+++ b/src/pyramids/dataset/engines/bands.py
@@ -1469,12 +1469,17 @@ def change_no_data_value(
             # leave a half-written GeoTIFF behind on the disk path: release every
             # handle to the file (both the wrapper and the local `dst` reference,
             # or GDAL keeps the file locked on Windows) and delete the partial file
-            # and its sidecar before re-raising.
+            # and its sidecar before re-raising. The unlinks are best-effort: if a
+            # GDAL/OS build still holds the file, swallow the OSError so the cleanup
+            # never masks the original exception (the file just lingers).
             if path is not None:
                 new_dataset.close()
                 dst = None
-                Path(target).unlink(missing_ok=True)
-                Path(f"{target}.aux.xml").unlink(missing_ok=True)
+                for leftover in (target, f"{target}.aux.xml"):
+                    try:
+                        Path(leftover).unlink()
+                    except OSError:
+                        pass
             raise
         # Flush the block cache so a disk-backed GeoTIFF has the swapped pixels on
         # disk before it is reopened (a no-op for the in-memory driver).

From b9f3a50021e4f854818c4ae77587ae28a68f25b1 Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 23:17:57 +0200
Subject: [PATCH 17/22] test(analysis): tighten decimation assertions and note
 sieve's metadata copy

Tighten the loose decimation checks: the decimated footprint now must track the covered
block's extent (not merely stay inside the raster) and its area within rel=0.3; the
histogram test asserts the decimated samples are real band values (nearest-neighbour), not
just a small count. Document that sieve's gdal.Translate seed also carries the band color
table / RAT / scale-offset onto the result.

Refs #969
---
 src/pyramids/dataset/engines/analysis.py         |  3 +++
 tests/dataset/analysis/test_dataset_vectorize.py | 13 ++++++++-----
 tests/dataset/plot/test_plot_dataset.py          | 10 +++++++++-
 3 files changed, 20 insertions(+), 6 deletions(-)

diff --git a/src/pyramids/dataset/engines/analysis.py b/src/pyramids/dataset/engines/analysis.py
index 9ac84df15c..07ee6817e2 100644
--- a/src/pyramids/dataset/engines/analysis.py
+++ b/src/pyramids/dataset/engines/analysis.py
@@ -1112,6 +1112,9 @@ def sieve(
         # (geotransform, CRS, dtype, and no-data carried across in the C layer)
         # instead of a full-band ``ReadAsArray`` -> ``WriteArray`` NumPy round
         # trip, so the whole band is never materialised as a NumPy array (#969).
+        # gdal.Translate also carries the band's color table / RAT / scale-offset
+        # onto the result (the old bare-MEM seed dropped them); the sieved pixels
+        # are unchanged either way, so this only preserves more metadata.
         out_ds = gdal.Translate("", self._ds.raster, format="MEM", bandList=[band + 1])
         dst_band = out_ds.GetRasterBand(1)
 
diff --git a/tests/dataset/analysis/test_dataset_vectorize.py b/tests/dataset/analysis/test_dataset_vectorize.py
index 262b351bf5..1258ff0abb 100644
--- a/tests/dataset/analysis/test_dataset_vectorize.py
+++ b/tests/dataset/analysis/test_dataset_vectorize.py
@@ -452,13 +452,16 @@ def test_max_samples_decimates_and_keeps_extent(self):
         approx = ds.footprint(exclude_values=[0], max_samples=400)
         assert approx is not None, "Decimated footprint must still be produced"
         minx, miny, maxx, maxy = approx.total_bounds
-        assert minx >= -1e-6 and miny >= -1e-6, (
-            "Footprint must stay within raster bounds"
+        # The covered block spans x 60..160, y 50..150; decimation (~10-unit cells)
+        # blurs the edges, so allow a tile-width tolerance but require the footprint
+        # to track the block, not merely stay inside the raster.
+        assert minx == pytest.approx(60, abs=12) and maxx == pytest.approx(160, abs=12), (
+            f"decimated footprint x-extent must track the block, got {minx}..{maxx}"
         )
-        assert maxx <= 200 + 1e-6 and maxy <= 200 + 1e-6, (
-            "Footprint must stay within bounds"
+        assert miny == pytest.approx(50, abs=12) and maxy == pytest.approx(150, abs=12), (
+            f"decimated footprint y-extent must track the block, got {miny}..{maxy}"
         )
-        assert float(approx.geometry.area.sum()) == pytest.approx(10000, rel=0.5), (
+        assert float(approx.geometry.area.sum()) == pytest.approx(10000, rel=0.3), (
             "Decimated footprint area should be roughly the block area"
         )
 
diff --git a/tests/dataset/plot/test_plot_dataset.py b/tests/dataset/plot/test_plot_dataset.py
index 088635a762..d7ff6f8f70 100644
--- a/tests/dataset/plot/test_plot_dataset.py
+++ b/tests/dataset/plot/test_plot_dataset.py
@@ -159,7 +159,9 @@ def filter_kwargs(kw):
                 return {}
 
             def __init__(self, values, ax=None, **kwargs):
-                captured["n"] = np.asarray(values).size
+                vals = np.asarray(values)
+                captured["n"] = vals.size
+                captured["vals"] = vals
 
             def histogram(self, bins=15):
                 return ("fig", "ax", {})
@@ -169,12 +171,18 @@ def histogram(self, bins=15):
         ):
             dataset.plot_histogram(band=0, max_samples=400)
             decimated = captured["n"]
+            decimated_vals = captured["vals"]
             dataset.plot_histogram(band=0)
             full = captured["n"]
         assert full == 10000, f"Exact read must see every pixel, got {full}"
         assert decimated <= 500, (
             f"max_samples=400 must decimate to ~400 samples, got {decimated}"
         )
+        # nearest-neighbour decimation must hand the glyph real band values, not
+        # interpolated or out-of-range garbage.
+        assert np.isin(decimated_vals, np.arange(10000)).all(), (
+            "every decimated sample must be an actual band value (nearest-neighbour)"
+        )
 
     @pytest.mark.plot
     def test_invalid_color_scale_raises(self, src: Dataset):

From a352aeb631491bb3603b817fd4ea58a33d314a07 Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 23:22:40 +0200
Subject: [PATCH 18/22] test(analysis): split composite assertions flagged by
 SonarCloud (S9073)

Split the `assert X and Y` checks added in this PR into separate asserts so a failure
points at the exact condition: the decimated-footprint extent bounds, the both-footprints
existence check, the multi-band non-matching-cell preservation, and the variable-view
apply both-paths check.

Refs #969
---
 tests/dataset/analysis/test_dataset_vectorize.py | 13 ++++++-------
 tests/dataset/unit/test_unit_nodata.py           |  3 ++-
 tests/netcdf/samples/test_inherited_arg_ops.py   |  5 ++---
 3 files changed, 10 insertions(+), 11 deletions(-)

diff --git a/tests/dataset/analysis/test_dataset_vectorize.py b/tests/dataset/analysis/test_dataset_vectorize.py
index 1258ff0abb..eac3ff43af 100644
--- a/tests/dataset/analysis/test_dataset_vectorize.py
+++ b/tests/dataset/analysis/test_dataset_vectorize.py
@@ -455,12 +455,10 @@ def test_max_samples_decimates_and_keeps_extent(self):
         # The covered block spans x 60..160, y 50..150; decimation (~10-unit cells)
         # blurs the edges, so allow a tile-width tolerance but require the footprint
         # to track the block, not merely stay inside the raster.
-        assert minx == pytest.approx(60, abs=12) and maxx == pytest.approx(160, abs=12), (
-            f"decimated footprint x-extent must track the block, got {minx}..{maxx}"
-        )
-        assert miny == pytest.approx(50, abs=12) and maxy == pytest.approx(150, abs=12), (
-            f"decimated footprint y-extent must track the block, got {miny}..{maxy}"
-        )
+        assert minx == pytest.approx(60, abs=12), f"x-min must track the block, got {minx}"
+        assert maxx == pytest.approx(160, abs=12), f"x-max must track the block, got {maxx}"
+        assert miny == pytest.approx(50, abs=12), f"y-min must track the block, got {miny}"
+        assert maxy == pytest.approx(150, abs=12), f"y-max must track the block, got {maxy}"
         assert float(approx.geometry.area.sum()) == pytest.approx(10000, rel=0.3), (
             "Decimated footprint area should be roughly the block area"
         )
@@ -498,7 +496,8 @@ def test_max_samples_preserves_extent_on_rotated_non_square_raster(self):
         )
         exact = ds.footprint()
         approx = ds.footprint(max_samples=400)
-        assert exact is not None and approx is not None, "Both footprints must exist"
+        assert exact is not None, "exact footprint must exist"
+        assert approx is not None, "decimated footprint must exist"
         np.testing.assert_allclose(
             approx.total_bounds,
             exact.total_bounds,
diff --git a/tests/dataset/unit/test_unit_nodata.py b/tests/dataset/unit/test_unit_nodata.py
index 6d2d8014c7..c7f0cd99e1 100644
--- a/tests/dataset/unit/test_unit_nodata.py
+++ b/tests/dataset/unit/test_unit_nodata.py
@@ -746,7 +746,8 @@ def test_multi_band_per_band_old_value_list(self):
         out = ds.change_no_data_value([-1.0, -2.0], old_value=[7.0, 8.0]).read_array()
         assert out[0][1, 0] == -1.0, "band 0's old no-data (7) must become -1"
         assert out[1][0, 1] == -2.0, "band 1's old no-data (8) must become -2"
-        assert out[0][0, 0] == 1.0 and out[1][0, 0] == 5.0, "non-matching cells preserved"
+        assert out[0][0, 0] == 1.0, "band 0 non-matching cell must be preserved"
+        assert out[1][0, 0] == 5.0, "band 1 non-matching cell must be preserved"
 
     def test_dtype_error_fires_on_band_with_no_matching_cells(self):
         """The dtype guard fires even when the old value matches zero cells in the band.
diff --git a/tests/netcdf/samples/test_inherited_arg_ops.py b/tests/netcdf/samples/test_inherited_arg_ops.py
index d3a6136313..bb6b54c437 100644
--- a/tests/netcdf/samples/test_inherited_arg_ops.py
+++ b/tests/netcdf/samples/test_inherited_arg_ops.py
@@ -169,9 +169,8 @@ def test_apply_on_variable_view(tos):
     """
     default = tos.apply(lambda a: a + 1)
     streamed = tos.apply(lambda a: a + 1, elementwise=True)
-    assert default is not None and streamed is not None, (
-        "apply must return a dataset on a variable view for both paths"
-    )
+    assert default is not None, "default apply must return a dataset on a variable view"
+    assert streamed is not None, "elementwise apply must return a dataset on a variable view"
     np.testing.assert_array_equal(
         np.asarray(default.read_array(band=0)),
         np.asarray(streamed.read_array(band=0)),

From e8e52d4d6c333ec10735fe3095ee87b131a35ab2 Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 23:28:30 +0200
Subject: [PATCH 19/22] style: wrap review-added test lines to ruff-format's
 88-col width

---
 tests/dataset/analysis/test_dataset_vectorize.py | 14 +++++++++-----
 tests/dataset/unit/test_unit_nodata.py           |  6 ++++--
 tests/netcdf/samples/test_inherited_arg_ops.py   |  4 ++--
 3 files changed, 15 insertions(+), 9 deletions(-)

diff --git a/tests/dataset/analysis/test_dataset_vectorize.py b/tests/dataset/analysis/test_dataset_vectorize.py
index eac3ff43af..68a925001d 100644
--- a/tests/dataset/analysis/test_dataset_vectorize.py
+++ b/tests/dataset/analysis/test_dataset_vectorize.py
@@ -455,10 +455,10 @@ def test_max_samples_decimates_and_keeps_extent(self):
         # The covered block spans x 60..160, y 50..150; decimation (~10-unit cells)
         # blurs the edges, so allow a tile-width tolerance but require the footprint
         # to track the block, not merely stay inside the raster.
-        assert minx == pytest.approx(60, abs=12), f"x-min must track the block, got {minx}"
-        assert maxx == pytest.approx(160, abs=12), f"x-max must track the block, got {maxx}"
-        assert miny == pytest.approx(50, abs=12), f"y-min must track the block, got {miny}"
-        assert maxy == pytest.approx(150, abs=12), f"y-max must track the block, got {maxy}"
+        assert minx == pytest.approx(60, abs=12), f"x-min off block: {minx}"
+        assert maxx == pytest.approx(160, abs=12), f"x-max off block: {maxx}"
+        assert miny == pytest.approx(50, abs=12), f"y-min off block: {miny}"
+        assert maxy == pytest.approx(150, abs=12), f"y-max off block: {maxy}"
         assert float(approx.geometry.area.sum()) == pytest.approx(10000, rel=0.3), (
             "Decimated footprint area should be roughly the block area"
         )
@@ -476,7 +476,11 @@ def test_max_samples_below_one_raises(self, bad):
         """
         arr = np.ones((8, 8), dtype="float32")
         ds = Dataset.create_from_array(
-            arr, top_left_corner=(0.0, 8.0), cell_size=1.0, epsg=3857, no_data_value=-9.0
+            arr,
+            top_left_corner=(0.0, 8.0),
+            cell_size=1.0,
+            epsg=3857,
+            no_data_value=-9.0,
         )
         with pytest.raises(ValueError, match="max_samples must be a positive integer"):
             ds.footprint(max_samples=bad)
diff --git a/tests/dataset/unit/test_unit_nodata.py b/tests/dataset/unit/test_unit_nodata.py
index c7f0cd99e1..a7736d9948 100644
--- a/tests/dataset/unit/test_unit_nodata.py
+++ b/tests/dataset/unit/test_unit_nodata.py
@@ -824,7 +824,9 @@ def test_inplace_with_path_updates_source_and_persists(self, tmp_path):
             no_data_value=-9999.0,
         )
         out = tmp_path / "inplace.tif"
-        result = ds.change_no_data_value(-1.0, old_value=-9999.0, inplace=True, path=out)
+        result = ds.change_no_data_value(
+            -1.0, old_value=-9999.0, inplace=True, path=out
+        )
         assert result is ds, "inplace should return the source dataset"
         assert out.exists(), "inplace + path must still write the file"
         assert ds.no_data_value[0] == -1.0, "source must report the new no-data"
@@ -863,7 +865,7 @@ def mock_read(band=None, window=None, **kwargs):
         with patch.object(ds, "read_array", mock_read):
             with pytest.raises(NoDataValueError):
                 ds.change_no_data_value(-1.0, old_value=-9999.0, path=out)
-        assert not out.exists(), "a failed disk-backed call must remove the partial file"
+        assert not out.exists(), "failed disk call must remove the partial file"
 
 
 class TestChangeNoDataAttrConversion:
diff --git a/tests/netcdf/samples/test_inherited_arg_ops.py b/tests/netcdf/samples/test_inherited_arg_ops.py
index bb6b54c437..aa71897b70 100644
--- a/tests/netcdf/samples/test_inherited_arg_ops.py
+++ b/tests/netcdf/samples/test_inherited_arg_ops.py
@@ -169,8 +169,8 @@ def test_apply_on_variable_view(tos):
     """
     default = tos.apply(lambda a: a + 1)
     streamed = tos.apply(lambda a: a + 1, elementwise=True)
-    assert default is not None, "default apply must return a dataset on a variable view"
-    assert streamed is not None, "elementwise apply must return a dataset on a variable view"
+    assert default is not None, "default apply must return a dataset on a view"
+    assert streamed is not None, "elementwise apply must return a dataset on a view"
     np.testing.assert_array_equal(
         np.asarray(default.read_array(band=0)),
         np.asarray(streamed.read_array(band=0)),

From 2aa5cb474fb3454ce5e8347a21d17a078b73e09f Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 23:41:53 +0200
Subject: [PATCH 20/22] refactor(bands): reduce change_no_data_value cognitive
 complexity (SonarCloud S3776)

The added path=/cleanup logic pushed change_no_data_value's cognitive complexity to 22.
Extract three helpers -- _normalize_no_data_arg (per-band arg validation),
_swap_all_bands (the per-band swap loop), and _discard_partial_output (the disk-path
failure cleanup) -- so the method reads as a linear clone -> set -> swap -> flush with a
single try/except. Behaviour is unchanged.

Refs #969
---
 src/pyramids/dataset/engines/bands.py | 107 ++++++++++++++++----------
 1 file changed, 68 insertions(+), 39 deletions(-)

diff --git a/src/pyramids/dataset/engines/bands.py b/src/pyramids/dataset/engines/bands.py
index 7028241918..bfb0c9e723 100644
--- a/src/pyramids/dataset/engines/bands.py
+++ b/src/pyramids/dataset/engines/bands.py
@@ -1349,6 +1349,64 @@ def _change_no_data_value_attr(self, band: int, no_data_value) -> None:
                 self._ds.raster.GetRasterBand(band + 1).SetNoDataValue(no_data_value)
         self._ds._no_data_value[band] = no_data_value
 
+    def _normalize_no_data_arg(self, value: Any, name: str) -> list:
+        """Normalize a scalar or per-band no-data value to a list of length band_count.
+
+        Args:
+            value: A scalar (broadcast to every band) or a per-band list.
+            name: The argument name, used in the error message.
+
+        Returns:
+            list: A per-band list of length `band_count`.
+
+        Raises:
+            NoDataValueError: `value` is a list whose length is not `band_count`.
+        """
+        if not isinstance(value, list):
+            return [value] * self._ds.band_count
+        if len(value) != self._ds.band_count:
+            raise NoDataValueError(
+                f"{name} must be a scalar or a list of length band_count "
+                f"({self._ds.band_count}); got a list of length {len(value)}."
+            )
+        return value
+
+    def _swap_all_bands(
+        self, new_dataset: Dataset, new_value: Any, old_value: list | None
+    ) -> None:
+        """Swap every band's old no-data cells to the new value, tile by tile.
+
+        Args:
+            new_dataset: The cloned destination written in place.
+            new_value: Per-band new no-data values (already dtype-coerced).
+            old_value: Per-band old no-data values, or `None` to match NaN.
+        """
+        for band in range(self._ds.band_count):
+            band_old_value = old_value[band] if old_value is not None else None
+            self._swap_no_data_tiled(
+                new_dataset, band, band_old_value, new_value[band]
+            )
+
+    @staticmethod
+    def _discard_partial_output(new_dataset: Dataset, target: str) -> None:
+        """Release the handle and delete a partially-written disk output (best-effort).
+
+        Closes the wrapper (needed with the caller dropping its own `dst` reference,
+        or GDAL keeps the file locked on Windows) and unlinks the partial file plus
+        its sidecar, swallowing `OSError` so a residual lock never masks the original
+        exception -- the file just lingers.
+
+        Args:
+            new_dataset: The destination wrapper to close.
+            target: The output file path whose partial file/sidecar to remove.
+        """
+        new_dataset.close()
+        for leftover in (target, f"{target}.aux.xml"):
+            try:
+                Path(leftover).unlink()
+            except OSError:
+                pass
+
     def change_no_data_value(
         self,
         new_value: Any,
@@ -1424,20 +1482,9 @@ def change_no_data_value(
 
               ```
         """
-        if not isinstance(new_value, list):
-            new_value = [new_value] * self._ds.band_count
-        if len(new_value) != self._ds.band_count:
-            raise NoDataValueError(
-                f"new_value must be a scalar or a list of length band_count "
-                f"({self._ds.band_count}); got a list of length {len(new_value)}."
-            )
-        if old_value is not None and not isinstance(old_value, list):
-            old_value = [old_value] * self._ds.band_count
-        if old_value is not None and len(old_value) != self._ds.band_count:
-            raise NoDataValueError(
-                f"old_value must be a scalar or a list of length band_count "
-                f"({self._ds.band_count}); got a list of length {len(old_value)}."
-            )
+        new_value = self._normalize_no_data_arg(new_value, "new_value")
+        if old_value is not None:
+            old_value = self._normalize_no_data_arg(old_value, "old_value")
         # Clone the full header + pixels with GDAL's block-based CreateCopy so
         # every band's colour table, description, scale/offset, RAT and metadata
         # survive exactly (no explicit per-property copy to drift out of sync).
@@ -1450,36 +1497,18 @@ def change_no_data_value(
         dst = gdal.GetDriverByName(driver).CreateCopy(target, self._ds.raster, 0)
         new_dataset = self._ds.__class__(dst, "write")
         try:
-            # the new_value could change inside the _set_no_data_value method before it is used to set the no_data_value
-            # attribute in the gdal object/pyramids object and to fill the band.
+            # _set_no_data_value may coerce new_value to each band's dtype; read it
+            # back from the object so the swap below uses the stored values.
             new_dataset._set_no_data_value(new_value)
-            # now we have to use the no_data_value value in the no_data_value attribute in the Dataset object as it is
-            # updated.
             new_value = new_dataset.no_data_value
-            for band in range(self._ds.band_count):
-                # old_value is normalized to a per-band list above (matching
-                # new_value); index it per-band here too instead of comparing
-                # against the whole list.
-                band_old_value = old_value[band] if old_value is not None else None
-                self._swap_no_data_tiled(
-                    new_dataset, band, band_old_value, new_value[band]
-                )
+            self._swap_all_bands(new_dataset, new_value, old_value)
         except Exception:
-            # A mid-stream failure (e.g. a dtype-mismatch NoDataValueError) must not
-            # leave a half-written GeoTIFF behind on the disk path: release every
-            # handle to the file (both the wrapper and the local `dst` reference,
-            # or GDAL keeps the file locked on Windows) and delete the partial file
-            # and its sidecar before re-raising. The unlinks are best-effort: if a
-            # GDAL/OS build still holds the file, swallow the OSError so the cleanup
-            # never masks the original exception (the file just lingers).
+            # A mid-stream failure must not leave a half-written GeoTIFF behind:
+            # drop the local handle and let the helper release + delete the partial
+            # file before re-raising (best-effort, so it never masks the error).
             if path is not None:
-                new_dataset.close()
                 dst = None
-                for leftover in (target, f"{target}.aux.xml"):
-                    try:
-                        Path(leftover).unlink()
-                    except OSError:
-                        pass
+                self._discard_partial_output(new_dataset, target)
             raise
         # Flush the block cache so a disk-backed GeoTIFF has the swapped pixels on
         # disk before it is reopened (a no-op for the in-memory driver).

From 776a774b4648dec04a2086375d321f6c1c8834a7 Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Wed, 12 Aug 2026 23:45:19 +0200
Subject: [PATCH 21/22] style: collapse a one-line call in _swap_all_bands to
 ruff-format width

---
 src/pyramids/dataset/engines/bands.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/src/pyramids/dataset/engines/bands.py b/src/pyramids/dataset/engines/bands.py
index bfb0c9e723..2efbb045a8 100644
--- a/src/pyramids/dataset/engines/bands.py
+++ b/src/pyramids/dataset/engines/bands.py
@@ -1383,9 +1383,7 @@ def _swap_all_bands(
         """
         for band in range(self._ds.band_count):
             band_old_value = old_value[band] if old_value is not None else None
-            self._swap_no_data_tiled(
-                new_dataset, band, band_old_value, new_value[band]
-            )
+            self._swap_no_data_tiled(new_dataset, band, band_old_value, new_value[band])
 
     @staticmethod
     def _discard_partial_output(new_dataset: Dataset, target: str) -> None:

From 58b2e953bd2ecfcb69c48b555ead2c06b4a7e9ac Mon Sep 17 00:00:00 2001
From: Mostafa Farrag <moah.farag@gmail.com>
Date: Thu, 13 Aug 2026 00:12:04 +0200
Subject: [PATCH 22/22] test(io): make stream_transform memory bound
 build-agnostic

test_peak_memory_is_bounded_by_the_tile asserted a hard-coded Python-heap budget
(dense_bytes // 4) that does not hold across GDAL builds: on the macOS cibuildwheel GDAL
the streamed pass traced ~3.9 MB (vs the 0.5 MB budget) even though it tiles correctly,
because tracemalloc only sees the Python heap, not GDAL's C buffers, and the wheel build's
read/write path allocates differently. Replace the absolute budget with a relative one:
measure a whole-array pass (same read -> transform -> write, materialised at once) in the
same process and assert the streamed peak stays below it. The baseline shares the streamed
path's machinery, so any per-build overhead cancels and only the full-array data streaming
avoids remains -- locally the streamed peak is ~2% of the whole-array peak.
---
 tests/dataset/io/test_stream_transform.py | 46 +++++++++++++++++------
 1 file changed, 34 insertions(+), 12 deletions(-)

diff --git a/tests/dataset/io/test_stream_transform.py b/tests/dataset/io/test_stream_transform.py
index 32d1af2a58..75c727764f 100644
--- a/tests/dataset/io/test_stream_transform.py
+++ b/tests/dataset/io/test_stream_transform.py
@@ -174,12 +174,16 @@ def test_writes_into_a_provided_out_dataset(self):
         )
 
     def test_peak_memory_is_bounded_by_the_tile(self, tmp_path):
-        """Streaming to disk peaks near one tile, far below the whole-array size.
+        """Streaming to disk peaks below a whole-array pass, proving the tiled read.
 
         Test scenario:
-            Transform a 1000x1000 int16 raster to a disk output with 128-pixel tiles;
-            the traced Python peak must be a fraction of the dense-array size, proving
-            the whole raster is never materialised at once.
+            Transform a 1000x1000 int16 raster to a disk output with 128-pixel tiles, and
+            compare the traced Python peak against a *whole-array* pass (the same read ->
+            transform -> write, but materialised at once) in the same process. The streamed
+            peak must stay below the whole-array peak. This is a build-agnostic check on
+            purpose: the absolute figures depend on the GDAL build (tracemalloc only sees the
+            Python heap, not GDAL's C buffers), but the same operation done all-at-once is
+            always an upper bound on the tiled version, whatever the build.
         """
         rows = cols = 1000
         src_path = tmp_path / "big.tif"
@@ -190,17 +194,35 @@ def test_peak_memory_is_bounded_by_the_tile(self, tmp_path):
             epsg=4326,
             path=str(src_path),
         ).close()
-        ds = Dataset.read_file(str(src_path))
-        dense_bytes = rows * cols * 2  # int16
+
+        def add_one(block):
+            return block + 1
+
+        # Whole-array baseline: the same read -> transform -> write, but the full
+        # source and result arrays are held at once (no tiling).
+        whole_ds = Dataset.read_file(str(src_path))
         tracemalloc.start()
-        ds.io.stream_transform(
-            lambda tile: tile + 1, tile_size=128, path=str(tmp_path / "big_out.tif")
+        source_arr = whole_ds.read_array()
+        result = add_one(source_arr)
+        whole_out = Dataset.empty_like(whole_ds, path=str(tmp_path / "whole_out.tif"))
+        whole_out.write_array(result)
+        whole_out.close()
+        _, whole_peak = tracemalloc.get_traced_memory()
+        tracemalloc.stop()
+        del source_arr, result
+
+        # Streamed: the same transform, tile by tile straight to disk.
+        streamed_ds = Dataset.read_file(str(src_path))
+        tracemalloc.start()
+        streamed_ds.io.stream_transform(
+            add_one, tile_size=128, path=str(tmp_path / "big_out.tif")
         )
-        _, peak = tracemalloc.get_traced_memory()
+        _, tiled_peak = tracemalloc.get_traced_memory()
         tracemalloc.stop()
-        assert peak < dense_bytes // 4, (
-            f"stream_transform peaked at {peak / 1e6:.1f} MB; a whole-array pass "
-            f"would need {dense_bytes / 1e6:.1f} MB — the read was not tiled"
+
+        assert tiled_peak < whole_peak, (
+            f"stream_transform peaked at {tiled_peak / 1e6:.1f} MB, not below the "
+            f"whole-array pass's {whole_peak / 1e6:.1f} MB — the read was not tiled"
         )