From fe86d7335f4a8d0167daa1344abaf59ce16a5e8c Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 11 May 2026 11:40:56 -0700 Subject: [PATCH 1/3] Fix open_geotiff(gpu=True) silently dropping window= and band= (#1605) The GPU dispatch path in open_geotiff forwarded dtype/overview_level/ name/chunks/max_pixels but not window or band, and read_geotiff_gpu did not declare either kwarg. Callers asking for a windowed or single-band GPU read got the full multi-band file back without warning -- a silent data-correctness drop that broke backend substitutability against the CPU eager path and read_geotiff_dask. Add window and band kwargs to read_geotiff_gpu and forward them through open_geotiff. The fix slices on device after the full-image GPU decode completes, which preserves correctness today; tile-level short-circuit is a future optimisation. Coords and attrs['transform'] match the eager path bit-for-bit so the only observable difference is the array backend. Regression tests cover direct read_geotiff_gpu calls, dispatch through open_geotiff(gpu=True), single-band-and-windowed combinations, and the out-of-bounds error paths. --- xrspatial/geotiff/__init__.py | 136 +++++++++++- .../tests/test_gpu_window_band_1605.py | 196 ++++++++++++++++++ 2 files changed, 324 insertions(+), 8 deletions(-) create mode 100644 xrspatial/geotiff/tests/test_gpu_window_band_1605.py diff --git a/xrspatial/geotiff/__init__.py b/xrspatial/geotiff/__init__.py index e66e6dc9..44bad5c8 100644 --- a/xrspatial/geotiff/__init__.py +++ b/xrspatial/geotiff/__init__.py @@ -498,6 +498,7 @@ def open_geotiff(source, *, dtype=None, window=None, if gpu: return read_geotiff_gpu(source, dtype=dtype, overview_level=overview_level, + window=window, band=band, name=name, chunks=chunks, max_pixels=max_pixels) @@ -1933,9 +1934,70 @@ def _apply_orientation_geo_info(geo_info, orientation: int, return geo_info +def _gpu_apply_window_band(arr_gpu, geo_info, *, window, band): + """Slice a fully-decoded GPU array down to a window and/or band. + + Used by ``read_geotiff_gpu`` to keep the public surface in line with + ``open_geotiff`` and ``read_geotiff_dask``: callers can pass ``window`` + and ``band``, and the returned DataArray covers exactly that subset. + + The current implementation slices on device after the full-image GPU + decode is complete. That preserves correctness but does no I/O + savings -- a future PR can short-circuit tile decode for partial + windows. For ``band`` selection, the savings are also post-decode + because the planar=1 (chunky) tile assembly returns all bands in a + single GPU buffer. + + Returns ``(arr_gpu, coords)`` where ``coords`` is a dict with + ``y`` / ``x`` numpy arrays sized to the output array. The caller is + responsible for setting ``attrs['transform']`` to the windowed origin + via ``_populate_attrs_from_geo_info(..., window=window)`` so the array + and the transform agree. + """ + if window is not None: + r0, c0, r1, c1 = window + arr_gpu = arr_gpu[r0:r1, c0:c1] + out_h = r1 - r0 + out_w = c1 - c0 + # Mirror the eager-numpy windowed coord computation in + # open_geotiff so the GPU-windowed coords carry the same + # absolute pixel-center values as the CPU path. + t = geo_info.transform + if t is None: + coords = { + 'y': np.arange(out_h, dtype=np.int64), + 'x': np.arange(out_w, dtype=np.int64), + } + else: + if geo_info.raster_type == RASTER_PIXEL_IS_POINT: + full_x = (np.arange(c0, c1, dtype=np.float64) + * t.pixel_width + t.origin_x) + full_y = (np.arange(r0, r1, dtype=np.float64) + * t.pixel_height + t.origin_y) + else: + full_x = (np.arange(c0, c1, dtype=np.float64) + * t.pixel_width + t.origin_x + + t.pixel_width * 0.5) + full_y = (np.arange(r0, r1, dtype=np.float64) + * t.pixel_height + t.origin_y + + t.pixel_height * 0.5) + coords = {'y': full_y, 'x': full_x} + else: + out_h = arr_gpu.shape[0] + out_w = arr_gpu.shape[1] + coords = _geo_to_coords(geo_info, out_h, out_w) + + if band is not None and arr_gpu.ndim == 3: + arr_gpu = arr_gpu[:, :, band] + + return arr_gpu, coords + + def read_geotiff_gpu(source: str, *, dtype=None, overview_level: int | None = None, + window: tuple | None = None, + band: int | None = None, name: str | None = None, chunks: int | tuple | None = None, max_pixels: int | None = None, @@ -1958,6 +2020,18 @@ def read_geotiff_gpu(source: str, *, File path. overview_level : int or None Overview level (0 = full resolution). + window : tuple or None + ``(row_start, col_start, row_stop, col_stop)`` for windowed + reading. None reads the full raster. The GPU pipeline currently + decodes all tiles and slices on device after assembly, so the + kwarg restores API parity with ``open_geotiff`` and + ``read_geotiff_dask`` but does not yet skip I/O for partial + windows. The returned coords, ``attrs['transform']``, and + shape match the eager numpy path. + band : int or None + Zero-based band index. None returns all bands (3D output for + multi-band files, 2D for single-band). Selecting a single band + yields a 2D DataArray. chunks : int, tuple, or None If set, return a Dask-chunked CuPy DataArray. int for square chunks, (row, col) tuple for rectangular. @@ -2051,6 +2125,19 @@ def read_geotiff_gpu(source: str, *, if max_pixels is None: max_pixels = MAX_PIXELS_DEFAULT + # Window basic shape check happens here; full bounds-vs-file validation + # runs after the IFD parse below so we can compare against the chosen + # overview level's actual height/width. ``band`` is similarly validated + # against ``ifd.samples_per_pixel`` after the header parse. + if window is not None: + if len(window) != 4: + raise ValueError( + f"window must be a 4-tuple (r0, c0, r1, c1), got {window!r}") + w_r0, w_c0, w_r1, w_c1 = window + if w_r0 >= w_r1 or w_c0 >= w_c1 or w_r0 < 0 or w_c0 < 0: + raise ValueError( + f"window={window} has non-positive size or negative origin.") + # Parse metadata on CPU (fast, <1ms) src = _FileSource(source) data = src.read_all() @@ -2075,6 +2162,29 @@ def read_geotiff_gpu(source: str, *, # only need to copy the post-flip geo_info back here. orientation = ifd.orientation + # Validate band against the selected IFD's sample count. + # ``samples_per_pixel`` is at least 1 for any valid TIFF; we treat + # ``band=0`` as "first band" for single-band files too so the + # behaviour mirrors ``read_geotiff_dask``. + ifd_samples = ifd.samples_per_pixel + if band is not None: + if ifd_samples <= 1: + if band != 0: + raise IndexError( + f"band={band} requested on a single-band file.") + elif not 0 <= band < ifd_samples: + raise IndexError( + f"band={band} out of range for {ifd_samples}-band file.") + + # Validate window upper bounds against the selected IFD's extent. + if window is not None: + w_r0, w_c0, w_r1, w_c1 = window + ifd_h, ifd_w = ifd.height, ifd.width + if w_r1 > ifd_h or w_c1 > ifd_w: + raise ValueError( + f"window={window} is outside the source extent " + f"({ifd_h}x{ifd_w}).") + if not ifd.is_tiled: # Fall back to CPU for stripped files. read_to_array remaps # the array but only updates geo_info.transform for orientations @@ -2089,12 +2199,11 @@ def read_geotiff_gpu(source: str, *, geo_info = _apply_orientation_geo_info( geo_info, orientation, file_h=ifd.height, file_w=ifd.width) - coords = _geo_to_coords(geo_info, arr_gpu.shape[0], arr_gpu.shape[1]) if name is None: import os name = os.path.splitext(os.path.basename(source))[0] attrs = {} - _populate_attrs_from_geo_info(attrs, geo_info) + _populate_attrs_from_geo_info(attrs, geo_info, window=window) # Apply nodata mask + record sentinel so the GPU read agrees # with the CPU eager path. Without this, integer rasters keep # the literal sentinel value and float rasters keep the @@ -2107,6 +2216,12 @@ def read_geotiff_gpu(source: str, *, target = np.dtype(dtype) _validate_dtype_cast(np.dtype(str(arr_gpu.dtype)), target) arr_gpu = arr_gpu.astype(target) + # Apply window/band slicing post-decode. The stripped CPU + # fallback already produces the full-image array; slice on the + # GPU so the result matches ``open_geotiff`` / + # ``read_geotiff_dask`` semantics. + arr_gpu, coords = _gpu_apply_window_band( + arr_gpu, geo_info, window=window, band=band) # Multi-band stripped reads come back as (y, x, band); mirror # the tiled branch so dims line up with ndim. Single-band stays # 2-D ('y', 'x'). @@ -2352,16 +2467,21 @@ def _read_once(): import os name = os.path.splitext(os.path.basename(source))[0] - # Use the post-orientation array shape so coords match the array. - out_h = arr_gpu.shape[0] - out_w = arr_gpu.shape[1] - coords = _geo_to_coords(geo_info, out_h, out_w) - attrs = {} - _populate_attrs_from_geo_info(attrs, geo_info) + _populate_attrs_from_geo_info(attrs, geo_info, window=window) if nodata is not None: attrs['nodata'] = nodata + # Apply window/band slicing post-decode. Coords are derived from the + # sliced array so the (y, x) labels line up with the user's requested + # subrectangle. This mirrors the ``open_geotiff`` / ``read_geotiff_dask`` + # contract: ``attrs['transform']`` always carries the full-source + # GeoTransform shifted to the window origin (via + # ``_populate_attrs_from_geo_info(..., window=window)``), while + # ``coords['y']`` / ``coords['x']`` cover only the windowed cells. + arr_gpu, coords = _gpu_apply_window_band( + arr_gpu, geo_info, window=window, band=band) + if arr_gpu.ndim == 3: dims = ['y', 'x', 'band'] coords['band'] = np.arange(arr_gpu.shape[2]) diff --git a/xrspatial/geotiff/tests/test_gpu_window_band_1605.py b/xrspatial/geotiff/tests/test_gpu_window_band_1605.py new file mode 100644 index 00000000..1837782d --- /dev/null +++ b/xrspatial/geotiff/tests/test_gpu_window_band_1605.py @@ -0,0 +1,196 @@ +"""Regression tests for issue #1605. + +``open_geotiff(source, gpu=True, window=..., band=...)`` used to silently +drop ``window`` and ``band``: the dispatcher did not forward them, and +``read_geotiff_gpu`` did not declare them. The fix adds both kwargs to +``read_geotiff_gpu`` and threads them through the dispatcher so the GPU +path returns the same subrectangle / band selection as the CPU eager +path and the dask path. + +These tests pin: + +1. Direct ``read_geotiff_gpu(..., window=..., band=...)`` accepts the + kwargs and produces the expected shape/values. +2. ``open_geotiff(..., gpu=True, window=..., band=...)`` no longer + silently returns the full file. +3. The GPU windowed coords / transform attr match the CPU eager path + bit-for-bit (so the only difference is the array backend). +4. Bounds validation for ``window`` and ``band`` raises the same errors + the CPU paths raise. +""" +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", +) + + +@pytest.fixture +def single_band_tiff(tmp_path): + """16x20 single-band tiled tiff with a non-trivial transform.""" + from xrspatial.geotiff import to_geotiff + + arr = np.arange(16 * 20, dtype=np.float32).reshape(16, 20) + da = xr.DataArray( + arr, + dims=['y', 'x'], + coords={ + 'y': np.arange(16, dtype=np.float64) * 0.25 + 100.0, + 'x': np.arange(20, dtype=np.float64) * 0.25 + 200.0, + }, + attrs={'crs': 4326}, + ) + p = tmp_path / 'window_band_1605_single.tif' + to_geotiff(da, str(p), tile_size=8) + return str(p), arr + + +@pytest.fixture +def multi_band_tiff(tmp_path): + """16x20x3 tiled tiff -- exercises the band-selection branch.""" + from xrspatial.geotiff import to_geotiff + + arr = np.arange(16 * 20 * 3, dtype=np.float32).reshape(16, 20, 3) + da = xr.DataArray( + arr, + dims=['y', 'x', 'band'], + coords={ + 'y': np.arange(16, dtype=np.float64) * 0.25 + 100.0, + 'x': np.arange(20, dtype=np.float64) * 0.25 + 200.0, + 'band': np.arange(3), + }, + attrs={'crs': 4326}, + ) + p = tmp_path / 'window_band_1605_multi.tif' + to_geotiff(da, str(p), tile_size=8) + return str(p), arr + + +@_gpu_only +def test_read_geotiff_gpu_window_matches_eager(single_band_tiff): + """Direct call: GPU window slice matches CPU eager window slice.""" + path, source_arr = single_band_tiff + from xrspatial.geotiff import open_geotiff, read_geotiff_gpu + + window = (2, 4, 12, 14) + + cpu = open_geotiff(path, window=window) + gpu = read_geotiff_gpu(path, window=window) + + assert gpu.shape == cpu.shape == (10, 10) + np.testing.assert_array_equal(gpu.data.get(), cpu.data) + np.testing.assert_array_equal(gpu['y'].values, cpu['y'].values) + np.testing.assert_array_equal(gpu['x'].values, cpu['x'].values) + assert gpu.attrs.get('transform') == cpu.attrs.get('transform') + + +@_gpu_only +def test_open_geotiff_gpu_window_no_longer_silently_dropped(single_band_tiff): + """Issue #1605: open_geotiff(gpu=True, window=...) honors the window.""" + path, source_arr = single_band_tiff + from xrspatial.geotiff import open_geotiff + + window = (3, 5, 9, 13) + gpu = open_geotiff(path, gpu=True, window=window) + + # Pre-fix shape was the full (16, 20). Post-fix it must shrink. + assert gpu.shape == (6, 8) + np.testing.assert_array_equal( + gpu.data.get(), + source_arr[3:9, 5:13], + ) + + +@_gpu_only +def test_read_geotiff_gpu_band_selection(multi_band_tiff): + """Direct call: band=k returns the kth band as a 2D DataArray.""" + path, source_arr = multi_band_tiff + from xrspatial.geotiff import open_geotiff, read_geotiff_gpu + + cpu = open_geotiff(path, band=1) + gpu = read_geotiff_gpu(path, band=1) + + assert gpu.shape == cpu.shape == (16, 20) + assert gpu.ndim == 2 + np.testing.assert_array_equal(gpu.data.get(), cpu.data) + + +@_gpu_only +def test_open_geotiff_gpu_band_no_longer_silently_dropped(multi_band_tiff): + """Issue #1605: open_geotiff(gpu=True, band=...) honors band.""" + path, source_arr = multi_band_tiff + from xrspatial.geotiff import open_geotiff + + gpu = open_geotiff(path, gpu=True, band=2) + + # Pre-fix would have returned a (16, 20, 3) DataArray (full file). + assert gpu.shape == (16, 20) + assert gpu.ndim == 2 + np.testing.assert_array_equal( + gpu.data.get(), + source_arr[:, :, 2], + ) + + +@_gpu_only +def test_read_geotiff_gpu_window_and_band(multi_band_tiff): + """window + band combine cleanly.""" + path, source_arr = multi_band_tiff + from xrspatial.geotiff import open_geotiff, read_geotiff_gpu + + window = (1, 2, 11, 17) + cpu = open_geotiff(path, window=window, band=0) + gpu = read_geotiff_gpu(path, window=window, band=0) + + assert gpu.shape == cpu.shape == (10, 15) + assert gpu.ndim == 2 + np.testing.assert_array_equal(gpu.data.get(), cpu.data) + + +@_gpu_only +def test_read_geotiff_gpu_window_bounds_validation(single_band_tiff): + """Out-of-bounds window raises ValueError, mirroring the dask path.""" + path, _ = single_band_tiff + from xrspatial.geotiff import read_geotiff_gpu + + with pytest.raises(ValueError, match="outside the source extent"): + read_geotiff_gpu(path, window=(0, 0, 100, 100)) + + with pytest.raises(ValueError, match="non-positive size"): + read_geotiff_gpu(path, window=(5, 0, 5, 10)) + + with pytest.raises(ValueError, match="must be a 4-tuple"): + read_geotiff_gpu(path, window=(0, 0, 5)) + + +@_gpu_only +def test_read_geotiff_gpu_band_bounds_validation(multi_band_tiff, + single_band_tiff): + """Out-of-range band raises IndexError.""" + multi_path, _ = multi_band_tiff + single_path, _ = single_band_tiff + from xrspatial.geotiff import read_geotiff_gpu + + with pytest.raises(IndexError, match="out of range"): + read_geotiff_gpu(multi_path, band=10) + + with pytest.raises(IndexError, match="single-band file"): + read_geotiff_gpu(single_path, band=1) From 50a8a8c74fb02dac4f7a99df0d294f3ff527b392 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 11 May 2026 11:41:20 -0700 Subject: [PATCH 2/3] Update api-consistency sweep state for geotiff (2026-05-11) --- .claude/sweep-api-consistency-state.csv | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.claude/sweep-api-consistency-state.csv b/.claude/sweep-api-consistency-state.csv index 6ec6c20b..3403bb13 100644 --- a/.claude/sweep-api-consistency-state.csv +++ b/.claude/sweep-api-consistency-state.csv @@ -1,3 +1,3 @@ module,last_inspected,issue,severity_max,categories_found,notes -geotiff,2026-05-10,1560;1561;1562,HIGH,1;3;5,"Filed gpu kwarg type drift (1560 HIGH), backend signature parity gap (1561 HIGH), writer docstring drift (1562 MEDIUM). Prior orphan-__all__ HIGH fixed in PR 1544. Cat 1 MEDIUM (write/dask_data internal asymmetry) and Cat 4 MEDIUM (chunks default) reviewed and judged defensible per private-helper / different-semantics scope." +geotiff,2026-05-11,1605,HIGH,1;5,"Filed open_geotiff(gpu=True) silently dropping window/band (HIGH, #1605) and fix PR -- adds window+band kwargs to read_geotiff_gpu and forwards them through open_geotiff GPU dispatch. Prior issues 1560/1561/1562 from 2026-05-10 audit all CLOSED. MEDIUM: read_geotiff_dask VRT defensive fallback drops window/band/max_pixels at line 1463 (acceptable since open_geotiff routes earlier and direct callers can switch to read_vrt). write_vrt wrapper uses **kwargs instead of explicit signature -- docs list the args but inspect.signature does not; cosmetic." reproject,2026-05-10,1570,HIGH,2;5,"Filed cross-module attrs['vertical_crs'] type collision (string vs EPSG int) vs xrspatial.geotiff. Fixed in PR (TBD): reproject now writes EPSG int and preserves friendly token under vertical_datum. MEDIUM kwarg-order drift (transform_precision vs chunk_size) and missing type hints vs geotiff documented but not fixed (cosmetic, kwarg-only)." From fec43eb9d1b9221dd5b57e8dc961cd977b5724a1 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 11 May 2026 12:00:49 -0700 Subject: [PATCH 3/3] Address review on #1608: window/orientation check + chunks for stripped TIFFs Two Copilot review notes on read_geotiff_gpu: 1. Orientation tag (274) != 1 combined with window= has ambiguous semantics (file pixel space vs. display pixel space). _reader.read_to_array rejects this combo with a specific error; read_geotiff_gpu previously validated window only against the pre-orientation IFD extent and then applied the slice after the remap, which could disagree with the CPU path. Raise the same ValueError -- reusing the CPU reader's wording -- so the GPU and CPU backends agree. 2. The stripped-file fallback returned a fully-realised CuPy DataArray without honouring chunks=, so callers asking for a Dask+CuPy result on a stripped TIFF got an eager array. Mirror the tiled branch's .chunk() step before returning so chunks= works consistently across layouts. Regression tests: - An Orientation=2 stripped TIFF + window= raises ValueError matching the CPU reader. - A stripped TIFF + chunks=N produces a Dask-backed DataArray whose blocks are CuPy arrays of the requested shape; tuple chunks accepted. --- xrspatial/geotiff/__init__.py | 32 ++++++++- .../tests/test_gpu_window_band_1605.py | 72 +++++++++++++++++++ 2 files changed, 102 insertions(+), 2 deletions(-) diff --git a/xrspatial/geotiff/__init__.py b/xrspatial/geotiff/__init__.py index 44bad5c8..d1c1c9c8 100644 --- a/xrspatial/geotiff/__init__.py +++ b/xrspatial/geotiff/__init__.py @@ -2162,6 +2162,22 @@ def read_geotiff_gpu(source: str, *, # only need to copy the post-flip geo_info back here. orientation = ifd.orientation + # Orientation tag (274): values 2-8 mean the stored pixel order + # differs from display order. A windowed read against a non-default + # orientation has ambiguous semantics (does the window refer to + # file pixels or display pixels?), so the CPU reader + # ``_reader.read_to_array`` rejects ``window=`` for orientation != 1. + # Mirror that here so the GPU path agrees with the CPU path and + # ``read_geotiff_dask``. Use the same error wording so the failure + # message is identical across backends. + if orientation != 1 and window is not None: + raise ValueError( + f"Orientation tag (274) is {orientation}; windowed reads " + f"(window=...) and dask-chunked reads (chunks=...) are not " + f"supported for non-default orientation. Read the full " + f"array first, then slice." + ) + # Validate band against the selected IFD's sample count. # ``samples_per_pixel`` is at least 1 for any valid TIFF; we treat # ``band=0`` as "first band" for single-band files too so the @@ -2230,8 +2246,20 @@ def read_geotiff_gpu(source: str, *, coords['band'] = np.arange(arr_gpu.shape[2]) else: dims = ['y', 'x'] - return xr.DataArray(arr_gpu, dims=dims, - coords=coords, name=name, attrs=attrs) + result = xr.DataArray(arr_gpu, dims=dims, + coords=coords, name=name, attrs=attrs) + # ``chunks`` was previously honoured only on the tiled path, + # so stripped TIFFs returned an unchunked DataArray even when + # the caller asked for a Dask+CuPy result. Mirror the tiled + # branch's chunking step so behaviour is consistent across + # layouts. + if chunks is not None: + if isinstance(chunks, int): + chunk_dict = {'y': chunks, 'x': chunks} + else: + chunk_dict = {'y': chunks[0], 'x': chunks[1]} + result = result.chunk(chunk_dict) + return result offsets = ifd.tile_offsets byte_counts = ifd.tile_byte_counts diff --git a/xrspatial/geotiff/tests/test_gpu_window_band_1605.py b/xrspatial/geotiff/tests/test_gpu_window_band_1605.py index 1837782d..f40c455b 100644 --- a/xrspatial/geotiff/tests/test_gpu_window_band_1605.py +++ b/xrspatial/geotiff/tests/test_gpu_window_band_1605.py @@ -194,3 +194,75 @@ def test_read_geotiff_gpu_band_bounds_validation(multi_band_tiff, with pytest.raises(IndexError, match="single-band file"): read_geotiff_gpu(single_path, band=1) + + +@_gpu_only +def test_read_geotiff_gpu_window_rejected_on_nondefault_orientation(tmp_path): + """Mirror the CPU reader: orientation != 1 + window= is ambiguous, raise. + + ``_reader.read_to_array`` rejects this combo with a specific message; + the GPU path must match so backend substitutability holds. The + stripped branch is the easy way to write an Orientation=2 file via + tifffile. + """ + tifffile = pytest.importorskip("tifffile") + from xrspatial.geotiff import read_geotiff_gpu + + arr = np.arange(16 * 20, dtype=np.float32).reshape(16, 20) + path = tmp_path / "orient2_stripped.tif" + # Orientation tag 274 = 2 (top-right, horizontal flip). + tifffile.imwrite( + str(path), arr, + extratags=[(274, 'H', 1, 2, True)], + ) + + with pytest.raises(ValueError, match=r"Orientation tag \(274\) is 2"): + read_geotiff_gpu(str(path), window=(0, 0, 8, 10)) + + +@_gpu_only +def test_read_geotiff_gpu_stripped_chunks_produces_dask(tmp_path): + """Stripped-file branch must honour ``chunks=`` like the tiled branch. + + Pre-fix, the stripped branch returned the eager CuPy DataArray without + calling ``.chunk(...)``, silently dropping ``chunks=`` for stripped + TIFFs. The tiled branch already chunked correctly; this test pins + the parity. + """ + tifffile = pytest.importorskip("tifffile") + import dask.array as dask_array + from xrspatial.geotiff import read_geotiff_gpu + + arr = np.arange(16 * 20, dtype=np.float32).reshape(16, 20) + path = tmp_path / "stripped_chunks.tif" + # No ``tile=`` argument -> stripped layout. Orientation defaults to 1. + tifffile.imwrite(str(path), arr) + + result = read_geotiff_gpu(str(path), chunks=8) + + # The result should be Dask-backed; .data is a dask Array wrapping CuPy. + assert isinstance(result.data, dask_array.Array), \ + f"stripped chunks=8 did not produce a dask-backed DataArray; got " \ + f"{type(result.data).__name__}" + assert result.chunks == ((8, 8), (8, 8, 4)) + # And the underlying chunks are CuPy arrays, not numpy. + block = result.data.blocks[0, 0].compute() + assert type(block).__module__.startswith("cupy"), \ + f"expected CuPy chunk, got {type(block).__module__}" + np.testing.assert_array_equal(result.data.compute().get(), arr) + + +@_gpu_only +def test_read_geotiff_gpu_stripped_chunks_tuple(tmp_path): + """Stripped branch accepts the (rh, cw) tuple chunks spec too.""" + tifffile = pytest.importorskip("tifffile") + import dask.array as dask_array + from xrspatial.geotiff import read_geotiff_gpu + + arr = np.arange(16 * 20, dtype=np.float32).reshape(16, 20) + path = tmp_path / "stripped_chunks_tuple.tif" + tifffile.imwrite(str(path), arr) + + result = read_geotiff_gpu(str(path), chunks=(4, 10)) + assert isinstance(result.data, dask_array.Array) + assert result.chunks == ((4, 4, 4, 4), (10, 10))