diff --git a/xrspatial/geotiff/_gpu_decode.py b/xrspatial/geotiff/_gpu_decode.py index 35f7bb835..3e9e8c66c 100644 --- a/xrspatial/geotiff/_gpu_decode.py +++ b/xrspatial/geotiff/_gpu_decode.py @@ -3151,7 +3151,7 @@ def _block_reduce_2d_gpu(arr2d, method, nodata=None): Odd inputs are NaN-padded along the trailing edge (float-promoted for integer dtypes) so the 2x2 block reshape works and the residual block is reduced via the same nan-aware aggregations (issue #2105). Mirrors - the CPU helper :func:`xrspatial.geotiff._writer._block_reduce_2d` so + the CPU helper :func:`xrspatial.geotiff._overview._block_reduce_2d` so the two backends produce identical overviews. When ``nodata`` is supplied and ``arr2d`` is a float dtype, cells that @@ -3178,7 +3178,7 @@ def _block_reduce_2d_gpu(arr2d, method, nodata=None): # Mode is expensive on GPU; fall back to CPU. The CPU helper # now handles odd-sized inputs natively. cpu_arr = arr2d.get() - from ._writer import _block_reduce_2d + from ._overview import _block_reduce_2d cpu_result = _block_reduce_2d(cpu_arr, 'mode', nodata=nodata) return cupy.asarray(cpu_result) @@ -3190,7 +3190,7 @@ def _block_reduce_2d_gpu(arr2d, method, nodata=None): # CPU writer and so the sentinel handling matches. The CPU # helper handles odd-sized inputs via edge-replicate padding. cpu_arr = arr2d.get() - from ._writer import _block_reduce_2d + from ._overview import _block_reduce_2d cpu_result = _block_reduce_2d(cpu_arr, 'cubic', nodata=nodata) return cupy.asarray(cpu_result) diff --git a/xrspatial/geotiff/_overview.py b/xrspatial/geotiff/_overview.py new file mode 100644 index 000000000..6360e5dee --- /dev/null +++ b/xrspatial/geotiff/_overview.py @@ -0,0 +1,474 @@ +"""Overview (COG pyramid) helpers extracted from :mod:`_writer`. + +This module is private to :mod:`xrspatial.geotiff`. It hosts the +CPU-side helpers that build the reduced-resolution overview IFDs for +Cloud-Optimized GeoTIFFs: the level-list validator, the 2x block +reducer, its replicate-pad helper, the integer-nodata resolver, and +the per-array overview entry point. + +Behavior is unchanged versus the original in ``_writer.py``; the +extraction is mechanical so the strip/tile / IFD / layout code in +``_writer.py`` no longer carries pyramid logic in the same file. +External callers (the ``_writers`` subpackage, GPU fallbacks, tests) +keep working via re-exports from ``_writer``. +""" +from __future__ import annotations + +import warnings + +import numpy as np + + +OVERVIEW_METHODS = ('mean', 'nearest', 'min', 'max', 'median', 'mode', 'cubic') + +#: Maximum number of overview levels generated by auto-overview mode in COG +#: writes. 8 halvings = 1/256 of the original resolution, which is enough +#: for any practical raster. Pass ``overview_levels=[...]`` explicitly to +#: override. +_MAX_OVERVIEW_LEVELS = 8 + + +def _validate_overview_levels(overview_levels, height=None, width=None): + """Validate and normalise an explicit ``overview_levels`` list. + + Each entry is a decimation factor relative to full resolution. + Factors must be strictly increasing integers >= 2, and each must + be a power of two because the underlying block reducer only does + 2x decimation per step (issue #1766 — prior to that fix, the values + were ignored and only the list length mattered). + + When ``height`` and ``width`` are supplied, each factor is also + checked against the input shape: a factor F is feasible only if + ``height // F >= 1`` and ``width // F >= 1``. Asking for a factor + that would reduce the source below 1 pixel in either dimension + raises ``ValueError`` instead of silently writing a zero-sized + overview IFD. + + Returns the validated list of ints. ``None`` passes through so the + caller can run its auto-generation path. + """ + if overview_levels is None: + return None + if not isinstance(overview_levels, (list, tuple)): + raise ValueError( + f"overview_levels must be a list or tuple of ints, got " + f"{type(overview_levels).__name__}.") + if len(overview_levels) == 0: + return [] + cleaned = [] + prev = 1 + for i, level in enumerate(overview_levels): + # Reject bools explicitly; ``bool`` is an ``int`` subclass and + # ``True``/``False`` would otherwise sneak through the integer + # check below. + if isinstance(level, bool) or not isinstance(level, (int, np.integer)): + raise ValueError( + f"overview_levels[{i}] must be an int >= 2, got " + f"{level!r} (type {type(level).__name__}).") + level = int(level) + if level < 2: + raise ValueError( + f"overview_levels[{i}] must be >= 2 (1 is the original " + f"full-resolution band), got {level}.") + if level <= prev: + raise ValueError( + f"overview_levels must be strictly increasing, got " + f"{list(overview_levels)} (entry at index {i} is " + f"{level}, previous was {prev}).") + # Power-of-two check: the underlying ``_make_overview`` only + # halves, so reaching factor N takes log2(N) halvings and N + # must be a power of two for the cumulative factor to land on + # the requested value exactly. + if (level & (level - 1)) != 0: + raise ValueError( + f"overview_levels[{i}]={level} is not a power of two. " + f"Only power-of-two decimation factors are supported " + f"(2, 4, 8, 16, ...).") + # Shape feasibility: refuse factors that would shrink the + # raster below 1 pixel. ``_block_reduce_2d`` halves via + # ``(dim // 2) * 2`` which produces a zero-sized array once + # ``dim < 2``, and chaining further halvings keeps it at zero. + # Without this check the writer silently emits zero-sized + # overview IFDs. + if height is not None and width is not None: + if height // level < 1 or width // level < 1: + raise ValueError( + f"overview_levels[{i}]={level} is too large for " + f"input shape ({height}, {width}); decimation " + f"would produce a zero-sized overview.") + cleaned.append(level) + prev = level + return cleaned + + +def _resolve_int_nodata(dtype, nodata): + """Return ``int(nodata)`` if it is representable as *dtype*, else None. + + Folds the three checks that gate the integer sentinel-to-NaN mask in + :func:`_block_reduce_2d` into one call: ``dtype`` is integer, the + sentinel is finite and integer-valued, and the integer fits the + dtype range. Out-of-range pairs like ``uint16`` + ``GDAL_NODATA=-9999`` + return None so the caller stays a no-op rather than tripping + ``OverflowError`` on the dtype cast. Mirrors + ``_int_nodata_in_range`` in ``_reader.py``. + """ + if nodata is None or dtype.kind not in ('i', 'u'): + return None + if not np.isfinite(nodata) or not float(nodata).is_integer(): + return None + nodata_int = int(nodata) + info = np.iinfo(dtype) + if info.min <= nodata_int <= info.max: + return nodata_int + return None + + +def _replicate_pad_2d(src, h2, w2): + """Pad ``src`` to ``(h2, w2)`` by replicating its trailing row/col. + + Used by the cubic overview branch to extend an odd-sized input to + even dimensions so a 0.5-zoom lands on the GDAL ceil shape and the + spline kernel sees a valid neighbourhood at the trailing edge. The + helper is dtype-preserving (no float promotion) so integer inputs + stay integer through the pad. Returns ``src`` unchanged when no + padding is needed. + """ + h, w = src.shape + if (h, w) == (h2, w2): + return src + padded = np.empty((h2, w2), dtype=src.dtype) + padded[:h, :w] = src + if h2 != h: + padded[h:, :w] = src[h - 1:h, :] + if w2 != w: + padded[:h, w:] = src[:, w - 1:w] + if h2 != h and w2 != w: + padded[h:, w:] = src[h - 1, w - 1] + return padded + + +def _block_reduce_2d(arr2d, method, nodata=None): + """2x block-reduce a single 2D plane using *method*. + + Output dimensions follow GDAL's ceil semantics (5x5 -> 3x3): every + base pixel contributes to the overview. Inputs with an odd height + or width are padded with NaN (or float-promoted and NaN-padded for + integer dtypes) along the trailing edge so the 2x2 block reshape + still works; the nan-aware aggregations then exclude the padded + cells from the residual blocks. ``nearest`` strides directly, + ``cubic`` replicates the last row/col so the spline kernel sees a + valid neighbourhood at the trailing edge, and ``mode`` ignores + NaN-padded cells when counting. Even-sized inputs take the + existing fast path with no padding overhead (issue #2105). + + 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). Cubic on integer dtypes follows the same path via a + float64 promotion so NaN can carry through the spline, with a + ``np.round(...).astype(arr2d.dtype)`` at the end to keep the cast + well-defined (issue #1975). + """ + h, w = arr2d.shape + # GDAL-style ceil semantics keep the trailing row/col when the input + # is odd-sized. Floor-and-crop drops those pixels and produces an + # overview that no longer covers the source extent (issue #2105). + oh, ow = (h + 1) // 2, (w + 1) // 2 + h2, w2 = 2 * oh, 2 * ow + need_pad = (h2, w2) != (h, w) + + if method == 'nearest': + # Top-left pixel of each 2x2 block. For odd-sized inputs the + # trailing 1x2 / 2x1 / 1x1 residual block's top-left is still + # the source pixel at (2*i, 2*j), so direct striding gives the + # correct ceil-shape output without padding. The trailing + # ``.copy()`` materialises the strided view so the caller owns + # a contiguous buffer (downstream nodata-mask rewrites mutate + # the result in place); without it those writes would also + # touch ``arr2d``. + return arr2d[::2, ::2].copy() + + if method == 'cubic': + try: + from scipy.ndimage import zoom + except ImportError: + raise ImportError( + "scipy is required for cubic overview resampling. " + "Install it with: pip install scipy") + # When ``nodata`` is supplied on a float array, the writer has + # already rewritten NaN to the sentinel value upstream. Feeding + # that sentinel-poisoned array straight into ``zoom`` blends the + # sentinel into neighbouring cells and produces ringing + # artefacts near nodata borders (issue #1623, same root cause + # as #1613 but for the cubic branch). + # + # Mask the sentinel back to NaN before the spline so the + # interpolation does not treat it as signal, run cubic with + # ``prefilter=False`` so a single NaN does not poison the entire + # row/column (the default B-spline prefilter is global), then + # rewrite any NaN in the result back to the sentinel so the + # on-disk overview keeps the same convention as the + # full-resolution band. The ``prefilter=False`` switch only + # fires when a sentinel was actually found in the input, so the + # default cubic semantics still apply to inputs without nodata. + # + # Integer rasters take the same path via a float64 promotion so + # NaN can carry through the spline; the result is rewritten + # back to the sentinel and rounded before casting to the source + # integer dtype (issue #1975, integer mirror of #1623). + # + # Odd-sized inputs are first edge-replicated to even dimensions + # so the 0.5 zoom factor lands on the GDAL ceil shape and the + # spline sees a valid neighbourhood at the trailing edge. + # Replicating sentinel values is fine: they are masked to NaN + # in the steps below alongside any in-bounds sentinel, and the + # post-zoom NaN->sentinel restore stamps the trailing overview + # pixels back to the sentinel as expected (issue #2105). The + # ``_replicate_pad_2d`` helper is module-level so the cubic + # path does not redefine a closure on every call. + + if (nodata is not None + and arr2d.dtype.kind == 'f' + and not np.isnan(nodata)): + try: + sentinel = arr2d.dtype.type(nodata) + except (OverflowError, ValueError): + sentinel = None + if sentinel is not None: + mask = arr2d == sentinel + if mask.any(): + src = _replicate_pad_2d(arr2d, h2, w2) + src_mask = src == sentinel + masked = np.where(src_mask, np.float64('nan'), src) + with warnings.catch_warnings(): + warnings.simplefilter('ignore', RuntimeWarning) + result = zoom(masked, 0.5, order=3, + prefilter=False) + nan_mask = np.isnan(result) + if nan_mask.any(): + result = result.copy() + result[nan_mask] = float(nodata) + return result.astype(arr2d.dtype) + # No sentinel pixels present; fall through to the + # default cubic path below so the spline runs without + # the prefilter=False NaN-safety branch. + nodata_int = _resolve_int_nodata(arr2d.dtype, nodata) + if nodata_int is not None: + sentinel = arr2d.dtype.type(nodata_int) + mask = arr2d == sentinel + if mask.any(): + src = _replicate_pad_2d(arr2d, h2, w2) + src_mask = src == sentinel + masked = np.where(src_mask, np.float64('nan'), + src.astype(np.float64)) + with warnings.catch_warnings(): + warnings.simplefilter('ignore', RuntimeWarning) + result = zoom(masked, 0.5, order=3, + prefilter=False) + nan_mask = np.isnan(result) + if nan_mask.any(): + result = np.where(nan_mask, float(nodata_int), + result) + return np.round(result).astype(arr2d.dtype) + # No sentinel pixels present in this integer array; fall + # through to the default cubic path below. + src = _replicate_pad_2d(arr2d, h2, w2) + return zoom(src, 0.5, order=3).astype(arr2d.dtype) + + if method == 'mode': + # Most-common value per 2x2 block (useful for classified rasters). + # Vectorized: sort each 4-cell block, then for each position count + # how many cells equal it. argmax picks the leftmost max-count + # position, which (post-sort) is the smallest tied value, matching + # the prior np.unique+argmax tie-break behavior ("lowest wins"). + # + # Odd-sized inputs are float-promoted and NaN-padded so the + # residual block has at least one valid cell plus NaN cells that + # never match themselves (NaN != NaN). Counts for the NaN cells + # stay at zero, so argmax picks a valid cell, and numpy's sort + # places NaN at the end so the "lowest valid wins" tie-break + # still holds (issue #2105). The result is cast back to the + # source dtype. + if need_pad: + padded = np.full((h2, w2), np.nan, dtype=np.float64) + padded[:h, :w] = arr2d + blocks = padded.reshape(oh, 2, ow, 2).transpose(0, 2, 1, 3).reshape(oh, ow, 4) + srt = np.sort(blocks, axis=-1) + counts = np.zeros(srt.shape, dtype=np.int8) + for i in range(4): + counts[..., i] = np.sum(srt == srt[..., i:i + 1], axis=-1) + pick = np.argmax(counts, axis=-1) + picked = np.take_along_axis(srt, pick[..., None], axis=-1).squeeze(-1) + if arr2d.dtype.kind in ('i', 'u'): + return np.round(picked).astype(arr2d.dtype) + return picked.astype(arr2d.dtype) + blocks = arr2d.reshape(oh, 2, ow, 2).transpose(0, 2, 1, 3).reshape(oh, ow, 4) + srt = np.sort(blocks, axis=-1) + counts = np.empty_like(srt, dtype=np.int8) + for i in range(4): + counts[..., i] = np.sum(srt == srt[..., i:i + 1], axis=-1) + pick = np.argmax(counts, axis=-1) + return np.take_along_axis(srt, pick[..., None], axis=-1).squeeze(-1) + + # Block reshape for mean/min/max/median. Odd-sized inputs are padded + # with NaN along the trailing edge so the reshape works on every + # input size; the padded cells are naturally excluded by the + # nan-aware aggregations below (issue #2105). The padded buffer + # keeps the source's float dtype (float32 inputs stay float32) so + # an odd-shape read does not pay a 2x intermediate-memory cost for + # a silent promote to float64. + if arr2d.dtype.kind == 'f': + if need_pad: + padded = np.full((h2, w2), np.nan, dtype=arr2d.dtype) + padded[:h, :w] = arr2d + blocks = padded.reshape(oh, 2, ow, 2) + else: + blocks = arr2d.reshape(oh, 2, ow, 2) + # When a sentinel was used in place of NaN by an upstream + # NaN-to-sentinel rewrite, mask it back to NaN here so nanmean / + # nanmin / nanmax / nanmedian honour the missing-data semantic. + # Without this the sentinel value participates in the reduction + # and poisons the overview (issue #1613). Match the upstream + # NaN->sentinel rewrite gate (``not np.isnan(nodata)``) so that + # ``nodata=+/-inf`` is masked here too. + if nodata is not None and not np.isnan(nodata): + try: + sentinel = arr2d.dtype.type(nodata) + except (OverflowError, ValueError): + sentinel = None + if sentinel is not None: + mask = blocks == sentinel + if mask.any(): + # ``np.where(mask, nan, blocks)`` produces a fresh + # array so the caller's input is not mutated. + blocks = np.where(mask, np.float64('nan'), blocks) + else: + if need_pad: + blocks = np.full((h2, w2), np.nan, dtype=np.float64) + blocks[:h, :w] = arr2d + blocks = blocks.reshape(oh, 2, ow, 2) + else: + blocks = arr2d.astype(np.float64).reshape(oh, 2, ow, 2) + # Integer rasters with a sentinel still need the NaN-mask the + # float branch above applies. The sentinel-equal cells in the + # original integer array are computed at the dtype's native + # width below (so a 64-bit sentinel near INT64_MAX still + # matches, which a float-cast comparison would miss) and the + # mask is then broadcast into the float64 ``blocks`` view. + # Without this, nanmean / nanmin / nanmax / nanmedian average + # the sentinel value into surrounding valid cells and produce + # overview pixels that the read side cannot remap because they + # no longer equal the sentinel. 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. + nodata_int = _resolve_int_nodata(arr2d.dtype, nodata) + if nodata_int is not None: + sentinel = arr2d.dtype.type(nodata_int) + if need_pad: + # Compute the sentinel mask at the integer's native width + # before padding, so a 64-bit sentinel near INT64_MAX (which + # is not exactly representable in float64) still matches. + # The mask is then padded with False to the same (h2, w2) + # shape as the float-promoted blocks view, so the residual + # padded cells (already NaN in ``blocks``) are not also + # rewritten via the where below. + int_mask = np.zeros((h2, w2), dtype=bool) + int_mask[:h, :w] = arr2d == sentinel + mask = int_mask.reshape(oh, 2, ow, 2) + else: + # 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 = arr2d.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 + # output is the desired signal that the caller rewrites to the + # sentinel, so suppress the warning locally to keep COG writes quiet. + with warnings.catch_warnings(): + warnings.simplefilter('ignore', RuntimeWarning) + if method == 'mean': + result = np.nanmean(blocks, axis=(1, 3)) + elif method == 'min': + result = np.nanmin(blocks, axis=(1, 3)) + elif method == 'max': + result = np.nanmax(blocks, axis=(1, 3)) + elif method == 'median': + flat = blocks.transpose(0, 2, 1, 3).reshape(oh, ow, 4) + result = np.nanmedian(flat, axis=2) + else: + raise ValueError( + f"Unknown overview resampling method: {method!r}. " + f"Use one of: {OVERVIEW_METHODS}") + + if arr2d.dtype.kind != 'f': + # 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. + nan_mask = np.isnan(result) + if nan_mask.any(): + nodata_int = _resolve_int_nodata(arr2d.dtype, nodata) + if nodata_int is not None: + result = np.where(nan_mask, float(nodata_int), result) + return np.round(result).astype(arr2d.dtype) + return result.astype(arr2d.dtype) + + +def _make_overview(arr: np.ndarray, method: str = 'mean', + nodata=None) -> np.ndarray: + """Generate a 2x decimated overview. + + Parameters + ---------- + arr : np.ndarray + 2D or 3D (height, width, bands) array. + method : str + Resampling method: 'mean' (default), 'nearest', 'min', 'max', + 'median', 'mode', or 'cubic'. + nodata : scalar or None + 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 + ------- + np.ndarray + Half-resolution array. + """ + if arr.ndim == 3: + bands = [_block_reduce_2d(arr[:, :, b], method, nodata=nodata) + for b in range(arr.shape[2])] + return np.stack(bands, axis=2) + return _block_reduce_2d(arr, method, nodata=nodata) diff --git a/xrspatial/geotiff/_writer.py b/xrspatial/geotiff/_writer.py index 0cbc11b6f..694da76c2 100644 --- a/xrspatial/geotiff/_writer.py +++ b/xrspatial/geotiff/_writer.py @@ -99,6 +99,22 @@ _serialize_tag_value, _should_use_bigtiff_streaming, ) +# Overview pyramid helpers (``_make_overview``, ``_block_reduce_2d``, +# ``_replicate_pad_2d``, ``_resolve_int_nodata``, +# ``_validate_overview_levels``) and the ``OVERVIEW_METHODS`` / +# ``_MAX_OVERVIEW_LEVELS`` constants live in ``_overview.py``. Re-export +# them here so internal call sites and external importers (the +# ``_writers`` subpackage, ``_gpu_decode``, tests) keep using the +# ``xrspatial.geotiff._writer`` import path. See issue #2259. +from ._overview import ( # noqa: F401 + OVERVIEW_METHODS, + _MAX_OVERVIEW_LEVELS, + _block_reduce_2d, + _make_overview, + _replicate_pad_2d, + _resolve_int_nodata, + _validate_overview_levels, +) # Tag IDs the writer must never accept from ``extra_tags``. NewSubfileType # (254) is a per-IFD status flag the writer emits on its own for overview @@ -349,14 +365,6 @@ def _compression_tag(compression_name: str) -> int: return _map[name] -OVERVIEW_METHODS = ('mean', 'nearest', 'min', 'max', 'median', 'mode', 'cubic') - -#: Maximum number of overview levels generated by auto-overview mode in COG -#: writes. 8 halvings = 1/256 of the original resolution, which is enough -#: for any practical raster. Pass ``overview_levels=[...]`` explicitly to -#: override. -_MAX_OVERVIEW_LEVELS = 8 - #: Total uncompressed payload (bytes) below which the strip and tile #: writers stay sequential. The thread-pool startup cost dominates on #: small rasters; above this size the per-block compression cost more @@ -366,452 +374,6 @@ def _compression_tag(compression_name: str) -> int: _PARALLEL_MIN_BYTES = 4 * 1024 * 1024 -def _validate_overview_levels(overview_levels, height=None, width=None): - """Validate and normalise an explicit ``overview_levels`` list. - - Each entry is a decimation factor relative to full resolution. - Factors must be strictly increasing integers >= 2, and each must - be a power of two because the underlying block reducer only does - 2x decimation per step (issue #1766 — prior to that fix, the values - were ignored and only the list length mattered). - - When ``height`` and ``width`` are supplied, each factor is also - checked against the input shape: a factor F is feasible only if - ``height // F >= 1`` and ``width // F >= 1``. Asking for a factor - that would reduce the source below 1 pixel in either dimension - raises ``ValueError`` instead of silently writing a zero-sized - overview IFD. - - Returns the validated list of ints. ``None`` passes through so the - caller can run its auto-generation path. - """ - if overview_levels is None: - return None - if not isinstance(overview_levels, (list, tuple)): - raise ValueError( - f"overview_levels must be a list or tuple of ints, got " - f"{type(overview_levels).__name__}.") - if len(overview_levels) == 0: - return [] - cleaned = [] - prev = 1 - for i, level in enumerate(overview_levels): - # Reject bools explicitly; ``bool`` is an ``int`` subclass and - # ``True``/``False`` would otherwise sneak through the integer - # check below. - if isinstance(level, bool) or not isinstance(level, (int, np.integer)): - raise ValueError( - f"overview_levels[{i}] must be an int >= 2, got " - f"{level!r} (type {type(level).__name__}).") - level = int(level) - if level < 2: - raise ValueError( - f"overview_levels[{i}] must be >= 2 (1 is the original " - f"full-resolution band), got {level}.") - if level <= prev: - raise ValueError( - f"overview_levels must be strictly increasing, got " - f"{list(overview_levels)} (entry at index {i} is " - f"{level}, previous was {prev}).") - # Power-of-two check: the underlying ``_make_overview`` only - # halves, so reaching factor N takes log2(N) halvings and N - # must be a power of two for the cumulative factor to land on - # the requested value exactly. - if (level & (level - 1)) != 0: - raise ValueError( - f"overview_levels[{i}]={level} is not a power of two. " - f"Only power-of-two decimation factors are supported " - f"(2, 4, 8, 16, ...).") - # Shape feasibility: refuse factors that would shrink the - # raster below 1 pixel. ``_block_reduce_2d`` halves via - # ``(dim // 2) * 2`` which produces a zero-sized array once - # ``dim < 2``, and chaining further halvings keeps it at zero. - # Without this check the writer silently emits zero-sized - # overview IFDs. - if height is not None and width is not None: - if height // level < 1 or width // level < 1: - raise ValueError( - f"overview_levels[{i}]={level} is too large for " - f"input shape ({height}, {width}); decimation " - f"would produce a zero-sized overview.") - cleaned.append(level) - prev = level - return cleaned - - -def _resolve_int_nodata(dtype, nodata): - """Return ``int(nodata)`` if it is representable as *dtype*, else None. - - Folds the three checks that gate the integer sentinel-to-NaN mask in - :func:`_block_reduce_2d` into one call: ``dtype`` is integer, the - sentinel is finite and integer-valued, and the integer fits the - dtype range. Out-of-range pairs like ``uint16`` + ``GDAL_NODATA=-9999`` - return None so the caller stays a no-op rather than tripping - ``OverflowError`` on the dtype cast. Mirrors - ``_int_nodata_in_range`` in ``_reader.py``. - """ - if nodata is None or dtype.kind not in ('i', 'u'): - return None - if not np.isfinite(nodata) or not float(nodata).is_integer(): - return None - nodata_int = int(nodata) - info = np.iinfo(dtype) - if info.min <= nodata_int <= info.max: - return nodata_int - return None - - -def _replicate_pad_2d(src, h2, w2): - """Pad ``src`` to ``(h2, w2)`` by replicating its trailing row/col. - - Used by the cubic overview branch to extend an odd-sized input to - even dimensions so a 0.5-zoom lands on the GDAL ceil shape and the - spline kernel sees a valid neighbourhood at the trailing edge. The - helper is dtype-preserving (no float promotion) so integer inputs - stay integer through the pad. Returns ``src`` unchanged when no - padding is needed. - """ - h, w = src.shape - if (h, w) == (h2, w2): - return src - padded = np.empty((h2, w2), dtype=src.dtype) - padded[:h, :w] = src - if h2 != h: - padded[h:, :w] = src[h - 1:h, :] - if w2 != w: - padded[:h, w:] = src[:, w - 1:w] - if h2 != h and w2 != w: - padded[h:, w:] = src[h - 1, w - 1] - return padded - - -def _block_reduce_2d(arr2d, method, nodata=None): - """2x block-reduce a single 2D plane using *method*. - - Output dimensions follow GDAL's ceil semantics (5x5 -> 3x3): every - base pixel contributes to the overview. Inputs with an odd height - or width are padded with NaN (or float-promoted and NaN-padded for - integer dtypes) along the trailing edge so the 2x2 block reshape - still works; the nan-aware aggregations then exclude the padded - cells from the residual blocks. ``nearest`` strides directly, - ``cubic`` replicates the last row/col so the spline kernel sees a - valid neighbourhood at the trailing edge, and ``mode`` ignores - NaN-padded cells when counting. Even-sized inputs take the - existing fast path with no padding overhead (issue #2105). - - 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). Cubic on integer dtypes follows the same path via a - float64 promotion so NaN can carry through the spline, with a - ``np.round(...).astype(arr2d.dtype)`` at the end to keep the cast - well-defined (issue #1975). - """ - h, w = arr2d.shape - # GDAL-style ceil semantics keep the trailing row/col when the input - # is odd-sized. Floor-and-crop drops those pixels and produces an - # overview that no longer covers the source extent (issue #2105). - oh, ow = (h + 1) // 2, (w + 1) // 2 - h2, w2 = 2 * oh, 2 * ow - need_pad = (h2, w2) != (h, w) - - if method == 'nearest': - # Top-left pixel of each 2x2 block. For odd-sized inputs the - # trailing 1x2 / 2x1 / 1x1 residual block's top-left is still - # the source pixel at (2*i, 2*j), so direct striding gives the - # correct ceil-shape output without padding. The trailing - # ``.copy()`` materialises the strided view so the caller owns - # a contiguous buffer (downstream nodata-mask rewrites mutate - # the result in place); without it those writes would also - # touch ``arr2d``. - return arr2d[::2, ::2].copy() - - if method == 'cubic': - try: - from scipy.ndimage import zoom - except ImportError: - raise ImportError( - "scipy is required for cubic overview resampling. " - "Install it with: pip install scipy") - # When ``nodata`` is supplied on a float array, the writer has - # already rewritten NaN to the sentinel value upstream. Feeding - # that sentinel-poisoned array straight into ``zoom`` blends the - # sentinel into neighbouring cells and produces ringing - # artefacts near nodata borders (issue #1623, same root cause - # as #1613 but for the cubic branch). - # - # Mask the sentinel back to NaN before the spline so the - # interpolation does not treat it as signal, run cubic with - # ``prefilter=False`` so a single NaN does not poison the entire - # row/column (the default B-spline prefilter is global), then - # rewrite any NaN in the result back to the sentinel so the - # on-disk overview keeps the same convention as the - # full-resolution band. The ``prefilter=False`` switch only - # fires when a sentinel was actually found in the input, so the - # default cubic semantics still apply to inputs without nodata. - # - # Integer rasters take the same path via a float64 promotion so - # NaN can carry through the spline; the result is rewritten - # back to the sentinel and rounded before casting to the source - # integer dtype (issue #1975, integer mirror of #1623). - # - # Odd-sized inputs are first edge-replicated to even dimensions - # so the 0.5 zoom factor lands on the GDAL ceil shape and the - # spline sees a valid neighbourhood at the trailing edge. - # Replicating sentinel values is fine: they are masked to NaN - # in the steps below alongside any in-bounds sentinel, and the - # post-zoom NaN->sentinel restore stamps the trailing overview - # pixels back to the sentinel as expected (issue #2105). The - # ``_replicate_pad_2d`` helper is module-level so the cubic - # path does not redefine a closure on every call. - - if (nodata is not None - and arr2d.dtype.kind == 'f' - and not np.isnan(nodata)): - try: - sentinel = arr2d.dtype.type(nodata) - except (OverflowError, ValueError): - sentinel = None - if sentinel is not None: - mask = arr2d == sentinel - if mask.any(): - src = _replicate_pad_2d(arr2d, h2, w2) - src_mask = src == sentinel - masked = np.where(src_mask, np.float64('nan'), src) - with warnings.catch_warnings(): - warnings.simplefilter('ignore', RuntimeWarning) - result = zoom(masked, 0.5, order=3, - prefilter=False) - nan_mask = np.isnan(result) - if nan_mask.any(): - result = result.copy() - result[nan_mask] = float(nodata) - return result.astype(arr2d.dtype) - # No sentinel pixels present; fall through to the - # default cubic path below so the spline runs without - # the prefilter=False NaN-safety branch. - nodata_int = _resolve_int_nodata(arr2d.dtype, nodata) - if nodata_int is not None: - sentinel = arr2d.dtype.type(nodata_int) - mask = arr2d == sentinel - if mask.any(): - src = _replicate_pad_2d(arr2d, h2, w2) - src_mask = src == sentinel - masked = np.where(src_mask, np.float64('nan'), - src.astype(np.float64)) - with warnings.catch_warnings(): - warnings.simplefilter('ignore', RuntimeWarning) - result = zoom(masked, 0.5, order=3, - prefilter=False) - nan_mask = np.isnan(result) - if nan_mask.any(): - result = np.where(nan_mask, float(nodata_int), - result) - return np.round(result).astype(arr2d.dtype) - # No sentinel pixels present in this integer array; fall - # through to the default cubic path below. - src = _replicate_pad_2d(arr2d, h2, w2) - return zoom(src, 0.5, order=3).astype(arr2d.dtype) - - if method == 'mode': - # Most-common value per 2x2 block (useful for classified rasters). - # Vectorized: sort each 4-cell block, then for each position count - # how many cells equal it. argmax picks the leftmost max-count - # position, which (post-sort) is the smallest tied value, matching - # the prior np.unique+argmax tie-break behavior ("lowest wins"). - # - # Odd-sized inputs are float-promoted and NaN-padded so the - # residual block has at least one valid cell plus NaN cells that - # never match themselves (NaN != NaN). Counts for the NaN cells - # stay at zero, so argmax picks a valid cell, and numpy's sort - # places NaN at the end so the "lowest valid wins" tie-break - # still holds (issue #2105). The result is cast back to the - # source dtype. - if need_pad: - padded = np.full((h2, w2), np.nan, dtype=np.float64) - padded[:h, :w] = arr2d - blocks = padded.reshape(oh, 2, ow, 2).transpose(0, 2, 1, 3).reshape(oh, ow, 4) - srt = np.sort(blocks, axis=-1) - counts = np.zeros(srt.shape, dtype=np.int8) - for i in range(4): - counts[..., i] = np.sum(srt == srt[..., i:i + 1], axis=-1) - pick = np.argmax(counts, axis=-1) - picked = np.take_along_axis(srt, pick[..., None], axis=-1).squeeze(-1) - if arr2d.dtype.kind in ('i', 'u'): - return np.round(picked).astype(arr2d.dtype) - return picked.astype(arr2d.dtype) - blocks = arr2d.reshape(oh, 2, ow, 2).transpose(0, 2, 1, 3).reshape(oh, ow, 4) - srt = np.sort(blocks, axis=-1) - counts = np.empty_like(srt, dtype=np.int8) - for i in range(4): - counts[..., i] = np.sum(srt == srt[..., i:i + 1], axis=-1) - pick = np.argmax(counts, axis=-1) - return np.take_along_axis(srt, pick[..., None], axis=-1).squeeze(-1) - - # Block reshape for mean/min/max/median. Odd-sized inputs are padded - # with NaN along the trailing edge so the reshape works on every - # input size; the padded cells are naturally excluded by the - # nan-aware aggregations below (issue #2105). The padded buffer - # keeps the source's float dtype (float32 inputs stay float32) so - # an odd-shape read does not pay a 2x intermediate-memory cost for - # a silent promote to float64. - if arr2d.dtype.kind == 'f': - if need_pad: - padded = np.full((h2, w2), np.nan, dtype=arr2d.dtype) - padded[:h, :w] = arr2d - blocks = padded.reshape(oh, 2, ow, 2) - else: - blocks = arr2d.reshape(oh, 2, ow, 2) - # When a sentinel was used in place of NaN by an upstream - # NaN-to-sentinel rewrite, mask it back to NaN here so nanmean / - # nanmin / nanmax / nanmedian honour the missing-data semantic. - # Without this the sentinel value participates in the reduction - # and poisons the overview (issue #1613). Match the upstream - # NaN->sentinel rewrite gate (``not np.isnan(nodata)``) so that - # ``nodata=+/-inf`` is masked here too. - if nodata is not None and not np.isnan(nodata): - try: - sentinel = arr2d.dtype.type(nodata) - except (OverflowError, ValueError): - sentinel = None - if sentinel is not None: - mask = blocks == sentinel - if mask.any(): - # ``np.where(mask, nan, blocks)`` produces a fresh - # array so the caller's input is not mutated. - blocks = np.where(mask, np.float64('nan'), blocks) - else: - if need_pad: - blocks = np.full((h2, w2), np.nan, dtype=np.float64) - blocks[:h, :w] = arr2d - blocks = blocks.reshape(oh, 2, ow, 2) - else: - blocks = arr2d.astype(np.float64).reshape(oh, 2, ow, 2) - # Integer rasters with a sentinel still need the NaN-mask the - # float branch above applies. The sentinel-equal cells in the - # original integer array are computed at the dtype's native - # width below (so a 64-bit sentinel near INT64_MAX still - # matches, which a float-cast comparison would miss) and the - # mask is then broadcast into the float64 ``blocks`` view. - # Without this, nanmean / nanmin / nanmax / nanmedian average - # the sentinel value into surrounding valid cells and produce - # overview pixels that the read side cannot remap because they - # no longer equal the sentinel. 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. - nodata_int = _resolve_int_nodata(arr2d.dtype, nodata) - if nodata_int is not None: - sentinel = arr2d.dtype.type(nodata_int) - if need_pad: - # Compute the sentinel mask at the integer's native width - # before padding, so a 64-bit sentinel near INT64_MAX (which - # is not exactly representable in float64) still matches. - # The mask is then padded with False to the same (h2, w2) - # shape as the float-promoted blocks view, so the residual - # padded cells (already NaN in ``blocks``) are not also - # rewritten via the where below. - int_mask = np.zeros((h2, w2), dtype=bool) - int_mask[:h, :w] = arr2d == sentinel - mask = int_mask.reshape(oh, 2, ow, 2) - else: - # 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 = arr2d.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 - # output is the desired signal that the caller rewrites to the - # sentinel, so suppress the warning locally to keep COG writes quiet. - with warnings.catch_warnings(): - warnings.simplefilter('ignore', RuntimeWarning) - if method == 'mean': - result = np.nanmean(blocks, axis=(1, 3)) - elif method == 'min': - result = np.nanmin(blocks, axis=(1, 3)) - elif method == 'max': - result = np.nanmax(blocks, axis=(1, 3)) - elif method == 'median': - flat = blocks.transpose(0, 2, 1, 3).reshape(oh, ow, 4) - result = np.nanmedian(flat, axis=2) - else: - raise ValueError( - f"Unknown overview resampling method: {method!r}. " - f"Use one of: {OVERVIEW_METHODS}") - - if arr2d.dtype.kind != 'f': - # 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. - nan_mask = np.isnan(result) - if nan_mask.any(): - nodata_int = _resolve_int_nodata(arr2d.dtype, nodata) - if nodata_int is not None: - result = np.where(nan_mask, float(nodata_int), result) - return np.round(result).astype(arr2d.dtype) - return result.astype(arr2d.dtype) - - -def _make_overview(arr: np.ndarray, method: str = 'mean', - nodata=None) -> np.ndarray: - """Generate a 2x decimated overview. - - Parameters - ---------- - arr : np.ndarray - 2D or 3D (height, width, bands) array. - method : str - Resampling method: 'mean' (default), 'nearest', 'min', 'max', - 'median', 'mode', or 'cubic'. - nodata : scalar or None - 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 - ------- - np.ndarray - Half-resolution array. - """ - if arr.ndim == 3: - bands = [_block_reduce_2d(arr[:, :, b], method, nodata=nodata) - for b in range(arr.shape[2])] - return np.stack(bands, axis=2) - return _block_reduce_2d(arr, method, nodata=nodata) - - # --------------------------------------------------------------------------- # Tag serialization helpers (``_float_to_rational``, ``_serialize_tag_value``, # ``_pack_tag_value``, ``_build_ifd``) live in ``_write_layout.py`` and are diff --git a/xrspatial/geotiff/_writers/gpu.py b/xrspatial/geotiff/_writers/gpu.py index b8d944b99..cb6cff32c 100644 --- a/xrspatial/geotiff/_writers/gpu.py +++ b/xrspatial/geotiff/_writers/gpu.py @@ -609,7 +609,7 @@ def _gpu_compress_to_part(gpu_arr, w, h, spp): # Overview generation -- mirrors the CPU writer's 8-level cap. if cog: if overview_levels is None: - from .._writer import _MAX_OVERVIEW_LEVELS + from .._overview import _MAX_OVERVIEW_LEVELS # Auto-generated lists hold actual decimation factors (2, # 4, 8, ...) so the loop below treats auto-generated and # user-supplied lists identically (issue #1766). @@ -628,7 +628,7 @@ def _gpu_compress_to_part(gpu_arr, w, h, spp): # strictly increasing, feasible for the input shape. # Previously the values were ignored and only the list # length mattered (issue #1766). - from .._writer import _validate_overview_levels + from .._overview import _validate_overview_levels overview_levels = _validate_overview_levels( overview_levels, height=height, width=width)