diff --git a/xrspatial/geotiff/_header.py b/xrspatial/geotiff/_header.py index 5f6f6a1b..e061d41c 100644 --- a/xrspatial/geotiff/_header.py +++ b/xrspatial/geotiff/_header.py @@ -23,6 +23,7 @@ TAG_COMPRESSION = 259 TAG_PHOTOMETRIC = 262 TAG_STRIP_OFFSETS = 273 +TAG_ORIENTATION = 274 TAG_SAMPLES_PER_PIXEL = 277 TAG_ROWS_PER_STRIP = 278 TAG_STRIP_BYTE_COUNTS = 279 @@ -162,6 +163,20 @@ def tile_byte_counts(self) -> tuple | None: def photometric(self) -> int: return self.get_value(TAG_PHOTOMETRIC, 1) + @property + def orientation(self) -> int: + """Orientation tag (274). Default 1 = top-left (no transform). + + Per TIFF 6.0 the eight valid values are: + 1=top-left, 2=top-right, 3=bottom-right, 4=bottom-left, + 5=left-top, 6=right-top, 7=right-bottom, 8=left-bottom. + Values 5-8 swap rows and columns relative to the stored layout. + """ + v = self.get_value(TAG_ORIENTATION, 1) + if isinstance(v, tuple): + v = v[0] + return int(v) + @property def planar_config(self) -> int: return self.get_value(TAG_PLANAR_CONFIG, 1) diff --git a/xrspatial/geotiff/_reader.py b/xrspatial/geotiff/_reader.py index e7fec46e..bb1eeed5 100644 --- a/xrspatial/geotiff/_reader.py +++ b/xrspatial/geotiff/_reader.py @@ -1085,6 +1085,65 @@ def _read_cog_http(url: str, overview_level: int | None = None, # Main read function # --------------------------------------------------------------------------- +def _apply_orientation(arr: np.ndarray, orientation: int) -> np.ndarray: + """Reorient a decoded TIFF array according to the Orientation tag (274). + + The TIFF 6.0 spec defines eight orientations describing where the + *first row* and *first column* of the stored data sit relative to the + visual top-left of the image: + + === ================= ======================================== + 1 top-left identity (default, no transform) + 2 top-right mirror horizontally (flip columns) + 3 bottom-right rotate 180 degrees + 4 bottom-left mirror vertically (flip rows) + 5 left-top transpose (rows<->columns) + 6 right-top rotate 90 clockwise + 7 right-bottom transverse (anti-transpose) + 8 left-bottom rotate 90 counter-clockwise + === ================= ======================================== + + Values 5-8 swap rows and columns: the file's stored width becomes the + output's height and vice versa. + + The input ``arr`` is shaped ``(height, width)`` or + ``(height, width, samples)``. Multi-band 3D arrays only have their + first two axes transformed; the sample axis is preserved. + """ + if orientation == 1: + return arr + if orientation == 2: + return np.ascontiguousarray(arr[:, ::-1]) + if orientation == 3: + return np.ascontiguousarray(arr[::-1, ::-1]) + if orientation == 4: + return np.ascontiguousarray(arr[::-1, :]) + # Orientations 5-8 swap rows and columns. + if arr.ndim == 3: + # Transpose only the spatial axes; keep the sample axis trailing. + if orientation == 5: + return np.ascontiguousarray(arr.transpose(1, 0, 2)) + if orientation == 6: + return np.ascontiguousarray(arr.transpose(1, 0, 2)[:, ::-1]) + if orientation == 7: + return np.ascontiguousarray(arr.transpose(1, 0, 2)[::-1, ::-1]) + if orientation == 8: + return np.ascontiguousarray(arr.transpose(1, 0, 2)[::-1, :]) + else: + if orientation == 5: + return np.ascontiguousarray(arr.T) + if orientation == 6: + return np.ascontiguousarray(arr.T[:, ::-1]) + if orientation == 7: + return np.ascontiguousarray(arr.T[::-1, ::-1]) + if orientation == 8: + return np.ascontiguousarray(arr.T[::-1, :]) + raise ValueError( + f"Invalid TIFF Orientation tag value: {orientation} " + f"(must be 1-8 per TIFF 6.0)" + ) + + def read_to_array(source, *, window=None, overview_level: int | None = None, band: int | None = None, max_pixels: int = MAX_PIXELS_DEFAULT, @@ -1143,6 +1202,17 @@ def read_to_array(source, *, window=None, overview_level: int | None = None, dtype = tiff_dtype_to_numpy(bps, ifd.sample_format) geo_info = extract_geo_info(ifd, data, header.byte_order) + # Orientation tag (274): values 2-8 mean the stored pixel order + # differs from display order. We need to remap the array post + # decode. A windowed read against a non-default orientation has + # ambiguous semantics (does the window refer to file pixels or + # display pixels?) so we reject that combo rather than guess. + orientation = ifd.orientation + if orientation != 1 and window is not None: + raise ValueError( + "orientation != 1 with window= is not supported" + ) + if ifd.is_tiled: arr = _read_tiles(data, ifd, header, dtype, window, max_pixels=max_pixels) @@ -1150,6 +1220,26 @@ def read_to_array(source, *, window=None, overview_level: int | None = None, arr = _read_strips(data, ifd, header, dtype, window, max_pixels=max_pixels) + if orientation != 1: + arr = _apply_orientation(arr, orientation) + # Orientations 5-8 swap rows and columns, so the file's stored + # pixel_width sits on the y-axis of the displayed array and + # vice versa. Replace the transform with one whose pixel_width + # comes from the file's pixel_height (and vice versa) so that + # the returned coords have the right number of entries on each + # axis. Orientation 5-8 with rotated geographies is rare in + # practice; the sign convention here keeps y decreasing + # downward, matching the standard top-left-origin convention. + if orientation in (5, 6, 7, 8): + t = geo_info.transform + from ._geotags import GeoTransform + geo_info.transform = GeoTransform( + origin_x=t.origin_x, + origin_y=t.origin_y, + pixel_width=abs(t.pixel_height), + pixel_height=-abs(t.pixel_width), + ) + # For multi-band with band selection, extract single band if arr.ndim == 3 and ifd.samples_per_pixel > 1 and band is not None: arr = arr[:, :, band] diff --git a/xrspatial/geotiff/tests/test_orientation.py b/xrspatial/geotiff/tests/test_orientation.py new file mode 100644 index 00000000..46a859da --- /dev/null +++ b/xrspatial/geotiff/tests/test_orientation.py @@ -0,0 +1,149 @@ +"""TIFF Orientation tag (274) decode tests for issue #1503. + +Before the fix, the reader silently ignored tag 274 and returned the file's +stored pixel order regardless of which corner the data was supposed to +start from. Files written with orientation 2-8 decoded to flipped or +rotated arrays, a silent geometric error. + +Each test writes a small array with a specific orientation tag, then +checks that ``open_geotiff`` produces the array remapped per TIFF 6.0 +spec. tifffile's ``imwrite`` is used to embed the Orientation tag via +``extratags``; tifffile's ``imread`` does *not* itself apply orientation +(it returns the stored buffer as-is) so the comparison target here is +computed from the source array rather than read back through tifffile. +""" +from __future__ import annotations + +import numpy as np +import pytest + +from xrspatial.geotiff import open_geotiff +from xrspatial.geotiff._reader import read_to_array + +tifffile = pytest.importorskip("tifffile") + + +# Eight orientation values defined by TIFF 6.0. +_ORIENTATIONS = [1, 2, 3, 4, 5, 6, 7, 8] + + +def _write_with_orientation(path, arr, orientation): + """Write *arr* to *path* with the given Orientation tag value. + + tifffile's ``imwrite`` does not expose Orientation as a kwarg, but + the ``extratags`` parameter accepts (tag_id, dtype_code, count, value, + write_once) tuples that get emitted into the IFD verbatim. ``H`` is + the unsigned short (TIFF type 3) struct code. + """ + tifffile.imwrite( + str(path), + arr, + extratags=[(274, 'H', 1, orientation, True)], + ) + + +def _expected_for_orientation(stored, orientation): + """Return what *stored* should look like after applying *orientation*. + + Mirrors the spec table in :func:`xrspatial.geotiff._reader._apply_orientation`. + """ + 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 orientation == 6: + return stored.T[:, ::-1] + if orientation == 7: + return stored.T[::-1, ::-1] + if orientation == 8: + return stored.T[::-1, :] + raise AssertionError(orientation) + + +@pytest.mark.parametrize("orientation", _ORIENTATIONS) +def test_orientation_matches_spec(tmp_path, orientation): + """open_geotiff applies the spec-defined transform for each orientation.""" + # Asymmetric data (different height and width, distinct row/column + # values) so any axis swap or flip shows up as a clear mismatch. + arr = np.arange(24, dtype=np.uint8).reshape(4, 6) + path = tmp_path / f"orient_{orientation}.tif" + _write_with_orientation(path, arr, orientation) + + expected = _expected_for_orientation(arr, orientation) + got = open_geotiff(str(path)) + + assert got.values.shape == expected.shape, ( + f"orientation={orientation}: shape mismatch " + f"got={got.values.shape} expected={expected.shape}" + ) + np.testing.assert_array_equal(got.values, expected) + + +@pytest.mark.parametrize("orientation", _ORIENTATIONS) +def test_orientation_coords_match_post_orientation_shape( + tmp_path, orientation +): + """y/x coordinate arrays size matches the post-orientation array shape.""" + arr = np.arange(24, dtype=np.uint8).reshape(4, 6) + path = tmp_path / f"orient_coords_{orientation}.tif" + _write_with_orientation(path, arr, orientation) + + da = open_geotiff(str(path)) + + h, w = da.values.shape + assert da.coords['y'].shape == (h,) + assert da.coords['x'].shape == (w,) + + +@pytest.mark.parametrize("orientation", [5, 6, 7, 8]) +def test_orientation_5_to_8_swap_dims(tmp_path, orientation): + """Orientations 5-8 swap rows and columns relative to the stored shape.""" + arr = np.arange(24, dtype=np.uint8).reshape(4, 6) # h=4, w=6 + path = tmp_path / f"orient_swap_{orientation}.tif" + _write_with_orientation(path, arr, orientation) + + da = open_geotiff(str(path)) + + # File stores h=4, w=6. After orientation 5-8 the displayed shape is + # (6, 4) -- width and height swap. + assert da.values.shape == (6, 4) + + +def test_orientation_default_unchanged(tmp_path): + """A file without an Orientation tag defaults to 1 (no transform).""" + arr = np.arange(24, dtype=np.uint8).reshape(4, 6) + path = tmp_path / "no_orient.tif" + tifffile.imwrite(str(path), arr) + + da = open_geotiff(str(path)) + np.testing.assert_array_equal(da.values, arr) + + +def test_orientation_with_window_raises(tmp_path): + """Windowed read on a non-default orientation raises ValueError.""" + arr = np.arange(24, dtype=np.uint8).reshape(4, 6) + path = tmp_path / "orient2_window.tif" + _write_with_orientation(path, arr, 2) + + with pytest.raises(ValueError, match="orientation"): + read_to_array(str(path), window=(0, 0, 2, 2)) + + with pytest.raises(ValueError, match="orientation"): + open_geotiff(str(path), window=(0, 0, 2, 2)) + + +def test_orientation_1_with_window_still_works(tmp_path): + """Default orientation (1) with window= keeps working as before.""" + arr = np.arange(24, dtype=np.uint8).reshape(4, 6) + path = tmp_path / "orient1_window.tif" + _write_with_orientation(path, arr, 1) + + da = open_geotiff(str(path), window=(0, 0, 2, 3)) + assert da.values.shape == (2, 3) + np.testing.assert_array_equal(da.values, arr[:2, :3])