diff --git a/xrspatial/geotiff/_geotags.py b/xrspatial/geotiff/_geotags.py index fe662403..998f5f14 100644 --- a/xrspatial/geotiff/_geotags.py +++ b/xrspatial/geotiff/_geotags.py @@ -8,7 +8,7 @@ IFD, TAG_IMAGE_WIDTH, TAG_IMAGE_LENGTH, TAG_BITS_PER_SAMPLE, TAG_COMPRESSION, TAG_PHOTOMETRIC, - TAG_STRIP_OFFSETS, TAG_SAMPLES_PER_PIXEL, + TAG_STRIP_OFFSETS, TAG_ORIENTATION, TAG_SAMPLES_PER_PIXEL, TAG_ROWS_PER_STRIP, TAG_STRIP_BYTE_COUNTS, TAG_X_RESOLUTION, TAG_Y_RESOLUTION, TAG_PLANAR_CONFIG, TAG_RESOLUTION_UNIT, @@ -34,7 +34,7 @@ _MANAGED_TAGS = frozenset({ TAG_IMAGE_WIDTH, TAG_IMAGE_LENGTH, TAG_BITS_PER_SAMPLE, TAG_COMPRESSION, TAG_PHOTOMETRIC, - TAG_STRIP_OFFSETS, TAG_SAMPLES_PER_PIXEL, + TAG_STRIP_OFFSETS, TAG_ORIENTATION, TAG_SAMPLES_PER_PIXEL, TAG_ROWS_PER_STRIP, TAG_STRIP_BYTE_COUNTS, TAG_X_RESOLUTION, TAG_Y_RESOLUTION, TAG_PLANAR_CONFIG, TAG_RESOLUTION_UNIT, 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..fed4338b 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,24 @@ 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. + # ``read_geotiff_dask`` chunks the file by issuing windowed reads, + # so this check also rejects ``chunks=`` for non-default + # orientation; the error mentions both so the failure is easy to + # diagnose if it surfaces under dask. + orientation = ifd.orientation + 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." + ) + if ifd.is_tiled: arr = _read_tiles(data, ifd, header, dtype, window, max_pixels=max_pixels) @@ -1150,10 +1227,47 @@ def read_to_array(source, *, window=None, overview_level: int | None = None, arr = _read_strips(data, ifd, header, dtype, window, max_pixels=max_pixels) - # For multi-band with band selection, extract single band + # Extract the requested band before reorienting so we work on a + # smaller 2D array rather than reorienting a full multi-band cube + # only to slice it afterwards. if arr.ndim == 3 and ifd.samples_per_pixel > 1 and band is not None: arr = arr[:, :, band] + 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. Swap the transform's pixel sizes so the coord + # arrays come out the right length. Signs are preserved + # rather than coerced to north-up, since some legitimate files + # use a non-standard sign convention (south-up, west-up). + # + # For orientations 6/7/8 (rotations + flips, not a pure + # transpose) the swap is geometrically inexact for georef'd + # files: a strict implementation would also adjust origin + # and re-sign per axis. Such files are vanishingly rare in + # practice (TIFF Orientation 5-8 with a meaningful + # ModelTransformation), and getting it right requires a + # design pass; we warn instead so the user knows to verify. + if orientation in (5, 6, 7, 8): + t = geo_info.transform + 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): + import warnings + 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=2, + ) + # MinIsWhite (photometric=0): invert single-band grayscale values if ifd.photometric == 0 and ifd.samples_per_pixel == 1: if arr.dtype.kind == 'u': diff --git a/xrspatial/geotiff/tests/test_orientation.py b/xrspatial/geotiff/tests/test_orientation.py new file mode 100644 index 00000000..da693002 --- /dev/null +++ b/xrspatial/geotiff/tests/test_orientation.py @@ -0,0 +1,257 @@ +"""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]) + + +@pytest.mark.parametrize("orientation", [2, 3, 4, 5, 6, 7, 8]) +def test_orientation_tag_not_passed_through_extra_tags(tmp_path, orientation): + """Tag 274 must not survive on the returned DataArray. + + Without this, a read+write round-trip would re-emit the original + Orientation tag even though the pixel buffer is already remapped -- + downstream readers would apply the orientation a second time. + """ + arr = np.arange(24, dtype=np.uint8).reshape(4, 6) + path = tmp_path / f"orient_passthrough_{orientation}.tif" + _write_with_orientation(path, arr, orientation) + + da = open_geotiff(str(path)) + extra = da.attrs.get('extra_tags') or [] + tag_ids = [t[0] if isinstance(t, tuple) else t for t in extra] + assert 274 not in tag_ids, ( + f"orientation={orientation}: tag 274 leaked into extra_tags={extra}" + ) + + +def test_orientation_round_trip_does_not_double_apply(tmp_path): + """open_geotiff -> to_geotiff -> open_geotiff returns the same array. + + Concretely: a file written with orientation=4 reads as flipped + (correct), the writer emits a normal file (no orientation tag), and + a second read returns the same array. If tag 274 leaked through + extra_tags, the second read would apply orientation=4 again. + """ + from xrspatial.geotiff import to_geotiff + + arr = np.arange(24, dtype=np.uint8).reshape(4, 6) + path1 = tmp_path / "orient4_in.tif" + _write_with_orientation(path1, arr, 4) + + da1 = open_geotiff(str(path1)) + + path2 = tmp_path / "orient4_out.tif" + to_geotiff(da1, str(path2)) + da2 = open_geotiff(str(path2)) + + np.testing.assert_array_equal(da2.values, da1.values) + + +@pytest.mark.parametrize("orientation", [5, 6, 7, 8]) +def test_orientation_5_to_8_warn_on_georef(tmp_path, orientation): + """Axis-swap orientations on georef'd files emit a warning. + + The transform swap is shape-correct but geometrically approximate + for orientations 6/7/8; the user should be told. + """ + import warnings as _warnings + + arr = np.arange(24, dtype=np.uint8).reshape(4, 6) + path = tmp_path / f"orient_georef_{orientation}.tif" + # tifffile lets us tag a CRS via the resolution + metadata; the + # simplest path that triggers our georef branch is to embed a + # ModelPixelScale + GeoKeyDirectory pair pointing at EPSG:4326. + tifffile.imwrite( + str(path), arr, + extratags=[ + (274, 'H', 1, orientation, True), + # ModelPixelScaleTag (33550): scale_x, scale_y, scale_z + (33550, 'd', 3, (1.0, 1.0, 0.0), True), + # ModelTiepointTag (33922): I, J, K, X, Y, Z + (33922, 'd', 6, (0.0, 0.0, 0.0, 0.0, 0.0, 0.0), True), + # GeoKeyDirectoryTag (34735): version 1.1.0, 1 key, + # GTModelType=2 (Geographic), GeographicType=4326 + (34735, 'H', 12, ( + 1, 1, 0, 2, + 1024, 0, 1, 2, # GTModelTypeGeoKey = 2 + 2048, 0, 1, 4326, # GeographicTypeGeoKey = 4326 + ), True), + ], + ) + + with _warnings.catch_warnings(record=True) as caught: + _warnings.simplefilter('always') + da = open_geotiff(str(path)) + + assert da is not None + msgs = [str(w.message) for w in caught if 'Orientation' in str(w.message)] + assert msgs, ( + f"orientation={orientation}: no warning emitted on georef file" + ) + + +def test_orientation_with_band_selection_returns_2d(tmp_path): + """band= followed by an orientation transpose returns a 2D array. + + Regression: an earlier draft applied orientation before slicing the + band, which wasted memory and produced confusing intermediates. + Now band slicing happens first. + """ + rgb = np.arange(4 * 6 * 3, dtype=np.uint8).reshape(4, 6, 3) + path = tmp_path / "orient5_rgb.tif" + tifffile.imwrite( + str(path), rgb, photometric='rgb', + extratags=[(274, 'H', 1, 5, True)], + ) + + # Orientation 5 transposes spatial axes, so output spatial shape is + # (6, 4). Band 1 returns just that channel. + da = open_geotiff(str(path), band=1) + assert da.values.shape == (6, 4) + expected = rgb[..., 1].T # band 1 then transpose + np.testing.assert_array_equal(da.values, expected)