From cd4f3301e43e1e0be5ec2a8079ecb4373aa738f4 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sat, 9 May 2026 18:21:48 -0700 Subject: [PATCH] Apply TIFF Orientation tag in read_geotiff_gpu (#1540) PR #1521 wired up the Orientation tag (274) on the CPU read path but left the GPU path on the raw stored buffer, so any TIFF written with orientation 2-8 decoded to different pixel values on CPU vs GPU. This adds two helpers in xrspatial/geotiff/__init__.py: - _apply_orientation_gpu mirrors the CPU _apply_orientation slicing logic against a cupy ndarray (eight orientations, 2-D and 3-D variants). - _apply_orientation_geo_info centralises the post-flip transform update so both read_to_array (CPU) and read_geotiff_gpu update GeoTransform consistently. read_geotiff_gpu now applies orientation post-decode for every pure-GPU code path, and reuses the geo_info returned by read_to_array when the CPU fallback path runs (since that path already applies the remap). A flag (arr_was_cpu_decoded) keeps us from double-applying. Test coverage: test_orientation_gpu.py exercises every orientation on single-band tiled, single-band stripped, and 3-band planar=2 tiled files, plus georef sel-fidelity for the mirror-flip cases. --- xrspatial/geotiff/__init__.py | 149 +++++++++++++- .../geotiff/tests/test_orientation_gpu.py | 182 ++++++++++++++++++ 2 files changed, 325 insertions(+), 6 deletions(-) create mode 100644 xrspatial/geotiff/tests/test_orientation_gpu.py diff --git a/xrspatial/geotiff/__init__.py b/xrspatial/geotiff/__init__.py index 4b27f245..bff7483e 100644 --- a/xrspatial/geotiff/__init__.py +++ b/xrspatial/geotiff/__init__.py @@ -1544,6 +1544,107 @@ def _gpu_decode_single_band_tiles( return arr_gpu +def _apply_orientation_gpu(arr_gpu, orientation: int): + """cupy-side mirror of :func:`._reader._apply_orientation`. + + The CPU reader applies the TIFF Orientation tag (274) post-decode so + pixel (0, 0) is always the visual top-left. The GPU read path used + to skip this remap, so reads of any file with orientation != 1 + returned different pixel buffers than the CPU reader (#1540). + + Same eight orientations the CPU helper handles. Operates on a cupy + ndarray and returns a cupy ndarray; ``cupy.ascontiguousarray`` is + applied so downstream views (DataArray.data) work without surprise + re-strides on the GPU. + """ + import cupy + if orientation == 1: + return arr_gpu + if orientation == 2: + return cupy.ascontiguousarray(arr_gpu[:, ::-1]) + if orientation == 3: + return cupy.ascontiguousarray(arr_gpu[::-1, ::-1]) + if orientation == 4: + return cupy.ascontiguousarray(arr_gpu[::-1, :]) + if arr_gpu.ndim == 3: + if orientation == 5: + return cupy.ascontiguousarray(arr_gpu.transpose(1, 0, 2)) + if orientation == 6: + return cupy.ascontiguousarray(arr_gpu.transpose(1, 0, 2)[:, ::-1]) + if orientation == 7: + return cupy.ascontiguousarray( + arr_gpu.transpose(1, 0, 2)[::-1, ::-1]) + if orientation == 8: + return cupy.ascontiguousarray(arr_gpu.transpose(1, 0, 2)[::-1, :]) + else: + if orientation == 5: + return cupy.ascontiguousarray(arr_gpu.T) + if orientation == 6: + return cupy.ascontiguousarray(arr_gpu.T[:, ::-1]) + if orientation == 7: + return cupy.ascontiguousarray(arr_gpu.T[::-1, ::-1]) + if orientation == 8: + return cupy.ascontiguousarray(arr_gpu.T[::-1, :]) + raise ValueError( + f"Invalid TIFF Orientation tag value: {orientation} " + f"(must be 1-8 per TIFF 6.0)" + ) + + +def _apply_orientation_geo_info(geo_info, orientation: int, + file_h: int, file_w: int): + """Mirror the transform updates `_reader.read_to_array` does post-flip. + + Centralised so both ``read_to_array`` (CPU) and ``read_geotiff_gpu`` + (this module) update the GeoTransform consistently. Operates only + on ``geo_info.transform``; the rest of the GeoInfo struct stays as + parsed. + """ + if orientation == 1 or geo_info is None or geo_info.transform is None: + return geo_info + t = geo_info.transform + if orientation in (2, 3, 4): + if geo_info.raster_type == RASTER_PIXEL_IS_POINT: + x_shift = file_w - 1 + y_shift = file_h - 1 + else: + x_shift = file_w + y_shift = file_h + new_origin_x = t.origin_x + new_origin_y = t.origin_y + new_px_w = t.pixel_width + new_px_h = t.pixel_height + if orientation in (2, 3): + new_origin_x = t.origin_x + x_shift * t.pixel_width + new_px_w = -t.pixel_width + if orientation in (3, 4): + new_origin_y = t.origin_y + y_shift * t.pixel_height + new_px_h = -t.pixel_height + geo_info.transform = GeoTransform( + origin_x=new_origin_x, + origin_y=new_origin_y, + pixel_width=new_px_w, + pixel_height=new_px_h, + ) + elif orientation in (5, 6, 7, 8): + geo_info.transform = GeoTransform( + origin_x=t.origin_x, + origin_y=t.origin_y, + pixel_width=t.pixel_height, + pixel_height=t.pixel_width, + ) + if (geo_info.crs_epsg is not None + or geo_info.crs_wkt is not None): + warnings.warn( + f"Orientation {orientation} swaps spatial axes on " + f"a georeferenced file; the returned coords are " + f"shape-correct but the geographic transform may " + f"need manual adjustment.", + stacklevel=3, + ) + return geo_info + + def read_geotiff_gpu(source: str, *, dtype=None, overview_level: int | None = None, @@ -1649,11 +1750,21 @@ def read_geotiff_gpu(source: str, *, bps = resolve_bits_per_sample(ifd.bits_per_sample) file_dtype = tiff_dtype_to_numpy(bps, ifd.sample_format) geo_info = extract_geo_info(ifd, data, header.byte_order) + # Capture the Orientation tag (274) once so the post-decode flip + # below picks it up for both the stripped fallback and the tiled + # GPU pipelines. CPU read_to_array applies the array remap + + # transform update for stripped reads, so for that branch we + # only need to copy the post-flip geo_info back here. + orientation = ifd.orientation if not ifd.is_tiled: - # Fall back to CPU for stripped files + # Fall back to CPU for stripped files. read_to_array already + # applies the orientation remap to both the array and the + # geo_info transform (#1537), so reuse its geo_info here + # rather than the pre-orientation one we just extracted. src.close() - arr_cpu, _ = read_to_array(source, overview_level=overview_level) + arr_cpu, geo_info = read_to_array( + source, overview_level=overview_level) arr_gpu = cupy.asarray(arr_cpu) coords = _geo_to_coords(geo_info, arr_gpu.shape[0], arr_gpu.shape[1]) if name is None: @@ -1722,6 +1833,12 @@ def read_geotiff_gpu(source: str, *, masked_fill = (_resolve_masked_fill(ifd.nodata_str, file_dtype) if compression == COMPRESSION_LERC else None) + # Track whether the array we end up with was already orientation-flipped + # by `read_to_array`. Any path that falls back to CPU decode picks up + # the orientation remap from PR #1521 + #1537 for free; the pure GPU + # paths still need the explicit remap added in #1540. + arr_was_cpu_decoded = False + # PlanarConfiguration=2 (separate bands): each band has its own list # of tiles back-to-back in TileOffsets / TileByteCounts. The GPU # tile-assembly kernel assumes a single chunky tile sequence with @@ -1784,8 +1901,10 @@ def _read_once(): break band_arrays.append(band_arr) if cpu_fallback_needed: - arr_cpu, _ = read_to_array(source, overview_level=overview_level) + arr_cpu, geo_info = read_to_array( + source, overview_level=overview_level) arr_gpu = cupy.asarray(arr_cpu) + arr_was_cpu_decoded = True else: arr_gpu = cupy.stack(band_arrays, axis=2) if arr_gpu.shape != (height, width, samples): @@ -1795,8 +1914,10 @@ def _read_once(): f"({height}, {width}, {samples})" ) elif has_sparse_tile: - arr_cpu, _ = read_to_array(source, overview_level=overview_level) + arr_cpu, geo_info = read_to_array( + source, overview_level=overview_level) arr_gpu = cupy.asarray(arr_cpu) + arr_was_cpu_decoded = True else: from ._gpu_decode import gpu_decode_tiles_from_file arr_gpu = None @@ -1850,8 +1971,10 @@ def _read_once(): RuntimeWarning, stacklevel=2, ) - arr_cpu, _ = read_to_array(source, overview_level=overview_level) + arr_cpu, geo_info = read_to_array( + source, overview_level=overview_level) arr_gpu = cupy.asarray(arr_cpu) + arr_was_cpu_decoded = True # Multi-band tiled output must be (H, W, samples) regardless of planar # config -- catch any shape regression in the kernels before we attach @@ -1866,6 +1989,17 @@ def _read_once(): f"({height}, {width}, {samples})" ) + # Apply the TIFF Orientation tag (274) to the GPU array. The CPU + # reader does this in `read_to_array` (#1521 + #1537), so any path + # that already went through the CPU fallback above is a no-op here. + # The pure GPU paths land at this point with a raw stored-order + # buffer, and without this remap they would silently disagree with + # the CPU reader (#1540). + if orientation != 1 and not arr_was_cpu_decoded: + arr_gpu = _apply_orientation_gpu(arr_gpu, orientation) + geo_info = _apply_orientation_geo_info( + geo_info, orientation, file_h=height, file_w=width) + if dtype is not None: target = np.dtype(dtype) _validate_dtype_cast(np.dtype(str(arr_gpu.dtype)), target) @@ -1876,7 +2010,10 @@ def _read_once(): import os name = os.path.splitext(os.path.basename(source))[0] - coords = _geo_to_coords(geo_info, height, width) + # 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 = {} if geo_info.crs_epsg is not None: diff --git a/xrspatial/geotiff/tests/test_orientation_gpu.py b/xrspatial/geotiff/tests/test_orientation_gpu.py new file mode 100644 index 00000000..1fc8dc70 --- /dev/null +++ b/xrspatial/geotiff/tests/test_orientation_gpu.py @@ -0,0 +1,182 @@ +"""GPU follow-up to PR #1521 (TIFF Orientation tag on read). + +The CPU reader applies the Orientation tag (274) post-decode so pixel +(0, 0) is always the visual top-left. The GPU read path used to skip +this remap, so reads of any file with orientation != 1 returned a +different pixel buffer than the CPU reader (#1540). +""" +from __future__ import annotations + +import importlib.util +import warnings as _warnings + +import numpy as np +import pytest + +tifffile = pytest.importorskip("tifffile") + + +def _gpu_available() -> bool: + """True if cupy is importable and CUDA is initialised.""" + 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", +) + + +_ORIENTATIONS = [1, 2, 3, 4, 5, 6, 7, 8] + + +def _expected_for_orientation(stored, orientation): + """Spec-defined remap so tests don't depend on the production helper.""" + if orientation == 1: + return stored + if orientation == 2: + return stored[:, ::-1] + if orientation == 3: + return stored[::-1, ::-1] + if orientation == 4: + return stored[::-1, :] + if orientation == 5: + return stored.T if stored.ndim == 2 else stored.transpose(1, 0, 2) + if orientation == 6: + return (stored.T[:, ::-1] if stored.ndim == 2 + else stored.transpose(1, 0, 2)[:, ::-1]) + if orientation == 7: + return (stored.T[::-1, ::-1] if stored.ndim == 2 + else stored.transpose(1, 0, 2)[::-1, ::-1]) + if orientation == 8: + return (stored.T[::-1, :] if stored.ndim == 2 + else stored.transpose(1, 0, 2)[::-1, :]) + raise AssertionError(orientation) + + +def _write_tiled(path, arr, orientation, tile=(16, 16), extra=None): + """Write *arr* tiled with the requested Orientation tag.""" + extras = [(274, 'H', 1, orientation, True)] + if extra: + extras.extend(extra) + tifffile.imwrite(str(path), arr, tile=tile, extratags=extras) + + +def _write_stripped(path, arr, orientation, extra=None): + """Write *arr* stripped (no tile=) with the requested Orientation tag.""" + extras = [(274, 'H', 1, orientation, True)] + if extra: + extras.extend(extra) + tifffile.imwrite(str(path), arr, extratags=extras) + + +@_gpu_only +@pytest.mark.parametrize("orientation", _ORIENTATIONS) +def test_gpu_tiled_matches_cpu(tmp_path, orientation): + """Tiled GPU read of every orientation matches the spec-remapped buffer.""" + from xrspatial.geotiff import read_geotiff_gpu + + arr = np.arange(256, dtype=np.uint8).reshape(16, 16) + path = tmp_path / f"gpu_orient_{orientation}.tif" + _write_tiled(path, arr, orientation) + + da = read_geotiff_gpu(str(path)) + expected = _expected_for_orientation(arr, orientation) + np.testing.assert_array_equal(da.data.get(), expected) + + +@_gpu_only +@pytest.mark.parametrize("orientation", _ORIENTATIONS) +def test_gpu_stripped_matches_cpu(tmp_path, orientation): + """Stripped GPU read also applies orientation (via CPU fallback path).""" + from xrspatial.geotiff import read_geotiff_gpu + + arr = np.arange(64, dtype=np.uint8).reshape(8, 8) + path = tmp_path / f"gpu_strip_orient_{orientation}.tif" + _write_stripped(path, arr, orientation) + + da = read_geotiff_gpu(str(path)) + expected = _expected_for_orientation(arr, orientation) + np.testing.assert_array_equal(da.data.get(), expected) + + +@_gpu_only +@pytest.mark.parametrize("orientation", _ORIENTATIONS) +def test_gpu_3band_tiled_matches_cpu(tmp_path, orientation): + """3-band tiled (planar=2) read also applies orientation per band.""" + from xrspatial.geotiff import read_geotiff_gpu + + rgb = np.arange(3 * 16 * 16, dtype=np.uint8).reshape(3, 16, 16) + path = tmp_path / f"gpu_rgb_orient_{orientation}.tif" + tifffile.imwrite( + str(path), rgb, photometric='rgb', tile=(16, 16), + extratags=[(274, 'H', 1, orientation, True)], + ) + + da = read_geotiff_gpu(str(path)) + # tifffile writes (bands, h, w) -> on disk planar=2; reader returns (y, x, band) + stored = np.transpose(rgb, (1, 2, 0)) + expected = _expected_for_orientation(stored, orientation) + np.testing.assert_array_equal(da.data.get(), expected) + + +@_gpu_only +@pytest.mark.parametrize("orientation", [2, 3, 4]) +def test_gpu_orient_2_3_4_coords_track_pixel_flip(tmp_path, orientation): + """For mirror-flip orientations, GPU coord array also flips.""" + from xrspatial.geotiff import read_geotiff_gpu + + arr = np.arange(256, dtype=np.uint8).reshape(16, 16) + path = tmp_path / f"gpu_orient_geo_{orientation}.tif" + _write_tiled( + path, arr, orientation, + extra=[ + (33550, 'd', 3, (1.0, 1.0, 0.0)), + (33922, 'd', 6, (0.0, 0.0, 0.0, 100.0, 50.0, 0.0)), + (34735, 'H', 12, ( + 1, 1, 0, 2, + 1024, 0, 1, 2, + 2048, 0, 1, 4326, + )), + ], + ) + + with _warnings.catch_warnings(): + _warnings.simplefilter('ignore') + da = read_geotiff_gpu(str(path)) + + # Pixel (0, 0) of the original file is value 0 at (x=100.5, y=49.5). + # Pixel (0, 15) is value 15 at (x=115.5, y=49.5). + # Pixel (15, 0) is value 240 at (x=100.5, y=34.5). + targets = [ + (100.5, 49.5, 0), + (115.5, 49.5, 15), + (100.5, 34.5, 240), + (115.5, 34.5, 255), + ] + for x, y, expected in targets: + got = int(da.sel(x=x, y=y).item()) + assert got == expected, ( + f"orient={orientation}: GPU sel(x={x}, y={y})={got}, " + f"expected {expected}" + ) + + +@_gpu_only +def test_gpu_default_orientation_unchanged(tmp_path): + """Files without Orientation tag still decode unchanged on GPU.""" + from xrspatial.geotiff import read_geotiff_gpu + + arr = np.arange(256, dtype=np.uint8).reshape(16, 16) + path = tmp_path / "gpu_no_orient.tif" + tifffile.imwrite(str(path), arr, tile=(16, 16)) + + da = read_geotiff_gpu(str(path)) + np.testing.assert_array_equal(da.data.get(), arr)