diff --git a/.claude/sweep-accuracy-state.csv b/.claude/sweep-accuracy-state.csv
index ea5b72b8f..e7a4074a1 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-12,1655,MEDIUM,2;5,"Pass 19 (2026-05-12): MEDIUM fixed -- issue #1655. read_vrt silently dropped 0 on a SimpleSource because of src.nodata or nodata at _vrt.py:370. Python treats 0.0 as falsy, so the per-source sentinel fell through to the band-level (or None when missing) and pixels equal to 0.0 in the source file survived as valid data. The in-code comment acknowledged the quirk as backward compat, but the resulting behaviour silently biased every NaN-aware aggregation on VRT mosaics whose sources used 0 as a sentinel (a common convention for unsigned remote-sensing imagery). Fix: src_nodata = src.nodata if src.nodata is not None else nodata. Five regression tests in test_vrt_source_nodata_zero_1655.py covering source NODATA=0, integer XML literal, non-zero unchanged, band-level NoDataValue=0 still honoured, and source-overrides-band precedence. All 100 vrt-related geotiff tests still pass; 3 pre-existing test_features.py matplotlib palette failures unrelated. Categories: Cat 2 (NaN propagation) + Cat 5 (backend inconsistency: read_geotiff masks 0 correctly when GDAL_NODATA tag is set; only VRT path was broken). | Pass 18 (2026-05-11): MEDIUM fixed -- issue #1642. PR #1641 (issue #1640) inherited level-0 georef on overview reads but kept the level-0 origin_x/origin_y unchanged. That is correct for PixelIsArea (origin = upper-left corner of pixel (0,0)) but wrong for PixelIsPoint (origin = center of pixel (0,0), GeoKey 1025 = 2). For a 1024x1024 PixelIsPoint COG with 10 m pixels and origin (0, 0), open_geotiff(overview_level=1) returned x[:3]=[0,20,40] instead of [5,25,45] (level-1 pixel 0 covers level-0 pixels 0-1 whose centers are 0 and 10, centroid 5); same for y. Downstream sel/interp/reproject silently snaps to the wrong pixel for any DEM-style PixelIsPoint COG (USGS, OpenTopography, Copernicus DEM). Categories: Cat 3 (off-by-one / boundary handling) + Cat 5 (raster_type-dependent backend convention). Fix: in extract_geo_info_with_overview_inheritance (_geotags.py), pick the effective raster_type first (overview-declared if non-default, otherwise inherited from parent), then when it is PixelIsPoint apply origin_shift = (scale - 1) * 0.5 * pixel_size_lvl0 along each axis before building the new GeoTransform. PixelIsArea path is byte-equivalent. 13 regression tests in test_overview_pixel_is_point_1642.py: centroid identity across all 4 backends, transform tuple across all 4 backends, uniform grid step, unit-level helper tests for both raster_types via stubbed extract_geo_info, own-geokeys-not-clobbered path on PixelIsPoint, and a PixelIsArea regression check. All 1397 existing non-network geotiff tests still pass (3 pre-existing matplotlib palette failures unrelated). Deferred LOW: non-power-of-two overview dimensions cause scale = base_w/ov_w to diverge from the true 2^level reduction (writer drops the right/bottom strip via h2=(h//2)*2; for h=1023 a level-1 overview has 511 rows so scale=2.0019 not 2.0). Fix would need to either (a) emit explicit geo tags on overview IFDs from the writer or (b) pass the level number into the inheritance helper; neither is a one-line change and the resulting coord error is sub-pixel of level 0. | Pass 17 (2026-05-11): MEDIUM fixed -- issue #1634. open_geotiff eager path windowed read produced confusing CoordinateValidationError when window extended past source extent. read_to_array clamped the window internally and returned a smaller array, but the eager code path used unclamped window indices for y/x coord generation (xrspatial/geotiff/__init__.py lines 562-572), so the coord array length differed from the data and xarray refused to construct the DataArray. Same bug affected the windowed transform shift in _populate_attrs_from_geo_info. The dask path (read_geotiff_dask) already validated up front since #1561, raising a clear ValueError with the format 'window=... is outside the source extent (HxW) or has non-positive size.' so the two backends diverged on the contract. Fix: validate the window up front in open_geotiff's eager branch via _read_geo_info (metadata-only read, no extra pixel cost) using the exact same condition the dask path uses, raising the same ValueError message format. Reproduction: 10x10 raster + window=(5,5,15,15) on eager raised CoordinateValidationError('conflicting sizes ... length 5 ... length 10'); now raises ValueError('window=(5, 5, 15, 15) is outside the source extent (10x10) or has non-positive size.'). Categories: Cat 3 (off-by-one / boundary handling) + Cat 5 (backend inconsistency). 12 regression tests in test_window_out_of_bounds_1634.py: negative start, past-right-edge, past-bottom-edge, past-both-edges, zero-size, inverted window, full-extent ok, interior subset, edge-aligned, eager-vs-dask parity, message-format parity, issue reproducer. All 1286 existing non-network geotiff tests still pass. | Pass 16 (2026-05-11): HIGH fixed -- issue #1623. to_geotiff(cog=True, overview_resampling='cubic', nodata=) on a float raster with NaN regions produced overview pixels with severe ringing artefacts near nodata borders. Same class of bug as #1613 but for the cubic branch: writer rewrites NaN to the sentinel upstream, then _block_reduce_2d(method=cubic) handed the sentinel-poisoned array straight to scipy.ndimage.zoom(order=3). The cubic spline blended the sentinel (e.g. -9999) into neighbouring cells, producing values like 1133.44, -10290.08 where the data was a constant 100. Repro on 16x16 float32 with a 4x4 NaN corner showed 18 polluted pixels in the 8x8 overview. Fix: when nodata is supplied on a float dtype and the sentinel is found, mask sentinel to NaN, run cubic with prefilter=False so a single NaN cannot poison the entire row/column (default B-spline prefilter is global), then rewrite any NaN in the result back to the sentinel. prefilter=False only fires when a sentinel is present so the non-nodata cubic semantics are unchanged. GPU side: _block_reduce_2d_gpu previously raised on method='cubic'; added a CPU fallback (same pattern as 'mode') so GPU writer produces byte-equivalent overviews. GPU_OVERVIEW_METHODS now includes 'cubic'. 12 regression tests in test_cog_cubic_overview_nodata_1623.py (helper no-ringing, poisoning repro, no-nodata unchanged, end-to-end round-trip, GPU fallback, CPU/GPU byte-match, +/-inf nodata mask, NaN-sentinel no-op, GPU_OVERVIEW_METHODS contract). All 1256 existing geotiff tests still pass (3 pre-existing matplotlib failures unrelated). | Pass 15 (2026-05-11): HIGH fixed -- issue #1613. to_geotiff(cog=True, nodata=) on a float raster with NaN produced a corrupted overview pyramid. The NaN-to-sentinel rewrite in __init__.py:1202 (CPU) and :2852 (GPU write_geotiff_gpu) ran BEFORE _make_overview / make_overview_gpu, so the nan-aware aggregations (np.nanmean/min/max/median, cupy.nanmean/min/max/median) saw the sentinel as a real number and biased every overview pixel. Reproduction with -9999 sentinel produced [[-4998.75,-4997.75],..] where np.nanmean gives [[1.5,3.5],..]. Both CPU and GPU paths affected; backend results matched each other but were both wrong (CAT 2 NaN propagation + CAT 5 documents the parity). Fix: _block_reduce_2d / _block_reduce_2d_gpu accept a nodata kwarg that masks the sentinel back to NaN for float dtypes before the reduction; the writer's overview loop passes nodata in, then rewrites all-sentinel reductions (which surface as NaN from the reducer) back to the sentinel for the on-disk pyramid. 11 regression tests in test_cog_overview_nodata_1613.py (CPU mean / partial-block / min/max/median / no-nodata passthrough / helper kwarg / all-sentinel block / GPU mean / GPU helper / CPU-GPU agreement). All 235 nodata/overview/cog tests still pass. | Pass 14 (2026-05-11): HIGH fixed -- issue #1611. read_vrt(band=None) on a multi-band integer VRT with per-band tags only masks band 0's sentinel. __init__.py lines 2795-2809 in read_vrt apply vrt.bands[0].nodata to the full ndim==3 array; bands 1+ keep their integer sentinels as literal finite values (e.g. 65000 surfaces as 65000.0 after the dtype=float64 cast, not NaN). Float-VRT path masks per-band correctly in _vrt._read_data lines 296-297 + 347-351. PR #1602 fixed the single-band band=N case for issue #1598; the band=None multi-band case is the same class of bug. Repro: 2-band uint16 VRT with NoDataValue 65535 / 65000 returns r.values[1,1,1] == 65000.0 instead of NaN; r.values[1,1,0] is NaN (band 0 sentinel masked). Fix scope: in read_vrt, when band is None, iterate over vrt.bands and mask each arr[..., i] slice against its own (gated by the same _int_nodata_in_range guard PR #1583 introduced). Severity HIGH (Cat 2 NaN propagation + Cat 5 backend inconsistency: identical input semantics produce different masking outcomes based on dtype, with finite garbage values where NaN expected). Fix in PR #1612: walks vrt.bands when band is None and ndim==3, masks each arr[..., i] slice against its own via the refactored _sentinel_for_dtype helper (reuses PR #1583's range guard so out-of-range/non-finite/fractional sentinels are a no-op). attrs['nodata'] still carries band 0's sentinel for band=None reads (documented contract). 7 regression tests in test_vrt_multiband_int_nodata_1611.py: uint16 per-band, int32 negative, mixed presence, dtype preservation when no sentinel hit, out-of-range gating, band=N non-regression, attrs contract. 135 existing vrt/nodata geotiff tests still pass. | Pass 13 (2026-05-11): HIGH fixed -- issue #1599. write_geotiff_gpu (and to_geotiff gpu=True) emitted raw NaN bytes for missing pixels even when nodata= was supplied, while the CPU writer substituted NaN with the sentinel before encoding. xrspatial-only round-trips were unaffected (the reader masks both NaN and the sentinel), but external readers (rasterio/GDAL/QGIS) that mask only on the GDAL_NODATA tag saw NaN pixels as valid data -- rasterio reported 100% valid pixels on a 25-NaN file vs CPU's 25-invalid report. Root cause: __init__.py lines 2579-2587 jumped from shape/dtype resolution straight to compression, missing the equivalent of the CPU writer's NaN-to-sentinel rewrite at to_geotiff line ~1156. Fix: cupy.isnan + masked write on a defensive copy of arr, gated on np_dtype.kind=='f' and not np.isnan(float(nodata)). Caller's CuPy buffer preserved (copy before mutate). 7 regression tests in test_gpu_writer_nan_sentinel_1599.py: substitution lands as sentinel, CPU/GPU byte-equivalent, caller buffer not mutated, no-NaN no-op, NaN sentinel skips substitution, rasterio sees identical invalid count on CPU/GPU, multiband 3D path. All other GPU writer tests still pass (50 passed across band-first, attrs, nodata, dask+cupy, writer, nodata aliases). | Pass 12 (2026-05-11): HIGH fixed -- issue #1581. Reading a uint TIFF with a negative GDAL_NODATA sentinel (e.g. uint16 + -9999) raised OverflowError on every backend because the nodata-mask code did arr.dtype.type(int(nodata)) with no range check. Three identical cast sites in __init__.py (numpy eager, _apply_nodata_mask_gpu, _delayed_read_window) plus _resolve_masked_fill and _sparse_fill_value in _reader.py. Fix: _int_nodata_in_range helper gates the cast; out-of-range sentinels are a no-op for value matching (the file can never contain that value), file dtype is preserved, attrs['nodata'] still surfaces the original sentinel so write round-trips keep the GDAL_NODATA tag intact. Matches rasterio behavior. 8 regression tests in test_nodata_out_of_range_1581.py cover the helper, both eager and dask read paths, in-range sentinel non-regression, and GPU helper (cupy-gated). | Pass 11 (2026-05-10): CLEAN. Audited the one additional commit since pass 10 -- #1559 (PR 1548, Centralise GeoTIFF attrs population across all read backends). Refactor extracts _populate_attrs_from_geo_info helper and routes eager numpy, dask, GPU stripped, GPU tiled read paths through it; before the fix dask only emitted crs/transform/raster_type/nodata while numpy emitted the full attrs set including x/y_resolution, resolution_unit, image_description, extra_samples, GDAL metadata, and the CRS-description fields. No data-path arithmetic touched; only attrs dict population. Windowed origin math (origin_x + c0*pixel_width, origin_y + r0*pixel_height) verified to produce -98.0 / 48.75 origin for window=(10,20,50,70) on a (0.1,-0.125) pixel-size raster, with PixelIsArea half-pixel offset preserved on coord lookups (-97.95, 48.6875). Cross-backend attrs parity re-verified: numpy/dask/cupy all emit identical key set on deflate+predictor3+nodata round-trip (crs, crs_wkt, nodata, transform, x_resolution, y_resolution). Data bit-parity re-verified across numpy/dask/cupy on same payload (np.array_equal with equal_nan=True). test_attrs_parity_1548.py (5 tests), test_reader.py/test_writer.py/test_dask_cupy_combined.py (25 tests), GPU orientation/predictor2-BE/LERC-mask/nodata/byteswap suites (65 tests) all green. No accuracy or backend-divergence findings. | Pass 10 (2026-05-10): CLEAN. Audited 5 recent commits: #1558 drop-defensive-copies (frombuffer path still .copy()s before in-place predictor decode at _reader.py:778), #1556 fp-predictor ngjit (writer pre-ravels so 1-D slice arg is correct, float32/64 LE+BE bit-exact), #1552 batched D2H (OOM guard fires before cupy.concatenate, host_buf offsets correct), #1551 parallel-decode gate (>= vs > sends 256x256 default to parallel path, no value diff confirmed via partial-tile parity), #1549 nvjpeg constants (gray + RGB GPU JPEG decode pixel-identical to Pillow CPU, max diff = 0). Cross-backend parity re-verified clean: numpy/dask+numpy/cupy/dask+cupy equal .data/.dtype/.coords/nodata/NaN-mask on deflate+predictor3+nodata; orientations 1-8 numpy==GPU; partial edge tiles 100x150, 257x383, 512x257 numpy==GPU==dask; predictor2 LE/BE round-trip uint8/int16/uint16/int32/uint32 pass; predictor3 LE/BE float32/64 pass. Deferred LOW (pre-existing, not opened): float16 (bps=16, SampleFormat=3) absent from tiff_dtype_to_numpy map - writer never emits, asymmetric but unreachable. | Pass 9 (2026-05-09): TWO HIGH fixed -- (a) PR #1539 closes #1537: TIFF Orientation tag 2/3/4 (mirror flips) on georeferenced files left y/x coords computed from the un-flipped transform, so xarray label lookups returned the wrong pixel even though _apply_orientation flipped the buffer. PR #1521 only updated the transform for the 5-8 axis-swap branch. Fix updates origin and pixel-scale signs along whichever axes were flipped, for both PixelIsArea (origin shifts by N*step) and PixelIsPoint (shifts by (N-1)*step). 10 new tests in test_orientation.py. (b) PR #1546 closes #1540: read_geotiff_gpu ignored Orientation tag completely; CPU correctly applied 2-8 (PR #1521) but GPU returned the raw stored buffer. Cross-backend disagreement on every non-default orientation. Fix adds _apply_orientation_gpu (cupy slicing mirror of the CPU helper) and _apply_orientation_geo_info, threads them into the tiled GPU pipeline, reuses CPU-fallback geo_info for the stripped path to avoid double-applying. 28 new tests in test_orientation_gpu.py (every orientation, single-band tiled, single-band stripped, 3-band tiled, mirror-flip sel-fidelity, default no-tag passthrough). Re-confirmed clean: HTTP coalesce_ranges with overlapping ranges and zero-length ranges, parallel streaming write thread-safety (each tile gets independent buffer via copy or padded zeros), planar=2 + chunky GPU LERC mask propagation matches CPU, IFD chain cap MAX_IFDS=256, max_z_error round-trip on tiled write, _resolve_masked_fill float vs integer dtype semantics. Deferred LOW: per-sample LERC mask (3D mask (h,w,samples)) collapsed to per-pixel ""any sample invalid"" on GPU while CPU honours per-sample; LERC implementations rarely emit 3D masks (verified: lerc.encode with 2D mask on 3-band returns 2D mask). Documented planar=2 + LERC + GPU silently drops mask (rare in practice, source comment acknowledges). | Pass 8 (2026-05-07): HIGH fixed in fix-jpeg-tiff-disable -- to_geotiff(compression='jpeg') wrote files that no external reader can decode. The writer tags compression=7 (new-style JPEG) but emits a self-contained JFIF stream per tile/strip and never writes the JPEGTables tag (347) that the TIFF spec requires for that codec. libtiff/GDAL/rasterio all reject the file with TIFFReadEncodedStrip() failed; our reader round-trips because Pillow decodes the standalone JFIF, hiding the break. Pass-4 notes flagged the read side of the same JPEGTables gap and deferred it; pass-8 covers the write side. Fix: reject compression='jpeg' at the to_geotiff entry with a clear ValueError pointing at deflate/zstd/lzw. The internal _writer.write is untouched so the existing self-decoding tests still cover the codec; re-enabling the public path needs a JPEGTables-aware encoder. PR diffs reviewed but not merged: #1512 (BytesIO source) and #1513 (LERC max_z_error) -- both look correct; #1512 file-like read path goes through read_all() once so the per-call BytesIOSource lock is theoretical, and #1513 forwards max_z_error through every overview/tile/strip/streaming path including _write_vrt_tiled and _compress_block. No regressions found in either open PR. Other surfaces audited clean: predictor=3 with float16 (writer auto-promotes to float32 on both eager and streaming paths, value-exact round-trip); planar=2 multi-tile read uses band_idx*tiles_per_band offset so no cross-contamination between planes; _header.py multi-byte tag parsing uses bo (byte_order) consistently; Pillow YCbCr-vs-tagged-RGB photometric mismatch becomes moot once JPEG is disabled. Deferred (LOW/MEDIUM, not filed): JPEG2000 writer accepts arbitrary dtype with no validation (rare codec, narrow risk); float16 dtype not in tiff_dtype_to_numpy decode map (writer never emits it - asymmetric but unreachable); Orientation tag (274) still ignored on read (pass-4 deferral). | Pass 7 (2026-05-07): HIGH fixed in fix-mmap-cache-refcount-after-replace -- _MmapCache.release() looked up the cache entry by realpath, so a holder that acquired the OLD mmap before an os.replace and released it AFTER another caller had acquired the post-replace entry would decrement the new holder's refcount. Subsequent eviction (cache full, or another acquire) closed the still-in-use mmap, breaking reads with 'mmap closed or invalid'. Real exposure: any concurrent reader/writer pattern where to_geotiff replaces a file that another reader had just opened via open_geotiff with chunks= or via _FileSource. PR #1506 added stale-replacement detection but did not fix the refcount confusion across the pop. Fix: acquire returns an opaque entry token; release takes the token and decrements that exact entry, regardless of cache state. Orphaned (popped) entries close their fh+mmap when their own refcount hits zero. _FileSource updated to pass the token. Regression test test_release_after_path_replacement_does_not_clobber_new_holder added. All 665 geotiff tests pass; GPU path verified. | Pass 6 (2026-05-07) PR #1507: BE pred2 numba TypingError. | Pass 5 (2026-05-06) PR #1506: mmap cache stale after file replace. | Pass 4 (2026-05-06) PR #1501: sparse COG tiles. | Pass 3 (2026-05-06) PR #1500: predictor=3 byte order. | Pass 2 (2026-05-05) PR #1498: predictor=2 sample-wise. | Pass 1 (2026-04-23) PR #1247. Re-confirmed clean over passes 2-7: items 2 (writer always emits LE TIFFs - hardcoded b'II'), 3 (RowsPerStrip default = height when missing), 4 (StripByteCounts missing raises clear ValueError), 5 (TileWidth without TileLength caught by 'tw <= 0 or th <= 0' check at _reader.py:688), 9 (read determinism on compressed+tiled+multiband), 11 (predictor=2 with awkward sample stride round-trips), 18 (compression_level=99 raises ValueError 'out of range for deflate (valid: 1-9)'), 21 (concurrent writes serialize correctly via mkstemp+os.replace), 24 (uint16 dtype preserved on numpy backend, dask honors chunks param), 26 (chunks rounds correctly with remainder chunk for non-tile-aligned). Deferred: item 8 (BytesIO/file-like sources are not supported, source.lower() error) - documented as 'str' parameter, not a bug; item 19 (LERC max_z_error not user-exposed by to_geotiff) - missing feature, not a bug."
+geotiff,2026-05-12,1691,HIGH,1;2;5,"Pass 20 (2026-05-12): HIGH fixed -- PR #1691 (no issue created; agent harness blocked gh issue create). Integer COG overview pyramid mixed sentinel into reduced pixels. _block_reduce_2d (_writer.py:258-264) and _block_reduce_2d_gpu (_gpu_decode.py:3027-3028) promoted integer blocks to float64 but never masked the sentinel to NaN before nanmean / nanmin / nanmax / nanmedian. The reduction averaged the sentinel into surrounding valid cells (e.g. (-9999 + 100 + 100 + 100)/4 = -2425 cast back to int16), producing overview pixels that the read-side int-to-NaN mask in open_geotiff couldn't recover because they didn't equal the sentinel. Silent garbage at every zoom above level 0 for to_geotiff(int_data, cog=True, nodata=N). Methods affected: mean, min, max, median; nearest/mode safe (no averaging). Fix: gate the sentinel-to-NaN mask on representability in the source integer dtype (mirrors _int_nodata_in_range in _reader.py) so uint16+GDAL_NODATA=""-9999"" stays a no-op; rewrite all-sentinel-block NaN back to sentinel before the integer dtype cast so the cast is well-defined (the caller's post-overview loop in write() only runs for floats). GPU mirror gets the same path with cupy.where + cupy.isnan for byte parity with CPU. 38 regression tests in test_cog_int_overview_nodata_2026_05_12.py: _block_reduce_2d per-dtype/per-method matrix (uint8/uint16/int16/int32 x mean/min/max/median), all-sentinel-block, no-nodata regression, out-of-range sentinel no-op, end-to-end uint16 + int16 round-trip, 3-band integer COG, GPU per-dtype/per-method matrix, CPU/GPU byte-match parity. All 1606 existing geotiff tests still pass. Categories: Cat 1 (precision/representation loss in nan-aware reduction) + Cat 2 (silent NaN-equivalent corruption from sentinel poisoning) + Cat 5 (backend parity between float and integer code paths within the same writer). Deferred LOW: HTTP COG path (_read_cog_http at _reader.py:1638) skips the band-range validation that local/dask/GPU added in #1673; band=-1 silently selects the last channel on HTTP while local raises IndexError. Cat 5, MEDIUM-leaning but separate concern from the overview fix; one-finding-per-PR per project policy. | Pass 19 (2026-05-12): MEDIUM fixed -- issue #1655. read_vrt silently dropped 0 on a SimpleSource because of src.nodata or nodata at _vrt.py:370. Python treats 0.0 as falsy, so the per-source sentinel fell through to the band-level (or None when missing) and pixels equal to 0.0 in the source file survived as valid data. The in-code comment acknowledged the quirk as backward compat, but the resulting behaviour silently biased every NaN-aware aggregation on VRT mosaics whose sources used 0 as a sentinel (a common convention for unsigned remote-sensing imagery). Fix: src_nodata = src.nodata if src.nodata is not None else nodata. Five regression tests in test_vrt_source_nodata_zero_1655.py covering source NODATA=0, integer XML literal, non-zero unchanged, band-level NoDataValue=0 still honoured, and source-overrides-band precedence. All 100 vrt-related geotiff tests still pass; 3 pre-existing test_features.py matplotlib palette failures unrelated. Categories: Cat 2 (NaN propagation) + Cat 5 (backend inconsistency: read_geotiff masks 0 correctly when GDAL_NODATA tag is set; only VRT path was broken). | Pass 18 (2026-05-11): MEDIUM fixed -- issue #1642. PR #1641 (issue #1640) inherited level-0 georef on overview reads but kept the level-0 origin_x/origin_y unchanged. That is correct for PixelIsArea (origin = upper-left corner of pixel (0,0)) but wrong for PixelIsPoint (origin = center of pixel (0,0), GeoKey 1025 = 2). For a 1024x1024 PixelIsPoint COG with 10 m pixels and origin (0, 0), open_geotiff(overview_level=1) returned x[:3]=[0,20,40] instead of [5,25,45] (level-1 pixel 0 covers level-0 pixels 0-1 whose centers are 0 and 10, centroid 5); same for y. Downstream sel/interp/reproject silently snaps to the wrong pixel for any DEM-style PixelIsPoint COG (USGS, OpenTopography, Copernicus DEM). Categories: Cat 3 (off-by-one / boundary handling) + Cat 5 (raster_type-dependent backend convention). Fix: in extract_geo_info_with_overview_inheritance (_geotags.py), pick the effective raster_type first (overview-declared if non-default, otherwise inherited from parent), then when it is PixelIsPoint apply origin_shift = (scale - 1) * 0.5 * pixel_size_lvl0 along each axis before building the new GeoTransform. PixelIsArea path is byte-equivalent. 13 regression tests in test_overview_pixel_is_point_1642.py: centroid identity across all 4 backends, transform tuple across all 4 backends, uniform grid step, unit-level helper tests for both raster_types via stubbed extract_geo_info, own-geokeys-not-clobbered path on PixelIsPoint, and a PixelIsArea regression check. All 1397 existing non-network geotiff tests still pass (3 pre-existing matplotlib palette failures unrelated). Deferred LOW: non-power-of-two overview dimensions cause scale = base_w/ov_w to diverge from the true 2^level reduction (writer drops the right/bottom strip via h2=(h//2)*2; for h=1023 a level-1 overview has 511 rows so scale=2.0019 not 2.0). Fix would need to either (a) emit explicit geo tags on overview IFDs from the writer or (b) pass the level number into the inheritance helper; neither is a one-line change and the resulting coord error is sub-pixel of level 0. | Pass 17 (2026-05-11): MEDIUM fixed -- issue #1634. open_geotiff eager path windowed read produced confusing CoordinateValidationError when window extended past source extent. read_to_array clamped the window internally and returned a smaller array, but the eager code path used unclamped window indices for y/x coord generation (xrspatial/geotiff/__init__.py lines 562-572), so the coord array length differed from the data and xarray refused to construct the DataArray. Same bug affected the windowed transform shift in _populate_attrs_from_geo_info. The dask path (read_geotiff_dask) already validated up front since #1561, raising a clear ValueError with the format 'window=... is outside the source extent (HxW) or has non-positive size.' so the two backends diverged on the contract. Fix: validate the window up front in open_geotiff's eager branch via _read_geo_info (metadata-only read, no extra pixel cost) using the exact same condition the dask path uses, raising the same ValueError message format. Reproduction: 10x10 raster + window=(5,5,15,15) on eager raised CoordinateValidationError('conflicting sizes ... length 5 ... length 10'); now raises ValueError('window=(5, 5, 15, 15) is outside the source extent (10x10) or has non-positive size.'). Categories: Cat 3 (off-by-one / boundary handling) + Cat 5 (backend inconsistency). 12 regression tests in test_window_out_of_bounds_1634.py: negative start, past-right-edge, past-bottom-edge, past-both-edges, zero-size, inverted window, full-extent ok, interior subset, edge-aligned, eager-vs-dask parity, message-format parity, issue reproducer. All 1286 existing non-network geotiff tests still pass. | Pass 16 (2026-05-11): HIGH fixed -- issue #1623. to_geotiff(cog=True, overview_resampling='cubic', nodata=) on a float raster with NaN regions produced overview pixels with severe ringing artefacts near nodata borders. Same class of bug as #1613 but for the cubic branch: writer rewrites NaN to the sentinel upstream, then _block_reduce_2d(method=cubic) handed the sentinel-poisoned array straight to scipy.ndimage.zoom(order=3). The cubic spline blended the sentinel (e.g. -9999) into neighbouring cells, producing values like 1133.44, -10290.08 where the data was a constant 100. Repro on 16x16 float32 with a 4x4 NaN corner showed 18 polluted pixels in the 8x8 overview. Fix: when nodata is supplied on a float dtype and the sentinel is found, mask sentinel to NaN, run cubic with prefilter=False so a single NaN cannot poison the entire row/column (default B-spline prefilter is global), then rewrite any NaN in the result back to the sentinel. prefilter=False only fires when a sentinel is present so the non-nodata cubic semantics are unchanged. GPU side: _block_reduce_2d_gpu previously raised on method='cubic'; added a CPU fallback (same pattern as 'mode') so GPU writer produces byte-equivalent overviews. GPU_OVERVIEW_METHODS now includes 'cubic'. 12 regression tests in test_cog_cubic_overview_nodata_1623.py (helper no-ringing, poisoning repro, no-nodata unchanged, end-to-end round-trip, GPU fallback, CPU/GPU byte-match, +/-inf nodata mask, NaN-sentinel no-op, GPU_OVERVIEW_METHODS contract). All 1256 existing geotiff tests still pass (3 pre-existing matplotlib failures unrelated). | Pass 15 (2026-05-11): HIGH fixed -- issue #1613. to_geotiff(cog=True, nodata=) on a float raster with NaN produced a corrupted overview pyramid. The NaN-to-sentinel rewrite in __init__.py:1202 (CPU) and :2852 (GPU write_geotiff_gpu) ran BEFORE _make_overview / make_overview_gpu, so the nan-aware aggregations (np.nanmean/min/max/median, cupy.nanmean/min/max/median) saw the sentinel as a real number and biased every overview pixel. Reproduction with -9999 sentinel produced [[-4998.75,-4997.75],..] where np.nanmean gives [[1.5,3.5],..]. Both CPU and GPU paths affected; backend results matched each other but were both wrong (CAT 2 NaN propagation + CAT 5 documents the parity). Fix: _block_reduce_2d / _block_reduce_2d_gpu accept a nodata kwarg that masks the sentinel back to NaN for float dtypes before the reduction; the writer's overview loop passes nodata in, then rewrites all-sentinel reductions (which surface as NaN from the reducer) back to the sentinel for the on-disk pyramid. 11 regression tests in test_cog_overview_nodata_1613.py (CPU mean / partial-block / min/max/median / no-nodata passthrough / helper kwarg / all-sentinel block / GPU mean / GPU helper / CPU-GPU agreement). All 235 nodata/overview/cog tests still pass. | Pass 14 (2026-05-11): HIGH fixed -- issue #1611. read_vrt(band=None) on a multi-band integer VRT with per-band tags only masks band 0's sentinel. __init__.py lines 2795-2809 in read_vrt apply vrt.bands[0].nodata to the full ndim==3 array; bands 1+ keep their integer sentinels as literal finite values (e.g. 65000 surfaces as 65000.0 after the dtype=float64 cast, not NaN). Float-VRT path masks per-band correctly in _vrt._read_data lines 296-297 + 347-351. PR #1602 fixed the single-band band=N case for issue #1598; the band=None multi-band case is the same class of bug. Repro: 2-band uint16 VRT with NoDataValue 65535 / 65000 returns r.values[1,1,1] == 65000.0 instead of NaN; r.values[1,1,0] is NaN (band 0 sentinel masked). Fix scope: in read_vrt, when band is None, iterate over vrt.bands and mask each arr[..., i] slice against its own (gated by the same _int_nodata_in_range guard PR #1583 introduced). Severity HIGH (Cat 2 NaN propagation + Cat 5 backend inconsistency: identical input semantics produce different masking outcomes based on dtype, with finite garbage values where NaN expected). Fix in PR #1612: walks vrt.bands when band is None and ndim==3, masks each arr[..., i] slice against its own via the refactored _sentinel_for_dtype helper (reuses PR #1583's range guard so out-of-range/non-finite/fractional sentinels are a no-op). attrs['nodata'] still carries band 0's sentinel for band=None reads (documented contract). 7 regression tests in test_vrt_multiband_int_nodata_1611.py: uint16 per-band, int32 negative, mixed presence, dtype preservation when no sentinel hit, out-of-range gating, band=N non-regression, attrs contract. 135 existing vrt/nodata geotiff tests still pass. | Pass 13 (2026-05-11): HIGH fixed -- issue #1599. write_geotiff_gpu (and to_geotiff gpu=True) emitted raw NaN bytes for missing pixels even when nodata= was supplied, while the CPU writer substituted NaN with the sentinel before encoding. xrspatial-only round-trips were unaffected (the reader masks both NaN and the sentinel), but external readers (rasterio/GDAL/QGIS) that mask only on the GDAL_NODATA tag saw NaN pixels as valid data -- rasterio reported 100% valid pixels on a 25-NaN file vs CPU's 25-invalid report. Root cause: __init__.py lines 2579-2587 jumped from shape/dtype resolution straight to compression, missing the equivalent of the CPU writer's NaN-to-sentinel rewrite at to_geotiff line ~1156. Fix: cupy.isnan + masked write on a defensive copy of arr, gated on np_dtype.kind=='f' and not np.isnan(float(nodata)). Caller's CuPy buffer preserved (copy before mutate). 7 regression tests in test_gpu_writer_nan_sentinel_1599.py: substitution lands as sentinel, CPU/GPU byte-equivalent, caller buffer not mutated, no-NaN no-op, NaN sentinel skips substitution, rasterio sees identical invalid count on CPU/GPU, multiband 3D path. All other GPU writer tests still pass (50 passed across band-first, attrs, nodata, dask+cupy, writer, nodata aliases). | Pass 12 (2026-05-11): HIGH fixed -- issue #1581. Reading a uint TIFF with a negative GDAL_NODATA sentinel (e.g. uint16 + -9999) raised OverflowError on every backend because the nodata-mask code did arr.dtype.type(int(nodata)) with no range check. Three identical cast sites in __init__.py (numpy eager, _apply_nodata_mask_gpu, _delayed_read_window) plus _resolve_masked_fill and _sparse_fill_value in _reader.py. Fix: _int_nodata_in_range helper gates the cast; out-of-range sentinels are a no-op for value matching (the file can never contain that value), file dtype is preserved, attrs['nodata'] still surfaces the original sentinel so write round-trips keep the GDAL_NODATA tag intact. Matches rasterio behavior. 8 regression tests in test_nodata_out_of_range_1581.py cover the helper, both eager and dask read paths, in-range sentinel non-regression, and GPU helper (cupy-gated). | Pass 11 (2026-05-10): CLEAN. Audited the one additional commit since pass 10 -- #1559 (PR 1548, Centralise GeoTIFF attrs population across all read backends). Refactor extracts _populate_attrs_from_geo_info helper and routes eager numpy, dask, GPU stripped, GPU tiled read paths through it; before the fix dask only emitted crs/transform/raster_type/nodata while numpy emitted the full attrs set including x/y_resolution, resolution_unit, image_description, extra_samples, GDAL metadata, and the CRS-description fields. No data-path arithmetic touched; only attrs dict population. Windowed origin math (origin_x + c0*pixel_width, origin_y + r0*pixel_height) verified to produce -98.0 / 48.75 origin for window=(10,20,50,70) on a (0.1,-0.125) pixel-size raster, with PixelIsArea half-pixel offset preserved on coord lookups (-97.95, 48.6875). Cross-backend attrs parity re-verified: numpy/dask/cupy all emit identical key set on deflate+predictor3+nodata round-trip (crs, crs_wkt, nodata, transform, x_resolution, y_resolution). Data bit-parity re-verified across numpy/dask/cupy on same payload (np.array_equal with equal_nan=True). test_attrs_parity_1548.py (5 tests), test_reader.py/test_writer.py/test_dask_cupy_combined.py (25 tests), GPU orientation/predictor2-BE/LERC-mask/nodata/byteswap suites (65 tests) all green. No accuracy or backend-divergence findings. | Pass 10 (2026-05-10): CLEAN. Audited 5 recent commits: #1558 drop-defensive-copies (frombuffer path still .copy()s before in-place predictor decode at _reader.py:778), #1556 fp-predictor ngjit (writer pre-ravels so 1-D slice arg is correct, float32/64 LE+BE bit-exact), #1552 batched D2H (OOM guard fires before cupy.concatenate, host_buf offsets correct), #1551 parallel-decode gate (>= vs > sends 256x256 default to parallel path, no value diff confirmed via partial-tile parity), #1549 nvjpeg constants (gray + RGB GPU JPEG decode pixel-identical to Pillow CPU, max diff = 0). Cross-backend parity re-verified clean: numpy/dask+numpy/cupy/dask+cupy equal .data/.dtype/.coords/nodata/NaN-mask on deflate+predictor3+nodata; orientations 1-8 numpy==GPU; partial edge tiles 100x150, 257x383, 512x257 numpy==GPU==dask; predictor2 LE/BE round-trip uint8/int16/uint16/int32/uint32 pass; predictor3 LE/BE float32/64 pass. Deferred LOW (pre-existing, not opened): float16 (bps=16, SampleFormat=3) absent from tiff_dtype_to_numpy map - writer never emits, asymmetric but unreachable. | Pass 9 (2026-05-09): TWO HIGH fixed -- (a) PR #1539 closes #1537: TIFF Orientation tag 2/3/4 (mirror flips) on georeferenced files left y/x coords computed from the un-flipped transform, so xarray label lookups returned the wrong pixel even though _apply_orientation flipped the buffer. PR #1521 only updated the transform for the 5-8 axis-swap branch. Fix updates origin and pixel-scale signs along whichever axes were flipped, for both PixelIsArea (origin shifts by N*step) and PixelIsPoint (shifts by (N-1)*step). 10 new tests in test_orientation.py. (b) PR #1546 closes #1540: read_geotiff_gpu ignored Orientation tag completely; CPU correctly applied 2-8 (PR #1521) but GPU returned the raw stored buffer. Cross-backend disagreement on every non-default orientation. Fix adds _apply_orientation_gpu (cupy slicing mirror of the CPU helper) and _apply_orientation_geo_info, threads them into the tiled GPU pipeline, reuses CPU-fallback geo_info for the stripped path to avoid double-applying. 28 new tests in test_orientation_gpu.py (every orientation, single-band tiled, single-band stripped, 3-band tiled, mirror-flip sel-fidelity, default no-tag passthrough). Re-confirmed clean: HTTP coalesce_ranges with overlapping ranges and zero-length ranges, parallel streaming write thread-safety (each tile gets independent buffer via copy or padded zeros), planar=2 + chunky GPU LERC mask propagation matches CPU, IFD chain cap MAX_IFDS=256, max_z_error round-trip on tiled write, _resolve_masked_fill float vs integer dtype semantics. Deferred LOW: per-sample LERC mask (3D mask (h,w,samples)) collapsed to per-pixel ""any sample invalid"" on GPU while CPU honours per-sample; LERC implementations rarely emit 3D masks (verified: lerc.encode with 2D mask on 3-band returns 2D mask). Documented planar=2 + LERC + GPU silently drops mask (rare in practice, source comment acknowledges). | Pass 8 (2026-05-07): HIGH fixed in fix-jpeg-tiff-disable -- to_geotiff(compression='jpeg') wrote files that no external reader can decode. The writer tags compression=7 (new-style JPEG) but emits a self-contained JFIF stream per tile/strip and never writes the JPEGTables tag (347) that the TIFF spec requires for that codec. libtiff/GDAL/rasterio all reject the file with TIFFReadEncodedStrip() failed; our reader round-trips because Pillow decodes the standalone JFIF, hiding the break. Pass-4 notes flagged the read side of the same JPEGTables gap and deferred it; pass-8 covers the write side. Fix: reject compression='jpeg' at the to_geotiff entry with a clear ValueError pointing at deflate/zstd/lzw. The internal _writer.write is untouched so the existing self-decoding tests still cover the codec; re-enabling the public path needs a JPEGTables-aware encoder. PR diffs reviewed but not merged: #1512 (BytesIO source) and #1513 (LERC max_z_error) -- both look correct; #1512 file-like read path goes through read_all() once so the per-call BytesIOSource lock is theoretical, and #1513 forwards max_z_error through every overview/tile/strip/streaming path including _write_vrt_tiled and _compress_block. No regressions found in either open PR. Other surfaces audited clean: predictor=3 with float16 (writer auto-promotes to float32 on both eager and streaming paths, value-exact round-trip); planar=2 multi-tile read uses band_idx*tiles_per_band offset so no cross-contamination between planes; _header.py multi-byte tag parsing uses bo (byte_order) consistently; Pillow YCbCr-vs-tagged-RGB photometric mismatch becomes moot once JPEG is disabled. Deferred (LOW/MEDIUM, not filed): JPEG2000 writer accepts arbitrary dtype with no validation (rare codec, narrow risk); float16 dtype not in tiff_dtype_to_numpy decode map (writer never emits it - asymmetric but unreachable); Orientation tag (274) still ignored on read (pass-4 deferral). | Pass 7 (2026-05-07): HIGH fixed in fix-mmap-cache-refcount-after-replace -- _MmapCache.release() looked up the cache entry by realpath, so a holder that acquired the OLD mmap before an os.replace and released it AFTER another caller had acquired the post-replace entry would decrement the new holder's refcount. Subsequent eviction (cache full, or another acquire) closed the still-in-use mmap, breaking reads with 'mmap closed or invalid'. Real exposure: any concurrent reader/writer pattern where to_geotiff replaces a file that another reader had just opened via open_geotiff with chunks= or via _FileSource. PR #1506 added stale-replacement detection but did not fix the refcount confusion across the pop. Fix: acquire returns an opaque entry token; release takes the token and decrements that exact entry, regardless of cache state. Orphaned (popped) entries close their fh+mmap when their own refcount hits zero. _FileSource updated to pass the token. Regression test test_release_after_path_replacement_does_not_clobber_new_holder added. All 665 geotiff tests pass; GPU path verified. | Pass 6 (2026-05-07) PR #1507: BE pred2 numba TypingError. | Pass 5 (2026-05-06) PR #1506: mmap cache stale after file replace. | Pass 4 (2026-05-06) PR #1501: sparse COG tiles. | Pass 3 (2026-05-06) PR #1500: predictor=3 byte order. | Pass 2 (2026-05-05) PR #1498: predictor=2 sample-wise. | Pass 1 (2026-04-23) PR #1247. Re-confirmed clean over passes 2-7: items 2 (writer always emits LE TIFFs - hardcoded b'II'), 3 (RowsPerStrip default = height when missing), 4 (StripByteCounts missing raises clear ValueError), 5 (TileWidth without TileLength caught by 'tw <= 0 or th <= 0' check at _reader.py:688), 9 (read determinism on compressed+tiled+multiband), 11 (predictor=2 with awkward sample stride round-trips), 18 (compression_level=99 raises ValueError 'out of range for deflate (valid: 1-9)'), 21 (concurrent writes serialize correctly via mkstemp+os.replace), 24 (uint16 dtype preserved on numpy backend, dask honors chunks param), 26 (chunks rounds correctly with remainder chunk for non-tile-aligned). Deferred: item 8 (BytesIO/file-like sources are not supported, source.lower() error) - documented as 'str' parameter, not a bug; item 19 (LERC max_z_error not user-exposed by to_geotiff) - missing feature, not a bug."
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/_gpu_decode.py b/xrspatial/geotiff/_gpu_decode.py
index e3a0f8055..c18c5dbad 100644
--- a/xrspatial/geotiff/_gpu_decode.py
+++ b/xrspatial/geotiff/_gpu_decode.py
@@ -3026,6 +3026,25 @@ def _block_reduce_2d_gpu(arr2d, method, nodata=None):
mask, cupy.float64('nan'), blocks)
else:
blocks = cropped.astype(cupy.float64).reshape(oh, 2, ow, 2)
+ # Integer GPU mirror of the CPU sentinel-mask: without it,
+ # cupy.nanmean / nanmin / nanmax / nanmedian average the sentinel
+ # value into surrounding valid cells and produce overview pixels
+ # that aren't masked on the read side because they don't equal
+ # the sentinel. Same root cause as the CPU bug fixed alongside
+ # this; the GPU writer needs byte parity with CPU (the contract
+ # from #1623).
+ if (nodata is not None
+ and np.isfinite(nodata)
+ and float(nodata).is_integer()):
+ nodata_int = int(nodata)
+ info = np.iinfo(arr2d.dtype)
+ if info.min <= nodata_int <= info.max:
+ sentinel = np.dtype(str(arr2d.dtype)).type(nodata_int)
+ int_blocks = cropped.reshape(oh, 2, ow, 2)
+ mask = int_blocks == sentinel
+ if bool(mask.any().item()):
+ blocks = cupy.where(
+ mask, cupy.float64('nan'), blocks)
if method == 'mean':
result = cupy.nanmean(blocks, axis=(1, 3))
@@ -3042,6 +3061,20 @@ def _block_reduce_2d_gpu(arr2d, method, nodata=None):
f"Use one of: {GPU_OVERVIEW_METHODS}")
if arr2d.dtype.kind != 'f':
+ # All-sentinel 2x2 blocks come back as NaN from cupy.nan*;
+ # casting NaN to integer is undefined on the GPU, so rewrite
+ # those back to the sentinel before the cast. Mirrors the CPU
+ # fallback for the same reason; matches the post-cast invariant
+ # the CPU writer's overview loop relies on.
+ if nodata is not None and np.isfinite(nodata):
+ nan_mask = cupy.isnan(result)
+ if bool(nan_mask.any().item()):
+ nodata_int = int(nodata) if float(nodata).is_integer() else None
+ if nodata_int is not None:
+ info = np.iinfo(arr2d.dtype)
+ if info.min <= nodata_int <= info.max:
+ result = cupy.where(
+ nan_mask, cupy.float64(nodata_int), result)
return cupy.around(result).astype(arr2d.dtype)
return result.astype(arr2d.dtype)
@@ -3059,11 +3092,13 @@ def make_overview_gpu(arr, method='mean', nodata=None):
implementation in :mod:`xrspatial.geotiff._writer` so the GPU
writer path produces the same overview bytes as the CPU writer.
nodata : scalar or None
- When supplied and ``arr`` is a float dtype, cells equal to the
- sentinel are masked back to NaN before the reduction so the
- sentinel does not bias the result. Required for COG output that
- sets ``nodata=...`` (issue #1613, extended to ``cubic`` in
- issue #1623). Ignored for integer arrays and for ``nearest``.
+ When supplied, cells equal to the sentinel are masked back to
+ NaN before the reduction so the sentinel does not bias the
+ result. Applies to float dtypes (issue #1613, extended to
+ ``cubic`` in #1623) and to integer dtypes (covers the
+ sentinel-poisoning case that mixed sentinel + valid pixels in
+ the nan-aware reduction). Ignored for ``nearest`` (no reduction
+ averaging occurs).
Returns
-------
diff --git a/xrspatial/geotiff/_writer.py b/xrspatial/geotiff/_writer.py
index 116f6228e..7758adcb7 100644
--- a/xrspatial/geotiff/_writer.py
+++ b/xrspatial/geotiff/_writer.py
@@ -153,17 +153,24 @@ def _compression_tag(compression_name: str) -> int:
def _block_reduce_2d(arr2d, method, nodata=None):
"""2x block-reduce a single 2D plane using *method*.
- When ``nodata`` is supplied and ``arr2d`` is a float dtype, cells that
- equal the sentinel are treated as NaN during the reduction so the
- ``nan*`` aggregation routines correctly skip them. For the nan-aware
- aggregation methods, the reduced output keeps NaN wherever every
- contributing input cell was the sentinel (so callers can rewrite
- that NaN back to the sentinel after the reduction). The sentinel is
- ignored entirely for integer dtypes and for ``nearest`` and ``mode``
- methods. The ``cubic`` branch honours ``nodata`` by masking the
- sentinel to NaN, running cubic with ``prefilter=False`` to keep the
- kernel local, and rewriting any NaN in the output back to the
- sentinel before returning (issue #1623).
+ When ``nodata`` is supplied, cells that equal the sentinel are
+ treated as NaN during the reduction so the ``nan*`` aggregation
+ routines correctly skip them. The float branch keeps any all-
+ sentinel block as NaN so the caller's post-overview loop can
+ rewrite it back to the sentinel; the integer branch rewrites NaN
+ back to the sentinel before the dtype cast so the cast is
+ well-defined (the caller's post-overview loop only handles the
+ float case). The ``nearest`` and ``mode`` methods do NOT mask the
+ sentinel: ``nearest`` returns the top-left pixel of each 2x2 block
+ and ``mode`` returns the most-frequent value, so the sentinel can
+ be selected as the overview pixel if it occupies that position
+ (``nearest``) or is the most frequent value in the block
+ (``mode``). Mean / median / min / max / cubic all mask the
+ sentinel before reduction. The ``cubic`` branch honours ``nodata``
+ by masking the sentinel to NaN, running cubic with
+ ``prefilter=False`` to keep the kernel local, and rewriting any
+ NaN in the output back to the sentinel before returning (issue
+ #1623).
"""
h, w = arr2d.shape
h2 = (h // 2) * 2
@@ -257,11 +264,34 @@ def _block_reduce_2d(arr2d, method, nodata=None):
blocks = np.where(mask, np.float64('nan'), blocks)
else:
blocks = cropped.astype(np.float64).reshape(oh, 2, ow, 2)
- # Integer rasters can also carry a sentinel that an upstream
- # promotion already converted to NaN; cropped is integer so no
- # masking is needed here. The blocks.astype(float64) cast above
- # would lose any NaN anyway -- integer sentinels are handled at
- # the call site by promoting to float64 before reduction.
+ # Integer rasters with a sentinel need the same NaN-mask the float
+ # branch above applies: without it, nanmean / nanmin / nanmax /
+ # nanmedian average the sentinel value into surrounding valid
+ # cells and produce overview pixels that are neither the sentinel
+ # nor any real measurement. The read-side int-to-NaN mask in
+ # ``open_geotiff`` only catches exact sentinel hits, so the
+ # poisoned values survive as silent garbage at every zoom level
+ # above 0. Gate on the sentinel being representable in the
+ # source integer dtype (mirrors ``_int_nodata_in_range`` in
+ # ``_reader.py``) so an out-of-range sentinel pair like
+ # ``uint16`` + ``GDAL_NODATA="-9999"`` stays a no-op rather than
+ # tripping ``OverflowError`` on the dtype cast.
+ if (nodata is not None
+ and np.isfinite(nodata)
+ and float(nodata).is_integer()):
+ nodata_int = int(nodata)
+ info = np.iinfo(arr2d.dtype)
+ if info.min <= nodata_int <= info.max:
+ sentinel = arr2d.dtype.type(nodata_int)
+ # Compare against the original integer block view so
+ # the equality runs at the integer's native width
+ # (avoids any float-cast rounding on adjacent values).
+ # The boolean mask broadcasts into the float64 block
+ # layout below.
+ int_blocks = cropped.reshape(oh, 2, ow, 2)
+ mask = int_blocks == sentinel
+ if mask.any():
+ blocks = np.where(mask, np.float64('nan'), blocks)
# nanmean / nanmin / nanmax / nanmedian emit RuntimeWarning when a
# 2x2 block is all-NaN (typical at nodata borders). The all-NaN
@@ -284,6 +314,23 @@ def _block_reduce_2d(arr2d, method, nodata=None):
f"Use one of: {OVERVIEW_METHODS}")
if arr2d.dtype.kind != 'f':
+ # All-sentinel 2x2 blocks come back as NaN from the nan-aware
+ # reduction; cast NaN to an integer dtype is undefined (varies
+ # between platforms / produces zero or INT_MIN). Rewrite those
+ # back to the sentinel before the cast so the integer overview
+ # pyramid carries the same masking convention as the
+ # full-resolution band. The float branch relies on the caller's
+ # post-overview rewrite in ``write()``; integer dtypes skip that
+ # branch because ``current.dtype.kind == 'f'`` is False, so we
+ # close the loop here.
+ if nodata is not None and np.isfinite(nodata):
+ nan_mask = np.isnan(result)
+ if nan_mask.any():
+ info = np.iinfo(arr2d.dtype)
+ nodata_int = int(nodata) if float(nodata).is_integer() else None
+ if (nodata_int is not None
+ and info.min <= nodata_int <= info.max):
+ result = np.where(nan_mask, float(nodata_int), result)
return np.round(result).astype(arr2d.dtype)
return result.astype(arr2d.dtype)
@@ -300,12 +347,14 @@ def _make_overview(arr: np.ndarray, method: str = 'mean',
Resampling method: 'mean' (default), 'nearest', 'min', 'max',
'median', 'mode', or 'cubic'.
nodata : scalar or None
- When supplied and ``arr`` is a float dtype, cells equal to the
- sentinel are masked back to NaN before the reduction so the
- sentinel does not bias the result. Required for COG output that
- sets ``nodata=...`` (issue #1613, extended to ``cubic`` in
- issue #1623). Ignored for integer arrays and for ``nearest`` /
- ``mode`` methods.
+ When supplied, cells equal to the sentinel are masked back to
+ NaN before the reduction so the sentinel does not bias the
+ result. Applies to both float dtypes (issue #1613, extended to
+ ``cubic`` in #1623) and integer dtypes (the mean / min / max /
+ median reductions used to average the sentinel into surrounding
+ valid pixels and produce overview values that the reader could
+ not mask). Ignored for ``nearest`` / ``mode`` methods (no
+ averaging occurs).
Returns
-------
diff --git a/xrspatial/geotiff/tests/test_cog_int_overview_nodata_2026_05_12.py b/xrspatial/geotiff/tests/test_cog_int_overview_nodata_2026_05_12.py
new file mode 100644
index 000000000..1ad627786
--- /dev/null
+++ b/xrspatial/geotiff/tests/test_cog_int_overview_nodata_2026_05_12.py
@@ -0,0 +1,315 @@
+"""COG overview generation respects the nodata sentinel for integer rasters.
+
+Companion to issue #1613 (float COG overview poisoning). Before this fix,
+``to_geotiff(int_data, cog=True, nodata=N)`` ran the overview reduction with
+the sentinel still present in the integer-cast float64 block. The nan-aware
+reduction (``np.nanmean`` / nanmin / nanmax / nanmedian) averaged the
+sentinel into surrounding valid pixels and produced overview values that
+the reader could not mask -- they did not equal the sentinel, so the
+int-to-NaN mask in ``open_geotiff`` left them as silent garbage.
+
+These tests pin the contract that the CPU writer (and the GPU mirror in
+``_block_reduce_2d_gpu``) skip the integer sentinel during overview
+reduction, so the resulting pyramid only contains real measurements and
+the sentinel value.
+"""
+from __future__ import annotations
+
+import importlib.util
+
+import numpy as np
+import pytest
+import xarray as xr
+
+
+def _gpu_available() -> bool:
+ if importlib.util.find_spec("cupy") is None:
+ return False
+ try:
+ import cupy
+ return bool(cupy.cuda.is_available())
+ except Exception:
+ return False
+
+
+_HAS_GPU = _gpu_available()
+_gpu_only = pytest.mark.skipif(
+ not _HAS_GPU,
+ reason="cupy + CUDA required",
+)
+
+
+# ---------------------------------------------------------------------------
+# Unit-level: _block_reduce_2d on integer dtypes
+# ---------------------------------------------------------------------------
+
+def _int_block_partial_sentinel(sentinel, dtype):
+ """4x4 integer raster where the right two columns of each row pair
+ mix valid and sentinel cells. Block (0, 1) has (100, 100, sentinel,
+ sentinel); block (1, 1) has (200, 200, sentinel, sentinel)."""
+ arr = np.array([
+ [100, 100, 100, 100],
+ [100, 100, sentinel, sentinel],
+ [200, 200, 200, 200],
+ [200, 200, sentinel, sentinel],
+ ], dtype=dtype)
+ return arr
+
+
+@pytest.mark.parametrize('method', ['mean', 'min', 'max', 'median'])
+@pytest.mark.parametrize('dtype,sentinel', [
+ (np.uint8, 255),
+ (np.uint16, 65535),
+ (np.int16, -9999),
+ (np.int32, -2_000_000_000),
+])
+def test_block_reduce_int_sentinel_masked(method, dtype, sentinel):
+ """Integer overview reductions must skip sentinel cells.
+
+ Before the fix, mean produced averages like ``(100+sentinel)/2`` cast
+ back to the integer dtype -- a non-sentinel value that the reader
+ leaves untouched. The fix masks the sentinel to NaN before the
+ reduction so nan-aware aggregation skips it.
+ """
+ from xrspatial.geotiff._writer import _block_reduce_2d
+
+ arr = _int_block_partial_sentinel(sentinel, dtype)
+ out = _block_reduce_2d(arr, method, nodata=sentinel)
+
+ # Every block now has at least one valid 100/200; result should equal
+ # the valid value (since for mean/min/max/median over {100, 100} is
+ # 100, and over {200, 200} is 200). Neither block has any cell that
+ # isn't 100, 200, or sentinel, so the output must be a subset of
+ # {100, 200}.
+ assert out.dtype == arr.dtype
+ out_vals = set(out.flatten().tolist())
+ assert out_vals.issubset({100, 200}), (
+ f"method={method} dtype={dtype} sentinel={sentinel} "
+ f"produced poisoned values: {out_vals - {100, 200}}"
+ )
+
+
+@pytest.mark.parametrize('dtype,sentinel', [
+ (np.uint16, 65535),
+ (np.int16, -9999),
+])
+def test_block_reduce_int_all_sentinel_block(dtype, sentinel):
+ """A 2x2 block that's entirely sentinel reduces to the sentinel.
+
+ Without the post-reduction NaN-to-sentinel rewrite in the integer
+ branch, the all-NaN block from nanmean would cast to undefined
+ integer behaviour (zero or INT_MIN depending on platform).
+ """
+ from xrspatial.geotiff._writer import _block_reduce_2d
+
+ arr = np.array([
+ [100, 100, sentinel, sentinel],
+ [100, 100, sentinel, sentinel],
+ [200, 200, 200, 200],
+ [200, 200, 200, 200],
+ ], dtype=dtype)
+
+ out = _block_reduce_2d(arr, 'mean', nodata=sentinel)
+ assert out.dtype == arr.dtype
+ # Top-right block is all-sentinel; output must be the sentinel
+ assert out[0, 1] == sentinel
+ # Other blocks contain only valid values
+ assert out[0, 0] == 100
+ assert out[1, 0] == 200
+ assert out[1, 1] == 200
+
+
+def test_block_reduce_int_no_nodata_unchanged():
+ """Without ``nodata``, the integer reduction code path stays unchanged.
+
+ Regression check: the fix must not alter the no-sentinel case.
+ """
+ from xrspatial.geotiff._writer import _block_reduce_2d
+
+ arr = np.array([
+ [1, 2, 3, 4],
+ [5, 6, 7, 8],
+ [9, 10, 11, 12],
+ [13, 14, 15, 16],
+ ], dtype=np.int16)
+
+ out = _block_reduce_2d(arr, 'mean')
+ # Block (0,0) = mean(1,2,5,6) = 3.5 -> round -> 4
+ # Block (0,1) = mean(3,4,7,8) = 5.5 -> round -> 6
+ # Block (1,0) = mean(9,10,13,14) = 11.5 -> round -> 12
+ # Block (1,1) = mean(11,12,15,16) = 13.5 -> round -> 14
+ expected = np.array([[4, 6], [12, 14]], dtype=np.int16)
+ np.testing.assert_array_equal(out, expected)
+
+
+def test_block_reduce_int_out_of_range_sentinel_noop():
+ """A sentinel outside the dtype's range is a no-op (no mask applied).
+
+ Mirrors the ``_int_nodata_in_range`` gating in ``_reader.py``: a
+ uint16 file with GDAL_NODATA="-9999" cannot match any decoded pixel,
+ so the reduction proceeds without the mask. This keeps the fix from
+ raising OverflowError on the dtype cast.
+ """
+ from xrspatial.geotiff._writer import _block_reduce_2d
+
+ # uint16 with nodata=-9999: out of range, no-op
+ arr = np.array([
+ [1, 2, 3, 4],
+ [5, 6, 7, 8],
+ ], dtype=np.uint16)
+ out = _block_reduce_2d(arr, 'mean', nodata=-9999)
+ # Should produce the same result as without the kwarg
+ expected = _block_reduce_2d(arr, 'mean')
+ np.testing.assert_array_equal(out, expected)
+
+
+# ---------------------------------------------------------------------------
+# End-to-end: to_geotiff + open_geotiff round trip
+# ---------------------------------------------------------------------------
+
+@pytest.fixture
+def _int_cog_inputs(tmp_path):
+ """uint16 raster, full of 100 with a 65x65 sentinel patch."""
+ H, W = 256, 256
+ data = np.full((H, W), 100, dtype=np.uint16)
+ data[64:129, 64:129] = 65535
+ da = xr.DataArray(
+ data,
+ dims=('y', 'x'),
+ coords={'y': np.arange(H, dtype=np.float64),
+ 'x': np.arange(W, dtype=np.float64)},
+ attrs={'crs': 4326},
+ )
+ return da, tmp_path
+
+
+@pytest.mark.parametrize('method', ['mean', 'min', 'max', 'median'])
+def test_cpu_int_cog_overview_not_poisoned(_int_cog_inputs, method):
+ """End-to-end: integer COG overview pyramid contains only valid values.
+
+ Before the fix, the level-1 read contained values like 16459 and
+ 32818 -- nan-aware-mean of (sentinel, 100, 100, 100) and (sentinel,
+ sentinel, 100, 100) cast back to uint16. The reader can't mask them
+ because they don't equal 65535.
+ """
+ from xrspatial.geotiff import to_geotiff, open_geotiff
+
+ da, tmp_path = _int_cog_inputs
+ path = str(tmp_path / f'int_overview_{method}_2026_05_12.tif')
+ to_geotiff(da, path, nodata=65535, cog=True,
+ overview_levels=[1], overview_resampling=method)
+
+ ov = open_geotiff(path, overview_level=1)
+ arr = np.asarray(ov.data)
+ unique = set(int(v) for v in np.unique(arr) if not np.isnan(v))
+ poisoned = unique - {100, 65535}
+ assert not poisoned, (
+ f"method={method} produced poisoned overview values: {poisoned}"
+ )
+
+
+def test_cpu_int_cog_overview_3band_not_poisoned(tmp_path):
+ """3-band integer COG: same fix applies via the 3D _make_overview branch."""
+ from xrspatial.geotiff import to_geotiff, open_geotiff
+
+ H, W = 256, 256
+ data = np.full((H, W, 3), 100, dtype=np.uint16)
+ data[64:129, 64:129, :] = 65535
+ da = xr.DataArray(
+ data,
+ dims=('y', 'x', 'band'),
+ coords={'y': np.arange(H, dtype=np.float64),
+ 'x': np.arange(W, dtype=np.float64),
+ 'band': [0, 1, 2]},
+ attrs={'crs': 4326},
+ )
+
+ path = str(tmp_path / 'int_overview_3band_2026_05_12.tif')
+ to_geotiff(da, path, nodata=65535, cog=True,
+ overview_levels=[1], overview_resampling='mean')
+
+ ov = open_geotiff(path, overview_level=1)
+ arr = np.asarray(ov.data)
+ unique = set(int(v) for v in np.unique(arr) if not np.isnan(v))
+ poisoned = unique - {100, 65535}
+ assert not poisoned, (
+ f"3-band integer overview produced poisoned values: {poisoned}"
+ )
+
+
+def test_cpu_int_cog_no_nodata_unchanged(tmp_path):
+ """No nodata kwarg: integer overview path stays as it was."""
+ from xrspatial.geotiff import to_geotiff, open_geotiff
+
+ H, W = 256, 256
+ data = np.full((H, W), 100, dtype=np.uint16)
+ data[100:200, 100:200] = 50
+ da = xr.DataArray(
+ data,
+ dims=('y', 'x'),
+ coords={'y': np.arange(H, dtype=np.float64),
+ 'x': np.arange(W, dtype=np.float64)},
+ attrs={'crs': 4326},
+ )
+
+ path = str(tmp_path / 'int_overview_no_nodata_2026_05_12.tif')
+ to_geotiff(da, path, cog=True,
+ overview_levels=[1], overview_resampling='mean')
+
+ ov = open_geotiff(path, overview_level=1)
+ arr = np.asarray(ov.data)
+ # No sentinel, so every overview pixel is a real average of 50 / 100.
+ # Block-boundary pixels are weighted means: (50,50,50,100)/4 = 62.5 -> 63
+ unique = set(int(v) for v in np.unique(arr))
+ # Must contain at least 50 and 100; boundary-mixing averages allowed.
+ assert 50 in unique
+ assert 100 in unique
+
+
+# ---------------------------------------------------------------------------
+# GPU mirror
+# ---------------------------------------------------------------------------
+
+@_gpu_only
+@pytest.mark.parametrize('method', ['mean', 'min', 'max', 'median'])
+@pytest.mark.parametrize('dtype,sentinel', [
+ (np.uint16, 65535),
+ (np.int16, -9999),
+])
+def test_gpu_block_reduce_int_sentinel_masked(method, dtype, sentinel):
+ """GPU mirror of the CPU integer sentinel-mask fix."""
+ import cupy
+ from xrspatial.geotiff._gpu_decode import _block_reduce_2d_gpu
+
+ arr = _int_block_partial_sentinel(sentinel, dtype)
+ cpu_arr = cupy.asarray(arr)
+ out_gpu = _block_reduce_2d_gpu(cpu_arr, method, nodata=sentinel)
+ out = out_gpu.get()
+
+ assert out.dtype == arr.dtype
+ out_vals = set(out.flatten().tolist())
+ assert out_vals.issubset({100, 200}), (
+ f"GPU method={method} dtype={dtype} produced poisoned values: "
+ f"{out_vals - {100, 200}}"
+ )
+
+
+@_gpu_only
+@pytest.mark.parametrize('method', ['mean', 'min', 'max', 'median'])
+def test_gpu_cpu_int_overview_byte_match(method):
+ """CPU and GPU integer overview reductions agree byte-for-byte.
+
+ Same parity contract as #1623 (cubic). Without the GPU fix, the GPU
+ pyramid would carry poisoned values while the CPU pyramid carried
+ sentinels -- two backends disagreeing on identical input.
+ """
+ import cupy
+ from xrspatial.geotiff._writer import _block_reduce_2d
+ from xrspatial.geotiff._gpu_decode import _block_reduce_2d_gpu
+
+ arr = _int_block_partial_sentinel(-9999, np.int16)
+ cpu_out = _block_reduce_2d(arr, method, nodata=-9999)
+ gpu_out = _block_reduce_2d_gpu(
+ cupy.asarray(arr), method, nodata=-9999).get()
+
+ np.testing.assert_array_equal(cpu_out, gpu_out)