Apply TIFF Orientation tag in read_geotiff_gpu (#1540)

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:

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)