Fill sparse COG tiles with nodata instead of crashing

.claude/sweep-accuracy-state.csv

@@ -12,7 +12,7 @@ emerging_hotspots,2026-04-30,,MEDIUM,2;3,MEDIUM: threshold_90 uses int() (trunca
 fire,2026-04-30,,,,All ops per-pixel (no accumulation/stencil/projected distance). NaN handled via x!=x; CUDA bounds use strict <; rdnbr and ros divisions guarded; CPU/GPU/dask paths algorithmically identical. No accuracy issues found.
 flood,2026-04-30,,MEDIUM,2;5,"MEDIUM (not fixed): dask backend preserves float32 input dtype while numpy promotes to float64 in flood_depth and curve_number_runoff; DataArray inputs for curve_number, mannings_n bypass scalar > 0 (and CN <= 100) range validation, silently producing NaN/garbage."
 focal,2026-03-30T13:00:00Z,1092,,,
-geotiff,2026-05-05,1498,HIGH,1;5,HIGH fixed: CPU fp_predictor_decode wrong byte-lane layout for multi-sample predictor=3 (GPU was correct). MEDIUMs also fixed on same PR: eager and streaming writers emit LONG8 strip/tile offsets in BigTIFF output (were LONG); VRT read honors AREA_OR_POINT=Point; VRT nodata cast uses source dtype instead of float32. LOW fixed: duplicate LERC codec block removed from _compression.py. 2026-05-05 deep pass: HIGH fixed in #1498 -- predictor=2 (horizontal differencing) did byte-wise differencing instead of sample-wise; uint8 was correct but uint16/int16/uint32/int32 silently corrupted pixel values both on read of GDAL/libtiff/tifffile/rasterio-written files and on write (their readers saw garbage). CPU + GPU + writer all reworked to dispatch per-dtype kernels at sample resolution with stride=samples. LOW deferred: ModelTiepointTag with multiple tiepoints uses only the first 6 doubles (no warning); rare in practice but silently produces wrong georeferencing for control-point-georeferenced files. No CRIT findings beyond the predictor=2 fix.
+geotiff,2026-05-06,1498,HIGH,1;5,HIGH fixed: CPU fp_predictor_decode wrong byte-lane layout for multi-sample predictor=3 (GPU was correct). MEDIUMs also fixed on same PR: eager and streaming writers emit LONG8 strip/tile offsets in BigTIFF output (were LONG); VRT read honors AREA_OR_POINT=Point; VRT nodata cast uses source dtype instead of float32. LOW fixed: duplicate LERC codec block removed from _compression.py. 2026-05-05 deep pass: HIGH fixed in #1498 -- predictor=2 (horizontal differencing) did byte-wise differencing instead of sample-wise; uint8 was correct but uint16/int16/uint32/int32 silently corrupted pixel values both on read of GDAL/libtiff/tifffile/rasterio-written files and on write (their readers saw garbage). CPU + GPU + writer all reworked to dispatch per-dtype kernels at sample resolution with stride=samples. LOW deferred: ModelTiepointTag with multiple tiepoints uses only the first 6 doubles (no warning); rare in practice but silently produces wrong georeferencing for control-point-georeferenced files. 2026-05-06 fourth pass: HIGH fixed in fix-sparse-cog-tile branch -- sparse COG tiles/strips (TileByteCounts/StripByteCounts == 0 from GDAL SPARSE_OK=TRUE) crashed local mmap reads ("incomplete or truncated stream"), returned uninitialised np.empty memory in the COG-HTTP path, and broke GPU reads. Fix: detect sparse blocks, pre-fill the result with the file's nodata sentinel (0 if unset), skip those entries in the decode loop. New tests cover tiled+stripped, with/without nodata, raw and accessor reads, plus GPU. Re-confirmed clean: PackBits decode against GDAL output, planar=2 + predictor=2 + uint16/uint32, non-tile-aligned dask streaming write, WhiteIsZero (PI=0) inversion. Deferred non-accuracy: tiled JPEG TIFFs from GDAL fail because reader does not inject JPEGTables (tag 347) into per-tile streams (interop crash, not numerical) and Orientation tag (274) is silently ignored. No CRIT findings.
 glcm,2026-05-01,1408,HIGH,2,"angle=None averaged NaN as 0, masking no-valid-pairs as zero texture; fixed via nanmean-style averaging"
 hillshade,2026-04-10T12:00:00Z,,,,"Horn's method correct. All backends consistent. NaN propagation correct. float32 adequate for [0,1] output."
 hydro,2026-04-30,,LOW,1,Only LOW: twi log(0)=-inf if fa=0 (out-of-contract); MFD weighted sum no Kahan (negligible). No CRIT/HIGH issues.

xrspatial/geotiff/__init__.py

@@ -1385,18 +1385,25 @@ def read_geotiff_gpu(source: str, *,
     finally:
         src.close()
 
-    # GPU decode: try GDS (SSD→GPU direct) first, then CPU mmap path
-    from ._gpu_decode import gpu_decode_tiles_from_file
-    arr_gpu = None
+    # GPU decode: try GDS (SSD→GPU direct) first, then CPU mmap path.
+    # Sparse tiles (byte_count == 0) are unsupported on the GPU pipeline;
+    # the CPU reader fills them with nodata and copies onto the GPU.
+    has_sparse_tile = any(bc == 0 for bc in byte_counts)
+    if has_sparse_tile:
+        arr_cpu, _ = read_to_array(source, overview_level=overview_level)
+        arr_gpu = cupy.asarray(arr_cpu)
+    else:
+        from ._gpu_decode import gpu_decode_tiles_from_file
+        arr_gpu = None
 
-    try:
-        arr_gpu = gpu_decode_tiles_from_file(
-            source, offsets, byte_counts,
-            tw, th, width, height,
-            compression, predictor, file_dtype, samples,
-        )
-    except Exception:
-        pass
+        try:
+            arr_gpu = gpu_decode_tiles_from_file(
+                source, offsets, byte_counts,
+                tw, th, width, height,
+                compression, predictor, file_dtype, samples,
+            )
+        except Exception:
+            pass
 
     if arr_gpu is None:
         # Fallback: extract tiles via CPU mmap, then GPU decode

xrspatial/geotiff/_reader.py

@@ -440,6 +440,38 @@ def _decode_strip_or_tile(data_slice, compression, width, height, samples,
 _NATIVE_ORDER = '<' if _sys.byteorder == 'little' else '>'
 
 
+def _sparse_fill_value(ifd: IFD, dtype: np.dtype):
+    """Resolve the fill value for sparse tiles/strips.
+
+    A sparse TIFF entry has TileByteCounts/StripByteCounts == 0 (and
+    typically the matching Offset == 0). GDAL emits these for SPARSE_OK
+    files where blocks containing only the nodata value are omitted.
+    The reader is expected to materialise such blocks as nodata, or
+    zero when nodata is unset (the default per the GDAL convention).
+    """
+    nodata_str = ifd.nodata_str
+    if nodata_str is not None:
+        try:
+            v = float(nodata_str)
+            if dtype.kind == 'f':
+                return dtype.type(v)
+            if not math.isnan(v) and not math.isinf(v):
+                return dtype.type(int(v))
+        except (TypeError, ValueError):
+            pass
+    return dtype.type(0)
+
+
+def _has_sparse(byte_counts) -> bool:
+    """Return True if any tile/strip is empty (byte_count == 0)."""
+    if byte_counts is None:
+        return False
+    for bc in byte_counts:
+        if bc == 0:
+            return True
+    return False
+
+
 # ---------------------------------------------------------------------------
 # Strip reader
 # ---------------------------------------------------------------------------
@@ -502,7 +534,17 @@ def _read_strips(data: bytes, ifd: IFD, header: TIFFHeader,
 
     _check_dimensions(out_w, out_h, samples, max_pixels)
 
-    if samples > 1:
+    # Sparse strips (StripByteCounts == 0) must materialise as nodata or 0
+    # rather than be decoded.  Pre-fill the result so any skipped strips
+    # land on a known fill value.
+    sparse = _has_sparse(byte_counts)
+    if sparse:
+        fill = _sparse_fill_value(ifd, dtype)
+        if samples > 1:
+            result = np.full((out_h, out_w, samples), fill, dtype=dtype)
+        else:
+            result = np.full((out_h, out_w), fill, dtype=dtype)
+    elif samples > 1:
         result = np.empty((out_h, out_w, samples), dtype=dtype)
     else:
         result = np.empty((out_h, out_w), dtype=dtype)
@@ -518,6 +560,9 @@ def _read_strips(data: bytes, ifd: IFD, header: TIFFHeader,
                 global_idx = band_offset + strip_idx
                 if global_idx >= len(offsets):
                     continue
+                if byte_counts[global_idx] == 0:
+                    # Sparse strip: result is already pre-filled.
+                    continue
                 strip_row = strip_idx * rps
                 strip_rows = min(rps, height - strip_row)
                 if strip_rows <= 0:
@@ -544,6 +589,9 @@ def _read_strips(data: bytes, ifd: IFD, header: TIFFHeader,
             strip_rows = min(rps, height - strip_row)
             if strip_rows <= 0:
                 continue
+            if byte_counts[strip_idx] == 0:
+                # Sparse strip: result is already pre-filled.
+                continue
 
             strip_data = data[offsets[strip_idx]:offsets[strip_idx] + byte_counts[strip_idx]]
             strip_pixels = _decode_strip_or_tile(
@@ -638,11 +686,24 @@ def _read_tiles(data: bytes, ifd: IFD, header: TIFFHeader,
     # would mask the problem, and the GPU path reads OOB. See issue #1219.
     validate_tile_layout(ifd)
 
-    _alloc = np.zeros if window is not None else np.empty
-    if samples > 1:
-        result = _alloc((out_h, out_w, samples), dtype=dtype)
+    # Sparse tiles (TileByteCounts == 0) must materialise as nodata or 0
+    # rather than be decoded.  Pre-fill the result so any skipped tiles
+    # land on a known fill value; otherwise sparse regions would leak
+    # uninitialised memory (full-image read) or stay zeroed regardless
+    # of the file's nodata setting (windowed read).
+    sparse = _has_sparse(byte_counts)
+    if sparse:
+        fill = _sparse_fill_value(ifd, dtype)
+        if samples > 1:
+            result = np.full((out_h, out_w, samples), fill, dtype=dtype)
+        else:
+            result = np.full((out_h, out_w), fill, dtype=dtype)
     else:
-        result = _alloc((out_h, out_w), dtype=dtype)
+        _alloc = np.zeros if window is not None else np.empty
+        if samples > 1:
+            result = _alloc((out_h, out_w, samples), dtype=dtype)
+        else:
+            result = _alloc((out_h, out_w), dtype=dtype)
 
     tile_row_start = r0 // th
     tile_row_end = min(math.ceil(r1 / th), tiles_down)
@@ -652,7 +713,8 @@ def _read_tiles(data: bytes, ifd: IFD, header: TIFFHeader,
     band_count = samples if (planar == 2 and samples > 1) else 1
     tiles_per_band = tiles_across * tiles_down
 
-    # Build list of tiles to decode
+    # Build list of tiles to decode.  Sparse tiles (byte_count==0) are
+    # skipped here -- the result is pre-filled with the sparse fill value.
     tile_jobs = []
     for band_idx in range(band_count):
         band_tile_offset = band_idx * tiles_per_band if band_count > 1 else 0
@@ -663,6 +725,8 @@ def _read_tiles(data: bytes, ifd: IFD, header: TIFFHeader,
                 tile_idx = band_tile_offset + tr * tiles_across + tc
                 if tile_idx >= len(offsets):
                     continue
+                if byte_counts[tile_idx] == 0:
+                    continue
                 tile_jobs.append((band_idx, tr, tc, tile_idx, tile_samples))
 
     # Decode tiles in parallel when the work per tile is large enough to
@@ -815,7 +879,17 @@ def _read_cog_http(url: str, overview_level: int | None = None,
     # TileOffsets length. See issue #1219.
     validate_tile_layout(ifd)
 
-    if samples > 1:
+    # Sparse tiles (TileByteCounts == 0) need to land on the file's nodata
+    # value (or 0 if unset) rather than uninitialised memory.  Detect them
+    # up front so the result buffer is pre-filled before tile placement.
+    sparse = _has_sparse(byte_counts)
+    if sparse:
+        fill = _sparse_fill_value(ifd, dtype)
+        if samples > 1:
+            result = np.full((height, width, samples), fill, dtype=dtype)
+        else:
+            result = np.full((height, width), fill, dtype=dtype)
+    elif samples > 1:
         result = np.empty((height, width, samples), dtype=dtype)
     else:
         result = np.empty((height, width), dtype=dtype)

xrspatial/geotiff/tests/test_sparse_cog.py

@@ -0,0 +1,114 @@
+"""Tests for sparse-block GeoTIFF reads.
+
+GDAL's ``SPARSE_OK=TRUE`` creation option records all-nodata tiles or
+strips with ``TileByteCounts == 0`` and ``TileOffsets == 0`` (same for
+strip layout).  The reader has to materialise such blocks as the file's
+nodata value (or zero when no nodata is set), not crash, and not return
+uninitialised memory.
+"""
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+rasterio = pytest.importorskip('rasterio')
+
+from xrspatial.geotiff import open_geotiff
+from xrspatial.geotiff._reader import read_to_array
+
+
+def _write_sparse_tiled(path, *, dtype='uint16', nodata=0,
+                        compress='DEFLATE', filled_value=100):
+    """Write a 128x128 raster where only the top-left 64x64 tile is filled."""
+    profile = {
+        'driver': 'GTiff', 'dtype': dtype,
+        'height': 128, 'width': 128, 'count': 1,
+        'tiled': True, 'blockxsize': 64, 'blockysize': 64,
+        'compress': compress, 'SPARSE_OK': 'TRUE',
+    }
+    if nodata is not None:
+        profile['nodata'] = nodata
+    fill = np.full((64, 64), filled_value, dtype=np.dtype(dtype))
+    with rasterio.open(path, 'w', **profile) as dst:
+        dst.write(fill, 1, window=rasterio.windows.Window(0, 0, 64, 64))
+
+
+def _write_sparse_stripped(path, *, dtype='uint16', nodata=0):
+    profile = {
+        'driver': 'GTiff', 'dtype': dtype,
+        'height': 128, 'width': 128, 'count': 1,
+        'tiled': False, 'blockysize': 16,
+        'compress': 'DEFLATE', 'SPARSE_OK': 'TRUE',
+    }
+    if nodata is not None:
+        profile['nodata'] = nodata
+    fill = np.full((32, 128), 200, dtype=np.dtype(dtype))
+    with rasterio.open(path, 'w', **profile) as dst:
+        dst.write(fill, 1, window=rasterio.windows.Window(0, 0, 128, 32))
+
+
+class TestSparseTiles:
+    def test_sparse_tile_with_nodata_round_trips(self, tmp_path):
+        path = str(tmp_path / 'sparse_nodata.tif')
+        _write_sparse_tiled(path, nodata=0)
+
+        arr = open_geotiff(path)
+        arr_np = np.asarray(arr)
+        # Filled tile carries the data, sparse tiles are NaN (nodata=0
+        # promoted by the accessor).
+        assert arr_np[:64, :64].sum() == 64 * 64 * 100
+        assert np.all(np.isnan(arr_np[:64, 64:]))
+        assert np.all(np.isnan(arr_np[64:, :]))
+
+    def test_sparse_tile_without_nodata_fills_zero(self, tmp_path):
+        path = str(tmp_path / 'sparse_no_nodata.tif')
+        _write_sparse_tiled(path, nodata=None)
+
+        arr = open_geotiff(path)
+        arr_np = np.asarray(arr)
+        assert arr_np.dtype == np.uint16
+        assert arr_np[:64, :64].sum() == 64 * 64 * 100
+        assert arr_np[:64, 64:].sum() == 0
+        assert arr_np[64:, :].sum() == 0
+
+    def test_sparse_tile_raw_read_uses_nodata_sentinel(self, tmp_path):
+        """``read_to_array`` returns raw values; sparse tiles == nodata."""
+        path = str(tmp_path / 'sparse_raw.tif')
+        _write_sparse_tiled(path, nodata=255)
+
+        arr, geo = read_to_array(path)
+        assert arr.dtype == np.uint16
+        assert arr[:64, :64].sum() == 64 * 64 * 100
+        assert np.all(arr[:64, 64:] == 255)
+        assert np.all(arr[64:, :] == 255)
+
+
+class TestSparseStrips:
+    def test_sparse_strip_with_nodata(self, tmp_path):
+        path = str(tmp_path / 'sparse_strips.tif')
+        _write_sparse_stripped(path, nodata=0)
+
+        arr = open_geotiff(path)
+        arr_np = np.asarray(arr)
+        assert arr_np[:32, :].sum() == 32 * 128 * 200
+        # Empty strips below row 32 land on nodata -> NaN via accessor
+        assert np.all(np.isnan(arr_np[32:, :]))
+
+
+@pytest.mark.skipif(
+    not pytest.importorskip('cupy', reason='cupy required'),
+    reason='cupy required',
+)
+class TestSparseTilesGPU:
+    def test_sparse_tile_gpu_round_trip(self, tmp_path):
+        path = str(tmp_path / 'sparse_gpu.tif')
+        _write_sparse_tiled(path, nodata=0)
+
+        arr = open_geotiff(path, gpu=True)
+        # GPU read returns raw values (no nodata->NaN promotion); sparse
+        # tiles become the nodata sentinel (0 here).
+        host = arr.data.get()
+        assert host.dtype == np.uint16
+        assert host[:64, :64].sum() == 64 * 64 * 100
+        assert host[:64, 64:].sum() == 0
+        assert host[64:, :].sum() == 0