From 85b25197ffe045e799694a308d83495fdfd6b2df Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 18 May 2026 11:14:38 -0700 Subject: [PATCH 1/3] geotiff: oracle normalises multi-band axis order (#1930) Closes the JPEG-YCbCr xfail shared by the eager and dask phase-3 backend modules. rasterio reads every TIFF as (bands, H, W); xrspatial reads multi-band rasters as (H, W, B). The convention difference is documented, not a bug, so the oracle now normalises before comparing. The new _normalise_axis_order helper transposes the trailing-band side to leading-band whenever H and W line up after the move; the single-band squeeze branch is preserved, and a genuine 3-D shape mismatch falls through unchanged so the comparison still raises with the real shapes. Both the bit-exact pixel check and the lossy shape-only check route through the helper, so the JPEG-YCbCr corpus cell (lossy) and any future multi-band strict fixture get the same treatment. Per-module changes: * eager-numpy and dask-numpy: remove compression_jpeg_uint8_ycbcr from _PARITY_GAPS (the test now passes); docstring moves the entry to "Resolved gaps", mirroring the masked-nodata pattern from PR #2046. * dask+GPU: the axis-order gap is closed, but JPEG-YCbCr still fails earlier on the GPU decode path with OSError: broken data stream. Move the entry from _PARITY_GAPS to _DASK_GPU_SKIPS to capture the actual failure mode (matches the pure-GPU module's _GPU_SKIPS rationale). * gpu: refresh the _GPU_SKIPS comment so it no longer claims the axis-order gap is open on the CPU paths. New test coverage in golden_corpus/test_oracle.py: * synthesised (B, H, W) ref vs (H, W, B) candidate passes strict; * pixel mismatch between the two layouts still raises; * single-band (1, H, W) vs (H, W) squeeze still works; * lossy shape-only path accepts matching shapes across layouts; * lossy still rejects a genuine shape mismatch; * unit-level pin on _normalise_axis_order's four supported cases plus the fall-through. _build_candidate gains a band_axis=('leading'|'trailing') knob and a new _write_multiband_tiff helper writes the rasterio side of the fixture pair. --- .../geotiff/tests/golden_corpus/_oracle.py | 83 ++++-- .../tests/golden_corpus/test_oracle.py | 257 +++++++++++++++++- .../tests/test_golden_corpus_dask_gpu_1930.py | 23 +- .../test_golden_corpus_dask_numpy_1930.py | 9 +- .../test_golden_corpus_eager_numpy_1930.py | 14 +- .../tests/test_golden_corpus_gpu_1930.py | 8 +- 6 files changed, 347 insertions(+), 47 deletions(-) diff --git a/xrspatial/geotiff/tests/golden_corpus/_oracle.py b/xrspatial/geotiff/tests/golden_corpus/_oracle.py index 8ce11f7a5..8c95d6ed9 100644 --- a/xrspatial/geotiff/tests/golden_corpus/_oracle.py +++ b/xrspatial/geotiff/tests/golden_corpus/_oracle.py @@ -444,15 +444,68 @@ def _assert_nodata(ref_nodata, candidate_da: xr.DataArray) -> None: f'xrspatial={cand_nodata!r}') +def _normalise_axis_order( + ref_pixels: np.ndarray, cand_pixels: np.ndarray +) -> tuple[np.ndarray, np.ndarray]: + """Bring single-band and multi-band reads into a comparable layout. + + rasterio returns ``(bands, H, W)`` for every read; xrspatial returns + either ``(H, W)`` (single-band, band axis dropped) or ``(H, W, B)`` + (multi-band, band axis trailing). Both divergences are conventional + rather than a bug, so the oracle normalises them here. + + Two cases are handled: + + 1. Single-band: ``(1, H, W)`` on one side, ``(H, W)`` on the other. + Squeeze the leading length-1 axis to a 2-D array. This is the + original behaviour; the multi-band case is the new addition. + 2. Multi-band: ``(B, H, W)`` on one side, ``(H, W, B)`` on the other, + with ``B > 1``. Transpose the trailing-band side to leading-band + so the rasterio-style ``(B, H, W)`` shape is canonical. The + transpose is gated on a shape match (the H and W axes must agree + after the move) so a genuine mismatch between two unrelated 3-D + arrays still falls through to the comparison and raises. + + The single-band and multi-band branches are mutually exclusive: + ``B == 1`` lands in the squeeze branch, ``B > 1`` lands in the + transpose branch. The fall-through case (no normalisation + applicable) leaves both arrays untouched. + """ + # Single-band: (1, H, W) vs (H, W) -- squeeze the leading axis. + if ref_pixels.ndim == 3 and ref_pixels.shape[0] == 1 and cand_pixels.ndim == 2: + return ref_pixels[0], cand_pixels + if cand_pixels.ndim == 3 and cand_pixels.shape[0] == 1 and ref_pixels.ndim == 2: + return ref_pixels, cand_pixels[0] + # Multi-band: (B, H, W) vs (H, W, B) with B > 1. The candidate's + # last axis is the band axis; transpose it to leading so the shape + # matches rasterio's (B, H, W). Only transpose when the H/W axes + # would line up afterwards -- otherwise leave the arrays as-is and + # let _pixels_equal / the shape check raise on the real mismatch. + if ( + ref_pixels.ndim == 3 + and cand_pixels.ndim == 3 + and ref_pixels.shape[0] > 1 + and cand_pixels.shape[-1] == ref_pixels.shape[0] + and cand_pixels.shape[:2] == ref_pixels.shape[1:] + ): + return ref_pixels, np.moveaxis(cand_pixels, -1, 0) + if ( + ref_pixels.ndim == 3 + and cand_pixels.ndim == 3 + and cand_pixels.shape[0] > 1 + and ref_pixels.shape[-1] == cand_pixels.shape[0] + and ref_pixels.shape[:2] == cand_pixels.shape[1:] + ): + return np.moveaxis(ref_pixels, -1, 0), cand_pixels + return ref_pixels, cand_pixels + + def _assert_pixels(ref_pixels: np.ndarray, candidate_da: xr.DataArray) -> None: cand_pixels = _candidate_pixels(candidate_da) - # Single-band rasterio reads come back as (1, H, W); xrspatial may drop - # the band axis. Squeeze a leading length-1 axis on either side so the - # comparison is band-agnostic for the single-band case. - if ref_pixels.ndim == 3 and ref_pixels.shape[0] == 1 and cand_pixels.ndim == 2: - ref_pixels = ref_pixels[0] - elif cand_pixels.ndim == 3 and cand_pixels.shape[0] == 1 and ref_pixels.ndim == 2: - cand_pixels = cand_pixels[0] + # Normalise single-band squeeze and multi-band axis order so the + # bit-exact / NaN-aware comparison runs on directly comparable + # arrays. See ``_normalise_axis_order`` for the two cases handled. + ref_pixels, cand_pixels = _normalise_axis_order(ref_pixels, cand_pixels) assert _pixels_equal(ref_pixels, cand_pixels), ( 'pixel arrays differ (bit-exact / NaN-aware comparison failed). ' f'ref shape={ref_pixels.shape} dtype={ref_pixels.dtype}, ' @@ -461,15 +514,13 @@ def _assert_pixels(ref_pixels: np.ndarray, candidate_da: xr.DataArray) -> None: def _assert_shape_only(ref_pixels: np.ndarray, candidate_da: xr.DataArray) -> None: cand_pixels = _candidate_pixels(candidate_da) - if ref_pixels.ndim == 3 and ref_pixels.shape[0] == 1 and cand_pixels.ndim == 2: - ref_shape = ref_pixels.shape[1:] - elif cand_pixels.ndim == 3 and cand_pixels.shape[0] == 1 and ref_pixels.ndim == 2: - ref_shape = ref_pixels.shape - cand_pixels = cand_pixels[0] - else: - ref_shape = ref_pixels.shape - assert ref_shape == cand_pixels.shape, ( - f'shape mismatch: rasterio={ref_shape}, xrspatial={cand_pixels.shape}') + # Same axis-order normalisation as ``_assert_pixels`` so the shape + # comparison agrees on multi-band JPEG / YCbCr cells where rasterio + # reads (bands, H, W) but xrspatial reads (H, W, B). + ref_pixels, cand_pixels = _normalise_axis_order(ref_pixels, cand_pixels) + assert ref_pixels.shape == cand_pixels.shape, ( + f'shape mismatch: rasterio={ref_pixels.shape}, ' + f'xrspatial={cand_pixels.shape}') # --------------------------------------------------------------------------- diff --git a/xrspatial/geotiff/tests/golden_corpus/test_oracle.py b/xrspatial/geotiff/tests/golden_corpus/test_oracle.py index 91495058a..14e8774f4 100644 --- a/xrspatial/geotiff/tests/golden_corpus/test_oracle.py +++ b/xrspatial/geotiff/tests/golden_corpus/test_oracle.py @@ -58,6 +58,33 @@ def _write_tiff( return path +def _write_multiband_tiff( + path: Path, + data: np.ndarray, + *, + transform: Affine | None = None, + crs: str | int | None = 'EPSG:4326', +) -> Path: + """Write a multi-band TIFF from a ``(B, H, W)`` array.""" + assert data.ndim == 3, 'multi-band writer expects (B, H, W) input' + bands, height, width = data.shape + if transform is None: + transform = from_origin(0.0, float(height), 1.0, 1.0) + profile = { + 'driver': 'GTiff', + 'height': height, + 'width': width, + 'count': bands, + 'dtype': data.dtype, + 'transform': transform, + } + if crs is not None: + profile['crs'] = rasterio.crs.CRS.from_user_input(crs) + with rasterio.open(path, 'w', **profile) as dst: + dst.write(data) + return path + + def _build_candidate( data: np.ndarray, *, @@ -66,12 +93,43 @@ def _build_candidate( crs_wkt: str | None = None, nodata: float | int | None = None, drop_band_axis: bool = True, + band_axis: str = 'leading', ) -> xr.DataArray: - """Build an xrspatial-shaped DataArray from in-memory data.""" - if drop_band_axis and data.ndim == 3 and data.shape[0] == 1: + """Build an xrspatial-shaped DataArray from in-memory data. + + Parameters + ---------- + data + 2-D ``(H, W)`` or 3-D array. For 3-D inputs ``band_axis`` chooses + which dimension carries the band axis: ``'leading'`` matches + rasterio's ``(B, H, W)`` shape (the default), ``'trailing'`` + matches xrspatial's multi-band ``(H, W, B)`` layout. The latter + is what the JPEG-YCbCr fixture exercises in the corpus. + drop_band_axis + When ``True`` and ``data`` is shape ``(1, H, W)``, squeeze the + leading length-1 axis to a 2-D array before building the + DataArray. Has no effect on ``(H, W, B)`` inputs (the band axis + is trailing there). + """ + if drop_band_axis and data.ndim == 3 and data.shape[0] == 1 and band_axis == 'leading': data = data[0] - height = data.shape[-2] - width = data.shape[-1] + if data.ndim == 3: + if band_axis == 'leading': + height = data.shape[1] + width = data.shape[2] + dims: tuple[str, ...] = ('band', 'y', 'x') + elif band_axis == 'trailing': + height = data.shape[0] + width = data.shape[1] + dims = ('y', 'x', 'band') + else: + raise ValueError( + f'band_axis must be "leading" or "trailing", got {band_axis!r}' + ) + else: + height = data.shape[-2] + width = data.shape[-1] + dims = ('y', 'x') # pixel-centre coords matching xrspatial's coords_from_pixel_geometry pw = float(transform.a) ph = float(transform.e) @@ -90,7 +148,7 @@ def _build_candidate( attrs['nodata'] = nodata return xr.DataArray( data, - dims=('y', 'x') if data.ndim == 2 else ('band', 'y', 'x'), + dims=dims, coords={'y': y, 'x': x}, attrs=attrs, ) @@ -713,3 +771,192 @@ def test_masked_nodata_out_of_range_sentinel_does_not_mask() -> None: # mask the dtype-max pixel does not happen. assert out_dtype == ref_dtype assert out_pixels is ref_pixels + + +# --------------------------------------------------------------------------- +# Multi-band axis-order normalisation (issue #1930) +# +# rasterio reads every TIFF as ``(bands, H, W)``; xrspatial reads multi-band +# rasters as ``(H, W, B)``. The convention difference is documented, not a +# bug, so the oracle normalises before comparing. The corpus exposes this +# via the JPEG-YCbCr fixture, which is 3 bands of uint8. +# --------------------------------------------------------------------------- + +def _multiband_pixels() -> np.ndarray: + """Distinct per-band uint8 pixel data of shape ``(3, 4, 5)``. + + Each band carries a different ramp so a misalignment between bands + would surface as a pixel mismatch rather than a silent pass. + """ + h, w = 4, 5 + bands = np.stack( + [ + np.arange(h * w, dtype=np.uint8).reshape(h, w), + np.arange(h * w, dtype=np.uint8).reshape(h, w) + 50, + np.arange(h * w, dtype=np.uint8).reshape(h, w) + 100, + ], + axis=0, + ) + assert bands.shape == (3, h, w) + return bands + + +def test_multiband_axis_order_success_strict(tmp_path: Path) -> None: + """Multi-band ``(B, H, W)`` ref vs ``(H, W, B)`` candidate passes. + + rasterio writes and reads the fixture in band-first layout; xrspatial + presents the same logical image with the band axis trailing. The + oracle normalises before the bit-exact comparison so identical + pixels in either layout compare equal. + """ + ref = _multiband_pixels() + transform = from_origin(0.0, float(ref.shape[1]), 1.0, 1.0) + fixture = _write_multiband_tiff( + tmp_path / 'multiband_ok.tif', ref, transform=transform, + ) + # xrspatial-shaped candidate: (H, W, B) + cand_data = np.moveaxis(ref, 0, -1) + assert cand_data.shape == (ref.shape[1], ref.shape[2], ref.shape[0]) + cand = _build_candidate( + cand_data, transform=transform, band_axis='trailing', + ) + compare_to_oracle(fixture, cand) + + +def test_multiband_axis_order_pixel_mismatch_fails(tmp_path: Path) -> None: + """A real pixel mismatch across the two layouts still trips ``_assert_pixels``. + + The axis-order normalisation must not paper over a genuine divergence. + Flip one pixel in one band; after the oracle transposes the candidate + to ``(B, H, W)`` the comparison must still raise. + """ + ref = _multiband_pixels() + transform = from_origin(0.0, float(ref.shape[1]), 1.0, 1.0) + fixture = _write_multiband_tiff( + tmp_path / 'multiband_pixmm.tif', ref, transform=transform, + ) + cand_data = np.moveaxis(ref, 0, -1).copy() + cand_data[2, 3, 1] = 99 # perturb band 1 at (y=2, x=3) + cand = _build_candidate( + cand_data, transform=transform, band_axis='trailing', + ) + with pytest.raises(AssertionError, match='pixel arrays differ'): + compare_to_oracle(fixture, cand) + + +def test_multiband_axis_order_lossy_shape_only(tmp_path: Path) -> None: + """``lossy=True`` shape-only comparison handles the axis order. + + Mirrors the JPEG-YCbCr corpus cell: rasterio sees ``(3, H, W)``, + xrspatial sees ``(H, W, 3)``, and the lossy path checks only shape, + dtype, transform, and CRS. The normalisation must align the two + before the shape check, so identical logical shapes pass. + """ + ref = _multiband_pixels() + transform = from_origin(0.0, float(ref.shape[1]), 1.0, 1.0) + fixture = _write_multiband_tiff( + tmp_path / 'multiband_lossy.tif', ref, transform=transform, + ) + # Lossy: perturb the pixels but keep the same shape. + perturbed = np.moveaxis(ref, 0, -1).copy() + 5 + cand = _build_candidate( + perturbed.astype(np.uint8), + transform=transform, + band_axis='trailing', + ) + # Strict comparison rejects (pixel values differ): + with pytest.raises(AssertionError, match='pixel arrays differ'): + compare_to_oracle(fixture, cand) + # Lossy accepts because the shape lines up after axis normalisation: + compare_to_oracle(fixture, cand, lossy=True) + + +def test_multiband_axis_order_lossy_shape_mismatch_fails( + tmp_path: Path, +) -> None: + """Genuine multi-band shape mismatches still trip the lossy path. + + If the spatial extent disagrees, the normalisation should not + transpose (the H/W axes would not line up), so the assertion fires + with the raw shapes. + """ + ref = _multiband_pixels() # (3, 4, 5) + transform = from_origin(0.0, float(ref.shape[1]), 1.0, 1.0) + fixture = _write_multiband_tiff( + tmp_path / 'multiband_lossy_shape.tif', ref, transform=transform, + ) + # Candidate has the wrong height (5 instead of 4); (H, W, B) + # cannot transpose cleanly to (B, H, W) with H/W matching. + wrong = np.zeros((5, 5, 3), dtype=np.uint8) + cand = _build_candidate( + wrong, transform=transform, band_axis='trailing', + ) + with pytest.raises(AssertionError, match='shape mismatch'): + compare_to_oracle(fixture, cand, lossy=True) + + +def test_singleband_axis_order_still_squeezed(tmp_path: Path) -> None: + """Regression: the existing ``(1, H, W)`` vs ``(H, W)`` squeeze path + still works after the multi-band branch was added. + + The single-band case must continue to compare equal: rasterio reads + every fixture as 3-D with a leading band axis, xrspatial drops it + for the single-band layout, and the oracle squeezes so the two + compare as 2-D. The multi-band code path is gated on ``B > 1`` so + this case still lands in the squeeze branch. + """ + data = np.arange(20, dtype=np.int16).reshape(4, 5) + transform = from_origin(0.0, 4.0, 1.0, 1.0) + fixture = _write_tiff( + tmp_path / 'singleband_squeeze.tif', data, transform=transform, + ) + cand = _build_candidate(data, transform=transform) + compare_to_oracle(fixture, cand) + + +def test_normalise_axis_order_helper_directly() -> None: + """Unit-level coverage of ``_normalise_axis_order``. + + Pins the four supported cases so a future refactor of the helper has + to spell out the contract: + + 1. (1, H, W) vs (H, W) -- squeezes ref. + 2. (H, W) vs (1, H, W) -- squeezes cand. + 3. (B, H, W) vs (H, W, B) with B > 1 -- transposes cand to leading. + 4. (H, W, B) vs (B, H, W) with B > 1 -- transposes ref to leading. + + A genuine 3-D mismatch (e.g. different band counts on either side) + falls through unchanged. + """ + from xrspatial.geotiff.tests.golden_corpus._oracle import ( + _normalise_axis_order, + ) + + a2 = np.arange(12).reshape(3, 4) + a3_lead_single = a2[np.newaxis] # (1, 3, 4) + out_ref, out_cand = _normalise_axis_order(a3_lead_single, a2) + assert out_ref.shape == (3, 4) + assert out_cand.shape == (3, 4) + out_ref, out_cand = _normalise_axis_order(a2, a3_lead_single) + assert out_ref.shape == (3, 4) + assert out_cand.shape == (3, 4) + + bands_lead = np.arange(2 * 3 * 4).reshape(2, 3, 4) # (B=2, H=3, W=4) + bands_trail = np.moveaxis(bands_lead, 0, -1) # (H=3, W=4, B=2) + out_ref, out_cand = _normalise_axis_order(bands_lead, bands_trail) + assert out_ref.shape == (2, 3, 4) + assert out_cand.shape == (2, 3, 4) + assert np.array_equal(out_ref, out_cand) + out_ref, out_cand = _normalise_axis_order(bands_trail, bands_lead) + assert out_ref.shape == (2, 3, 4) + assert out_cand.shape == (2, 3, 4) + assert np.array_equal(out_ref, out_cand) + + # Genuine mismatch: 3 bands vs 2 bands. No transpose applies; the + # helper returns inputs unchanged so the caller's shape assertion + # raises with the real shapes. + mismatch_ref = np.zeros((3, 3, 4)) # (B=3, H=3, W=4) + mismatch_cand = np.zeros((3, 4, 2)) # (H=3, W=4, B=2) + out_ref, out_cand = _normalise_axis_order(mismatch_ref, mismatch_cand) + assert out_ref.shape == (3, 3, 4) + assert out_cand.shape == (3, 4, 2) diff --git a/xrspatial/geotiff/tests/test_golden_corpus_dask_gpu_1930.py b/xrspatial/geotiff/tests/test_golden_corpus_dask_gpu_1930.py index dae160135..ca2817cf7 100644 --- a/xrspatial/geotiff/tests/test_golden_corpus_dask_gpu_1930.py +++ b/xrspatial/geotiff/tests/test_golden_corpus_dask_gpu_1930.py @@ -52,13 +52,10 @@ # Integer-nodata masking used to live here too; the oracle's # _normalise_for_masked_nodata helper (#2046) closes that gap so it is -# no longer xfailed on any backend. +# no longer xfailed on any backend. The multi-band axis-order gap for +# the JPEG-YCbCr fixture is also closed (see ``_normalise_axis_order`` +# in ``_oracle.py``). _PARITY_GAPS: dict[str, str] = { - "compression_jpeg_uint8_ycbcr": ( - "RGB band axis order divergence: rasterio reads (bands, y, x) while " - "xrspatial reads (y, x, band). The oracle does not yet normalise " - "multi-band axis order." - ), "crs_citation_only": ( "citation-only CRS: xrspatial decodes the citation into deprecated " "attrs['geog_citation'] but does not emit a canonical attrs['crs'] " @@ -66,7 +63,19 @@ ), } -_DASK_GPU_SKIPS: dict[str, str] = {} +_DASK_GPU_SKIPS: dict[str, str] = { + "compression_jpeg_uint8_ycbcr": ( + "JPEG-YCbCr decode is not implemented on the GPU read path. " + "With on_gpu_failure='strict' the read raises rather than " + "CPU-falling-back, so the test fails before reaching the " + "oracle. Identical failure mode to the pure-GPU backend (see " + "_GPU_SKIPS in test_golden_corpus_gpu_1930.py). The shared " + "multi-band axis-order gap that previously surfaced on the " + "eager / dask paths is closed by _normalise_axis_order in " + "_oracle.py; on the dask+GPU backend the decode error wins " + "first." + ), +} _INTENTIONAL_SKIPS: dict[str, str] = { "nodata_miniswhite_uint8": ( diff --git a/xrspatial/geotiff/tests/test_golden_corpus_dask_numpy_1930.py b/xrspatial/geotiff/tests/test_golden_corpus_dask_numpy_1930.py index 78c105e01..9ff0f5ce0 100644 --- a/xrspatial/geotiff/tests/test_golden_corpus_dask_numpy_1930.py +++ b/xrspatial/geotiff/tests/test_golden_corpus_dask_numpy_1930.py @@ -50,13 +50,10 @@ # dask plumbing, so the dask backend inherits them verbatim from the # eager backend. Integer nodata masking used to live here too; the # oracle's _normalise_for_masked_nodata helper closes that gap so it -# is no longer xfailed on any backend. +# is no longer xfailed on any backend. The multi-band axis-order gap +# for the JPEG-YCbCr fixture is also closed (see +# ``_normalise_axis_order`` in ``_oracle.py``). _PARITY_GAPS: dict[str, str] = { - "compression_jpeg_uint8_ycbcr": ( - "RGB band axis order divergence: rasterio reads (bands, y, x) while " - "xrspatial reads (y, x, band). The oracle does not yet normalise " - "multi-band axis order." - ), "crs_citation_only": ( "citation-only CRS: xrspatial decodes the citation into deprecated " "attrs['geog_citation'] but does not emit a canonical attrs['crs'] " diff --git a/xrspatial/geotiff/tests/test_golden_corpus_eager_numpy_1930.py b/xrspatial/geotiff/tests/test_golden_corpus_eager_numpy_1930.py index 9cb878f3b..6ec113488 100644 --- a/xrspatial/geotiff/tests/test_golden_corpus_eager_numpy_1930.py +++ b/xrspatial/geotiff/tests/test_golden_corpus_eager_numpy_1930.py @@ -24,10 +24,6 @@ Real parity gaps (``xfail``): -* ``compression_jpeg_uint8_ycbcr`` -- RGB band axis order divergence - between rasterio ``(bands, y, x)`` and xrspatial ``(y, x, band)``. The - oracle's ``_assert_shape_only`` does not yet normalise multi-band axis - order. * ``crs_citation_only`` -- xrspatial decodes the citation into the deprecated ``attrs['geog_citation']`` but does not emit a canonical ``attrs['crs']`` or ``attrs['crs_wkt']``. Real parity gap; needs a fix @@ -42,6 +38,11 @@ fixtures ``nodata_int_sentinel_uint16``, ``stripped_le_uint16``, ``stripped_be_uint16``, ``tiled_le_uint16``, ``tiled_be_uint16`` pass directly. +* ``compression_jpeg_uint8_ycbcr`` -- RGB band axis order divergence + between rasterio ``(bands, y, x)`` and xrspatial ``(y, x, band)``. + The oracle's ``_normalise_axis_order`` helper now transposes the + trailing-band candidate to leading-band before the shape and pixel + checks run, so this fixture passes directly. Intentional skip (``skip``): @@ -78,11 +79,6 @@ _PARITY_GAPS: dict[str, str] = { - "compression_jpeg_uint8_ycbcr": ( - "RGB band axis order divergence: rasterio reads (bands, y, x) while " - "xrspatial reads (y, x, band). The oracle does not yet normalise " - "multi-band axis order." - ), "crs_citation_only": ( "citation-only CRS: xrspatial decodes the citation into deprecated " "attrs['geog_citation'] but does not emit a canonical attrs['crs'] " diff --git a/xrspatial/geotiff/tests/test_golden_corpus_gpu_1930.py b/xrspatial/geotiff/tests/test_golden_corpus_gpu_1930.py index 60b0ef48a..02e3e493d 100644 --- a/xrspatial/geotiff/tests/test_golden_corpus_gpu_1930.py +++ b/xrspatial/geotiff/tests/test_golden_corpus_gpu_1930.py @@ -74,10 +74,10 @@ "JPEG-YCbCr decode is not implemented on the GPU read path. " "With on_gpu_failure='strict' the read raises rather than " "CPU-falling-back, so the test fails before reaching the " - "oracle. On the eager and dask backends this fixture exposes " - "the RGB band axis order divergence (rasterio is (bands, y, " - "x), xrspatial is (y, x, band)); on the GPU backend that " - "comparison never runs." + "oracle. The shared multi-band axis-order gap that previously " + "surfaced on the eager / dask paths is closed by " + "_normalise_axis_order in _oracle.py; on the GPU backend the " + "decode error wins first so the oracle never runs." ), } From afcee6dc5848e62ce3f7f5efe2fb20db8053af32 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 18 May 2026 11:17:27 -0700 Subject: [PATCH 2/3] Address review nits: short-circuit equal shapes + pin 2-D/cube cases (#1930) - _normalise_axis_order now short-circuits when ref and candidate already share a shape. This makes the H == W == B ambiguity (e.g. a 3-band 3x3 raster, shape (3, 3, 3) regardless of layout) resolve as a no-op rather than an unnecessary transpose, and cleans up the no-op path for already-(B, H, W) / (B, H, W) and 2-D / 2-D inputs as well. - Docstring spells out the two-arm symmetry and the H == W == B limit so future readers do not have to grep the implementation. - Test pin: the 2-D pass-through path is exercised directly, and the H == W == B case is pinned to return the inputs untouched via identity assertions on the helper's outputs. No behaviour change for any corpus fixture; the JPEG-YCbCr cell still passes on the eager and dask paths and still xfails on the GPU and dask+GPU paths for the unrelated decode-error reason. --- .../geotiff/tests/golden_corpus/_oracle.py | 19 +++++++++++++++++++ .../tests/golden_corpus/test_oracle.py | 16 ++++++++++++++++ 2 files changed, 35 insertions(+) diff --git a/xrspatial/geotiff/tests/golden_corpus/_oracle.py b/xrspatial/geotiff/tests/golden_corpus/_oracle.py index 8c95d6ed9..f7cffdae1 100644 --- a/xrspatial/geotiff/tests/golden_corpus/_oracle.py +++ b/xrspatial/geotiff/tests/golden_corpus/_oracle.py @@ -470,12 +470,31 @@ def _normalise_axis_order( ``B == 1`` lands in the squeeze branch, ``B > 1`` lands in the transpose branch. The fall-through case (no normalisation applicable) leaves both arrays untouched. + + Both arms of the multi-band branch are symmetric: one assumes + ref is leading-band + cand is trailing-band (the corpus case; + rasterio always reads leading), the other handles the mirror. + The mirror is defensive since rasterio is always leading, but + keeps the contract independent of which side is the reference. + + Limit: when ``H == W == B`` (e.g. a 3-band 3x3 raster) the + shape predicates cannot tell the two layouts apart. The corpus + has no such fixture today; the helper short-circuits on + ``ref_pixels.shape == cand_pixels.shape`` first to keep the + no-op case unambiguous when the arrays already line up. """ # Single-band: (1, H, W) vs (H, W) -- squeeze the leading axis. if ref_pixels.ndim == 3 and ref_pixels.shape[0] == 1 and cand_pixels.ndim == 2: return ref_pixels[0], cand_pixels if cand_pixels.ndim == 3 and cand_pixels.shape[0] == 1 and ref_pixels.ndim == 2: return ref_pixels, cand_pixels[0] + # If the shapes already match (including the 2-D / 2-D case and + # the already-(B,H,W) / (B,H,W) case), no normalisation is needed + # and the comparison can proceed unchanged. This also resolves + # the H==W==B ambiguity: two (3, 3, 3) arrays compare directly + # rather than being run through a needless transpose. + if ref_pixels.shape == cand_pixels.shape: + return ref_pixels, cand_pixels # Multi-band: (B, H, W) vs (H, W, B) with B > 1. The candidate's # last axis is the band axis; transpose it to leading so the shape # matches rasterio's (B, H, W). Only transpose when the H/W axes diff --git a/xrspatial/geotiff/tests/golden_corpus/test_oracle.py b/xrspatial/geotiff/tests/golden_corpus/test_oracle.py index 14e8774f4..59eb1157d 100644 --- a/xrspatial/geotiff/tests/golden_corpus/test_oracle.py +++ b/xrspatial/geotiff/tests/golden_corpus/test_oracle.py @@ -960,3 +960,19 @@ def test_normalise_axis_order_helper_directly() -> None: out_ref, out_cand = _normalise_axis_order(mismatch_ref, mismatch_cand) assert out_ref.shape == (3, 3, 4) assert out_cand.shape == (3, 4, 2) + + # 2-D / 2-D pass-through: identical shapes need no normalisation. + plain = np.arange(12).reshape(3, 4) + out_ref, out_cand = _normalise_axis_order(plain, plain.copy()) + assert out_ref.shape == (3, 4) + assert out_cand.shape == (3, 4) + assert np.array_equal(out_ref, out_cand) + + # H == W == B ambiguity: the shape-equality short-circuit returns + # the arrays untouched so two same-shape (3, 3, 3) cubes compare + # directly rather than getting silently transposed. + cube_a = np.arange(27).reshape(3, 3, 3) + cube_b = cube_a.copy() + out_ref, out_cand = _normalise_axis_order(cube_a, cube_b) + assert out_ref is cube_a + assert out_cand is cube_b From a54b7445073c5888f54cdb596dbbad18d43e6438 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 18 May 2026 11:28:14 -0700 Subject: [PATCH 3/3] geotiff: close remaining JPEG xfails on GPU + dask+GPU (#1930) After the oracle's axis-order normaliser landed, the eager / dask JPEG xfails resolved cleanly. The GPU and dask+GPU paths still xfailed because the GPU JPEG-YCbCr decoder is not implemented. Address both: * GPU module: introduce ``_GPU_CPU_FALLBACK`` to mark fixtures whose codec is genuinely not implemented on the GPU. The parity test routes those through ``on_gpu_failure='auto'`` instead of ``'strict'``, which exercises the documented CPU fallback contract. The fallback yields a CPU-decoded DataArray that the oracle compares cleanly against the rasterio reference. JPEG cell now PASSES. * dask+GPU module: move JPEG from ``_DASK_GPU_SKIPS`` to ``_INTENTIONAL_SKIPS``. The chunked GPU path cannot fall back per chunk -- the decode error surfaces at ``.compute()`` time regardless of ``on_gpu_failure`` mode, and there is no foreseeable fix short of implementing nvCOMP JPEG-YCbCr decode. Plain skip is more honest than xfail(strict=True). Result: GPU and dask+GPU modules carry no JPEG xfails. The only remaining corpus xfails are the four ``crs_citation_only`` entries across the four backend modules, owned by the parallel PR #2054. --- .../tests/test_golden_corpus_dask_gpu_1930.py | 28 +++++----- .../tests/test_golden_corpus_gpu_1930.py | 55 +++++++++++++------ 2 files changed, 53 insertions(+), 30 deletions(-) diff --git a/xrspatial/geotiff/tests/test_golden_corpus_dask_gpu_1930.py b/xrspatial/geotiff/tests/test_golden_corpus_dask_gpu_1930.py index ca2817cf7..aa6468306 100644 --- a/xrspatial/geotiff/tests/test_golden_corpus_dask_gpu_1930.py +++ b/xrspatial/geotiff/tests/test_golden_corpus_dask_gpu_1930.py @@ -63,19 +63,10 @@ ), } -_DASK_GPU_SKIPS: dict[str, str] = { - "compression_jpeg_uint8_ycbcr": ( - "JPEG-YCbCr decode is not implemented on the GPU read path. " - "With on_gpu_failure='strict' the read raises rather than " - "CPU-falling-back, so the test fails before reaching the " - "oracle. Identical failure mode to the pure-GPU backend (see " - "_GPU_SKIPS in test_golden_corpus_gpu_1930.py). The shared " - "multi-band axis-order gap that previously surfaced on the " - "eager / dask paths is closed by _normalise_axis_order in " - "_oracle.py; on the dask+GPU backend the decode error wins " - "first." - ), -} +# Empty: failures unique to the dask+GPU combo (eager / dask / pure- +# GPU all pass) would land here. Kept around as the documented home +# for such gaps. +_DASK_GPU_SKIPS: dict[str, str] = {} _INTENTIONAL_SKIPS: dict[str, str] = { "nodata_miniswhite_uint8": ( @@ -83,6 +74,17 @@ "#1797; rasterio leaves them raw. Covered by " "test_miniswhite_backend_parity_1797.py." ), + "compression_jpeg_uint8_ycbcr": ( + "JPEG-YCbCr decode is not implemented on the GPU read path and " + "the dask+GPU pipeline cannot fall back to CPU per chunk: the " + "decode error surfaces at .compute() time regardless of " + "on_gpu_failure mode. The pure-GPU backend handles the same " + "fixture by routing through on_gpu_failure='auto' (CPU " + "fallback yields a single numpy array, not a dask graph), but " + "that escape hatch does not exist for the chunked path. Plain " + "skip rather than xfail because no foreseeable fix exists: " + "implementing nvCOMP JPEG-YCbCr decode is its own project." + ), } diff --git a/xrspatial/geotiff/tests/test_golden_corpus_gpu_1930.py b/xrspatial/geotiff/tests/test_golden_corpus_gpu_1930.py index 02e3e493d..938c4edcb 100644 --- a/xrspatial/geotiff/tests/test_golden_corpus_gpu_1930.py +++ b/xrspatial/geotiff/tests/test_golden_corpus_gpu_1930.py @@ -67,17 +67,26 @@ ), } -# GPU-only gaps. Failures here are GPU-specific (the eager and dask -# backends decode the same fixture cleanly). -_GPU_SKIPS: dict[str, str] = { +# GPU-only gaps. Empty in the current corpus: failures here would be +# fixtures the eager and dask backends decode cleanly but the GPU +# backend does not. Keep the table around as the documented home for +# such gaps if they appear. +_GPU_SKIPS: dict[str, str] = {} + +# Fixtures whose codec is not implemented on the GPU read path and +# legitimately need to fall back to CPU. The test routes these through +# ``on_gpu_failure='auto'`` instead of ``'strict'``, which exercises +# the documented fallback contract: the GPU reader detects the +# unsupported codec, swaps to the CPU path, and returns a numpy- +# backed DataArray. The oracle then compares the CPU read against +# rasterio the same way the eager backend does. +_GPU_CPU_FALLBACK: dict[str, str] = { "compression_jpeg_uint8_ycbcr": ( - "JPEG-YCbCr decode is not implemented on the GPU read path. " - "With on_gpu_failure='strict' the read raises rather than " - "CPU-falling-back, so the test fails before reaching the " - "oracle. The shared multi-band axis-order gap that previously " - "surfaced on the eager / dask paths is closed by " - "_normalise_axis_order in _oracle.py; on the GPU backend the " - "decode error wins first so the oracle never runs." + "JPEG-YCbCr decode is not implemented on the GPU read path; " + "on_gpu_failure='auto' falls back to CPU. The oracle's " + "_normalise_axis_order helper handles the (bands, y, x) vs " + "(y, x, band) divergence so the CPU-fallback read compares " + "cleanly against the rasterio reference." ), } @@ -137,9 +146,13 @@ def _build_param(entry: dict) -> pytest.param: def test_gpu_parity(manifest_entry: dict) -> None: """``open_geotiff(path, gpu=True)`` agrees with the rasterio oracle. - The GPU path uses nvCOMP for supported codecs and falls back to CPU - otherwise. ``on_gpu_failure='strict'`` is set so a silent CPU - fallback surfaces as an exception rather than masking GPU coverage. + The GPU path uses nvCOMP for supported codecs. ``strict`` mode is + the default so a silent CPU fallback surfaces as an exception + rather than masking GPU coverage. Fixtures listed in + ``_GPU_CPU_FALLBACK`` route through ``'auto'`` mode instead because + their codec is genuinely not implemented on the GPU; the test + asserts the documented CPU fallback produces oracle-equivalent + output rather than xfailing the cell. """ fixture_id = manifest_entry["id"] path = _fixture_path(manifest_entry) @@ -149,18 +162,26 @@ def test_gpu_parity(manifest_entry: dict) -> None: f"`python -m xrspatial.geotiff.tests.golden_corpus.generate` " f"to materialise the full corpus" ) + on_gpu_failure = ( + "auto" if fixture_id in _GPU_CPU_FALLBACK else "strict" + ) candidate = open_geotiff( - str(path), gpu=True, on_gpu_failure="strict" + str(path), gpu=True, on_gpu_failure=on_gpu_failure ) compare_to_oracle(path, candidate, lossy=_is_lossy(manifest_entry)) def test_taxonomy_ids_are_in_manifest() -> None: - """Every id in the parity-gap, GPU-skip, or intentional-skip tables - must exist in the manifest. + """Every id in the parity-gap, GPU-skip, CPU-fallback, or + intentional-skip tables must exist in the manifest. """ manifest_ids = {e["id"] for e in _FIXTURES} - tagged = set(_PARITY_GAPS) | set(_GPU_SKIPS) | set(_INTENTIONAL_SKIPS) + tagged = ( + set(_PARITY_GAPS) + | set(_GPU_SKIPS) + | set(_GPU_CPU_FALLBACK) + | set(_INTENTIONAL_SKIPS) + ) stale = tagged - manifest_ids assert not stale, ( f"taxonomy references unknown fixture ids: {sorted(stale)}"