diff --git a/.claude/sweep-accuracy-state.csv b/.claude/sweep-accuracy-state.csv index 939ccd0d..4bebbdc8 100644 --- a/.claude/sweep-accuracy-state.csv +++ b/.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. diff --git a/xrspatial/geotiff/__init__.py b/xrspatial/geotiff/__init__.py index 95b1ce15..dcb382e1 100644 --- a/xrspatial/geotiff/__init__.py +++ b/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 diff --git a/xrspatial/geotiff/_reader.py b/xrspatial/geotiff/_reader.py index 2770c250..aacf89a7 100644 --- a/xrspatial/geotiff/_reader.py +++ b/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) diff --git a/xrspatial/geotiff/tests/test_sparse_cog.py b/xrspatial/geotiff/tests/test_sparse_cog.py new file mode 100644 index 00000000..68fd7aba --- /dev/null +++ b/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