diff --git a/docs/source/user_guide/attrs_contract.rst b/docs/source/user_guide/attrs_contract.rst index 265ecd8b0..26b2e1df3 100644 --- a/docs/source/user_guide/attrs_contract.rst +++ b/docs/source/user_guide/attrs_contract.rst @@ -63,6 +63,19 @@ write. affine transform tuple matching the rasterio ``Affine`` ordering. Omitted for files with no ``ModelTransformation`` / ``ModelPixelScale`` / ``ModelTiepoint`` tags. + * - ``rotated_affine`` + - tuple + - Full ``(a, b, c, d, e, f)`` rasterio-style 6-tuple for files + opened with ``allow_rotated=True`` whose source carried a + rotated / sheared ``ModelTransformationTag``. The axis-aligned + ``transform`` cannot express the rotation terms; this attr + surfaces the rotated mapping so downstream code (custom warps, + visualisation) can recover it. Only emitted on the rotated + opt-in path; absent on axis-aligned reads and on plain + no-georef files. Read-only -- ``to_geotiff`` drops the rotation + on the way out until the writer learns to emit + ``ModelTransformationTag`` (issue #2115 follow-up). See issue + #2129. * - ``nodata`` - scalar - Numeric NoData sentinel. Emitted by readers when the file @@ -124,7 +137,7 @@ write. ``ResolutionUnit`` ids 1, 2, 3). * - ``_xrspatial_geotiff_contract`` - int - - Contract version. Currently ``2``. See `Versioning`_. + - Contract version. Currently ``4``. See `Versioning`_. * - ``_xrspatial_no_georef`` - bool - Stamped ``True`` on reads of files with no GeoTIFF transform @@ -328,7 +341,7 @@ Versioning ========== The contract is versioned through ``attrs['_xrspatial_geotiff_contract']``. -The current value is ``2``. Future revisions that add canonical keys, +The current value is ``4``. Future revisions that add canonical keys, move keys between tiers, or change a key's semantics will bump the integer. Callers that depend on a specific layout can branch on the version, and writers will emit the version they were built against. @@ -344,3 +357,17 @@ and matplotlib-colormap attrs that v1 emitted on read under a ``attrs[key]`` will now see ``KeyError``; switch to ``attrs.get(key)`` or migrate to the canonical ``crs`` / ``crs_wkt`` plus :mod:`pyproj` recipe documented in `Removed in contract v2`_. + +Contract v3 (issue #2136) added the ``georef_status`` attr to the +canonical tier, encoding the five distinct states the reader can +land in (``full``, ``transform_only``, ``crs_only``, ``none``, +``rotated_dropped``) so downstream code can branch on a single +value instead of reconstructing the state from the union of ``crs``, +``crs_wkt``, ``transform``, and ``_xrspatial_no_georef``. + +Contract v4 (issue #2129) added the ``rotated_affine`` attr to the +canonical tier. The attr surfaces the rotated 6-tuple from +``ModelTransformationTag`` on the ``allow_rotated=True`` opt-in path +so callers can recover the rotated mapping. The writer drops it on +round-trip until ``to_geotiff`` learns to emit +``ModelTransformationTag`` (issue #2115 follow-up). diff --git a/xrspatial/geotiff/__init__.py b/xrspatial/geotiff/__init__.py index 63e4561f1..4009ad275 100644 --- a/xrspatial/geotiff/__init__.py +++ b/xrspatial/geotiff/__init__.py @@ -464,10 +464,14 @@ def open_geotiff(source: str | BinaryIO, *, transform because keeping them while the axis-aligned transform is gone misleads downstream code that gates on ``"crs" in da.attrs`` to mean the array is spatially usable - (issue #2126). The contract is read-only -- ``to_geotiff`` does - not currently emit rotated transforms, so a read-then-write - round-trip writes an identity-affine output and silently drops - the rotation (issue #2115). + (issue #2126). The rotated 6-tuple itself is surfaced on + ``attrs['rotated_affine']`` as ``(a, b, c, d, e, f)`` (rasterio + ``Affine`` ordering) so consumers that know how to handle + rotated rasters can recover the mapping (issue #2129). The + contract is read-only -- ``to_geotiff`` does not currently + emit rotated transforms, so a read-then-write round-trip + writes an identity-affine output and silently drops the + rotation (issue #2115). Returns ------- diff --git a/xrspatial/geotiff/_attrs.py b/xrspatial/geotiff/_attrs.py index 17cc4dc7f..6a9d3f63d 100644 --- a/xrspatial/geotiff/_attrs.py +++ b/xrspatial/geotiff/_attrs.py @@ -42,6 +42,14 @@ ModelPixelScale, or ModelTiepoint), and for rotated reads opened with ``allow_rotated=True`` (axis-aligned 6-tuple would silently drop the rotation terms). +- ``rotated_affine`` (#2129): rasterio-style 6-tuple + ``(a, b, c, d, e, f)`` capturing the full ``ModelTransformationTag`` + on the ``allow_rotated=True`` opt-in path. Only emitted when the + source carried a rotated / sheared transform; absent on plain + no-georef reads and on axis-aligned reads (which already round-trip + via ``transform``). Read-only -- the writer drops it on the way out + until ``to_geotiff`` learns to emit ``ModelTransformationTag`` + (issue #2115 follow-up). - ``nodata``: declared file sentinel as stored in the GDAL_NODATA tag. Set whenever the source declares one, as a scalar of the source dtype, regardless of whether the in-memory array is float-with-NaN @@ -258,7 +266,13 @@ # code that branches on them still works; the new attr is additive and # disambiguates ``crs_only`` from ``none`` and ``rotated_dropped`` from # the truly-no-transform case. -_ATTRS_CONTRACT_VERSION = 3 +# +# Version 4 (issue #2129) adds ``attrs['rotated_affine']`` for the +# ``allow_rotated=True`` opt-in path. The 6-tuple is read-only -- the +# writer drops it on round-trip until ``to_geotiff`` grows a +# ``ModelTransformationTag`` emit path (#2115 follow-up). Existing keys +# keep their pre-v4 shape. +_ATTRS_CONTRACT_VERSION = 4 # Canonical ``attrs['georef_status']`` values (issue #2136). One attr @@ -347,6 +361,16 @@ class GeoTIFFMetadata: # branching on attrs after the dict has been built. georef_status: str | None = None + # Rotated 6-tuple from ``ModelTransformationTag`` on the + # ``allow_rotated=True`` opt-in path (issue #2129). Carried on the + # record so the eager / dask / GPU / VRT read paths emit + # ``attrs['rotated_affine']`` through the same marshalling step. + # Read-only: :func:`attrs_to_metadata` intentionally does NOT + # populate this field from incoming attrs so the writer keeps + # dropping the rotation on round-trip until ``to_geotiff`` learns to + # emit ``ModelTransformationTag`` (#2115 follow-up). + rotated_affine: tuple | None = None + # Contract version stamped on read contract_version: int = _ATTRS_CONTRACT_VERSION @@ -399,6 +423,15 @@ def geo_info_to_metadata(geo_info, *, window=None) -> GeoTIFFMetadata: crs_epsg = None if rotated_optin else geo_info.crs_epsg crs_wkt = None if rotated_optin else geo_info.crs_wkt + # Surface the rotated 6-tuple on the public attrs (issue #2129) so + # downstream code that knows how to handle rotated rasters can read + # it without diving into the internal ``GeoInfo`` / ``GeoTransform`` + # objects. Tuple cast normalises lists or numpy sequences coming + # from the parser into the documented ``tuple`` shape. + rotated_affine_tuple = ( + tuple(src_t.rotated_affine) if rotated_optin else None + ) + raster_type = ( 'point' if geo_info.raster_type == RASTER_PIXEL_IS_POINT else 'area') @@ -436,6 +469,7 @@ def geo_info_to_metadata(geo_info, *, window=None) -> GeoTIFFMetadata: # ``_compute_georef_status_from_parts`` to fill this field # without synthesising a ``GeoInfo``. georef_status=_compute_georef_status(geo_info), + rotated_affine=rotated_affine_tuple, contract_version=_ATTRS_CONTRACT_VERSION, ) @@ -482,6 +516,15 @@ def metadata_to_attrs(md: GeoTIFFMetadata) -> dict: elif not md.has_georef: attrs[_NO_GEOREF_KEY] = True + # ``rotated_affine`` (issue #2129) rides alongside the + # ``_xrspatial_no_georef`` marker on the ``allow_rotated=True`` path + # so callers can recover the rotated mapping. Only set on read; the + # writer-side :func:`attrs_to_metadata` deliberately does not parse + # it back, so a read-then-write round-trip drops the rotation until + # the writer grows ``ModelTransformationTag`` emit support (#2115). + if md.rotated_affine is not None: + attrs['rotated_affine'] = md.rotated_affine + if md.nodata is not None: attrs['nodata'] = md.nodata attrs['masked_nodata'] = bool(md.masked_nodata) diff --git a/xrspatial/geotiff/_backends/vrt.py b/xrspatial/geotiff/_backends/vrt.py index 58828cf70..aa71776c5 100644 --- a/xrspatial/geotiff/_backends/vrt.py +++ b/xrspatial/geotiff/_backends/vrt.py @@ -333,6 +333,16 @@ def read_vrt(source: str, *, # ``rotated_dropped`` bucket as a rotated ``ModelTransformationTag``. _vrt_keep_crs = bool(vrt.crs_wkt) and not _vrt_is_rotated _vrt_epsg = _wkt_to_epsg(vrt.crs_wkt) if _vrt_keep_crs else None + # Surface the rotated 6-tuple alongside ``georef_status='rotated_dropped'`` + # so the VRT path matches the non-VRT ``ModelTransformationTag`` path + # introduced by issue #2129. Without this the rotated VRT would land + # in the same bucket as the rotated TIFF but offer no way to recover + # the mapping. The GDAL geo_transform is already 6-tuple ordered; + # ``_gdal_geotransform_to_affine_tuple`` converts to rasterio + # ``Affine`` ordering (a, b, c, d, e, f). + _vrt_rotated_affine = ( + _gdal_geotransform_to_affine_tuple(gt) if _vrt_is_rotated else None + ) _vrt_md = GeoTIFFMetadata( crs_epsg=_vrt_epsg, crs_wkt=vrt.crs_wkt if _vrt_keep_crs else None, @@ -347,6 +357,7 @@ def read_vrt(source: str, *, has_crs=vrt.crs_wkt is not None and not _vrt_is_rotated, rotated_dropped=_vrt_is_rotated, ), + rotated_affine=_vrt_rotated_affine, ) attrs = metadata_to_attrs(_vrt_md) # When a specific band is selected, source its nodata from that @@ -836,6 +847,12 @@ def _read_vrt_chunked(source, *, window, band, name, chunks, gpu, dtype, # Rotated VRTs drop CRS attrs alongside the transform (#2122). _vrt_keep_crs = bool(vrt.crs_wkt) and not _vrt_is_rotated _vrt_epsg = _wkt_to_epsg(vrt.crs_wkt) if _vrt_keep_crs else None + # See the eager VRT branch for the rationale; issue #2129 carries + # the rotated 6-tuple onto the chunked VRT path too so dask reads + # of rotated VRTs match the eager-VRT and non-VRT TIFF surface. + _vrt_rotated_affine = ( + _gdal_geotransform_to_affine_tuple(gt) if _vrt_is_rotated else None + ) # ``georef_status`` (issue #2136). See the eager VRT branch above # for the rationale; the rotated VRT path lands the array in the # ``rotated_dropped`` bucket so consumers can branch on it. @@ -850,6 +867,7 @@ def _read_vrt_chunked(source, *, window, band, name, chunks, gpu, dtype, has_crs=vrt.crs_wkt is not None and not _vrt_is_rotated, rotated_dropped=_vrt_is_rotated, ), + rotated_affine=_vrt_rotated_affine, ) attrs = metadata_to_attrs(_vrt_md) diff --git a/xrspatial/geotiff/tests/test_attrs_contract_version_1984.py b/xrspatial/geotiff/tests/test_attrs_contract_version_1984.py index 90bd53102..c2976294d 100644 --- a/xrspatial/geotiff/tests/test_attrs_contract_version_1984.py +++ b/xrspatial/geotiff/tests/test_attrs_contract_version_1984.py @@ -82,10 +82,12 @@ def test_attrs_contract_version_constant_is_current(): """Pin the integer value so a careless bump shows up here first. Contract v3 (issue #2136) added ``attrs['georef_status']`` to the - canonical tier. Bumping past 3 should be paired with a docs update - and a sibling test for the new key. + canonical tier. Contract v4 (issue #2129) added + ``attrs['rotated_affine']`` for the ``allow_rotated=True`` opt-in + path. Bumping past 4 should be paired with a docs update and a + sibling test for the new key. """ - assert _ATTRS_CONTRACT_VERSION == 3 + assert _ATTRS_CONTRACT_VERSION == 4 def test_eager_numpy_stamps_contract_version(tmp_path): diff --git a/xrspatial/geotiff/tests/test_georef_status_2136.py b/xrspatial/geotiff/tests/test_georef_status_2136.py index 311cba08d..e0b6e0561 100644 --- a/xrspatial/geotiff/tests/test_georef_status_2136.py +++ b/xrspatial/geotiff/tests/test_georef_status_2136.py @@ -75,9 +75,11 @@ def _gpu_available() -> bool: # --------------------------------------------------------------------------- -def test_contract_version_is_three(): - """The new attr lands in v3; pin the bump alongside the new key.""" - assert _ATTRS_CONTRACT_VERSION == 3 +def test_contract_version_is_at_least_three(): + """The ``georef_status`` attr lands in v3; pin a lower bound so future + contract bumps that keep the attr (e.g. ``rotated_affine`` in v4 / + issue #2129) do not regress this test.""" + assert _ATTRS_CONTRACT_VERSION >= 3 def test_public_constants_reexported(): diff --git a/xrspatial/geotiff/tests/test_geotiff_metadata_2139.py b/xrspatial/geotiff/tests/test_geotiff_metadata_2139.py index b64e3fb8e..63bc7bd24 100644 --- a/xrspatial/geotiff/tests/test_geotiff_metadata_2139.py +++ b/xrspatial/geotiff/tests/test_geotiff_metadata_2139.py @@ -244,6 +244,14 @@ def _representative_attrs_dicts(): * no-georef file * user-defined CRS (WKT only, no EPSG) * VRT with holes + + The ``rotated_affine`` attr (issue #2129) is deliberately excluded: + it is emitted on read but :func:`attrs_to_metadata` does NOT parse + it back (the writer drops the rotation until #2115 ships + ``ModelTransformationTag`` support). Adding it here would assert a + symmetry the contract intentionally breaks. See + ``test_rotated_affine_attr_2129.py::test_attrs_to_metadata_drops_rotated_affine`` + for the pin on the read-only direction. """ yield { '_xrspatial_geotiff_contract': _ATTRS_CONTRACT_VERSION, diff --git a/xrspatial/geotiff/tests/test_rotated_affine_attr_2129.py b/xrspatial/geotiff/tests/test_rotated_affine_attr_2129.py new file mode 100644 index 000000000..c6ee4dc20 --- /dev/null +++ b/xrspatial/geotiff/tests/test_rotated_affine_attr_2129.py @@ -0,0 +1,337 @@ +"""``allow_rotated=True`` surfaces ``attrs['rotated_affine']`` (#2129). + +Issue #2126 fixed the CRS-drop side of the ``allow_rotated=True`` +contract, but the rotated 6-tuple itself was still unreachable from +public code: it lived on ``geo_info.transform.rotated_affine`` while +``geo_info`` was not on the returned DataArray. This file pins the +follow-up that emits ``attrs['rotated_affine']`` so callers can read +the rotated mapping without reaching into reader internals. + +The attr: + +* appears only when the source carried a rotated ``ModelTransformationTag`` + AND the caller passed ``allow_rotated=True``; +* is a rasterio-style 6-tuple ``(a, b, c, d, e, f)`` matching the + ordering already documented on ``GeoTransform.rotated_affine``; +* is dropped on the way back through ``to_geotiff`` until the writer + grows ``ModelTransformationTag`` support (#2115 follow-up). The + drop is verified by feeding a synthesised attrs dict through + ``attrs_to_metadata`` and asserting the parsed record does not + carry the rotated 6-tuple forward. +""" +from __future__ import annotations + +import numpy as np +import pytest + +from xrspatial.geotiff import open_geotiff +from xrspatial.geotiff._attrs import ( + _ATTRS_CONTRACT_VERSION, + _populate_attrs_from_geo_info, + attrs_to_metadata, + geo_info_to_metadata, +) +from xrspatial.geotiff._geotags import GeoInfo, GeoTransform + + +_ROTATED_TUPLE = (8.66, -5.0, 100.0, 5.0, 8.66, 200.0) + + +def _rotated_geo_info(*, with_crs: bool = True) -> GeoInfo: + """Build a GeoInfo that mimics the ``allow_rotated=True`` parser output.""" + t = GeoTransform( + origin_x=0.0, origin_y=0.0, + pixel_width=1.0, pixel_height=-1.0, + rotated_affine=_ROTATED_TUPLE, + ) + return GeoInfo( + transform=t, + has_georef=False, + crs_epsg=4326 if with_crs else None, + crs_wkt='GEOGCS["WGS 84"]' if with_crs else None, + ) + + +# --------------------------------------------------------------------------- +# Unit tests on the helpers (no file I/O) +# --------------------------------------------------------------------------- + + +def test_rotated_optin_emits_rotated_affine_tuple(): + gi = _rotated_geo_info() + attrs: dict = {} + _populate_attrs_from_geo_info(attrs, gi) + + assert attrs.get('rotated_affine') == _ROTATED_TUPLE + # Sanity: rotated path still drops crs / transform (existing #2126 + # contract). Re-checked here so a regression to either branch + # surfaces in the same test file as the new attr. + assert 'crs' not in attrs + assert 'transform' not in attrs + + +def test_rotated_optin_emits_rotated_affine_without_crs(): + gi = _rotated_geo_info(with_crs=False) + attrs: dict = {} + _populate_attrs_from_geo_info(attrs, gi) + + assert attrs.get('rotated_affine') == _ROTATED_TUPLE + + +def test_plain_no_georef_omits_rotated_affine(): + # Plain no-georef file: no transform tag, no rotated matrix. The new + # attr must NOT appear on this path -- it is reserved for the + # ``allow_rotated=True`` opt-in. + gi = GeoInfo(transform=GeoTransform(), has_georef=False, crs_epsg=4326) + attrs: dict = {} + _populate_attrs_from_geo_info(attrs, gi) + + assert 'rotated_affine' not in attrs + + +def test_axis_aligned_read_omits_rotated_affine(): + gi = GeoInfo( + transform=GeoTransform( + origin_x=100.0, origin_y=200.0, + pixel_width=10.0, pixel_height=-10.0, + ), + has_georef=True, + crs_epsg=4326, + ) + attrs: dict = {} + _populate_attrs_from_geo_info(attrs, gi) + + assert 'rotated_affine' not in attrs + + +def test_rotated_affine_is_tuple_not_list(): + # ``GeoTransform.rotated_affine`` is documented as a tuple; the + # public attr should respect the same type so downstream code can + # rely on ``isinstance(attrs['rotated_affine'], tuple)``. The + # parser tuple-casts the field even if a future change stores it + # as a list or numpy sequence. + gi = _rotated_geo_info() + # Replace with a list to simulate a parser change. + gi.transform.rotated_affine = list(_ROTATED_TUPLE) + md = geo_info_to_metadata(gi) + + assert isinstance(md.rotated_affine, tuple) + assert md.rotated_affine == _ROTATED_TUPLE + + +# --------------------------------------------------------------------------- +# Round-trip contract: the writer parser must drop ``rotated_affine``. +# --------------------------------------------------------------------------- + + +def test_attrs_to_metadata_drops_rotated_affine(): + """The write-side boundary parser intentionally does not carry the + rotated 6-tuple forward; the writer would otherwise need a + ``ModelTransformationTag`` emit path (#2115). Keeping it off the + record ensures ``to_geotiff`` keeps writing a plain no-georef file + until that follow-up lands.""" + attrs = { + 'rotated_affine': _ROTATED_TUPLE, + '_xrspatial_no_georef': True, + '_xrspatial_geotiff_contract': _ATTRS_CONTRACT_VERSION, + } + md = attrs_to_metadata(attrs) + + assert md.rotated_affine is None + assert md.has_georef is False + assert md.transform is None + + +# --------------------------------------------------------------------------- +# End-to-end via open_geotiff against a synthesised rotated GeoTIFF. +# --------------------------------------------------------------------------- + + +def _write_rotated_tiff(path, arr, *, epsg=None): + """Write a rotated GeoTIFF with an optional GeoKey-encoded CRS. + + The 4x4 ``ModelTransformationTag`` matrix encodes a 30-degree + rotation with 10-unit pixel spacing. Mirrors the helper in + ``test_allow_rotated_crs_drop_2126.py`` so the two suites share the + same synthetic file shape. + """ + tifffile = pytest.importorskip("tifffile") + cos30 = 0.8660254037844387 + sin30 = 0.5 + m = ( + 10.0 * cos30, -10.0 * sin30, 0.0, 100.0, + 10.0 * sin30, 10.0 * cos30, 0.0, 200.0, + 0.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, 1.0, + ) + extratags = [ + # ModelTransformationTag (34264) -- DOUBLE, count=16. + (34264, 12, 16, m, False), + ] + if epsg is not None: + # GeoKeyDirectory: header (4 shorts) + one key for GeographicTypeGeoKey. + geo_key_directory = ( + 1, 1, 0, 1, + 2048, 0, 1, int(epsg), + ) + extratags.append((34735, 3, 8, geo_key_directory, False)) + tifffile.imwrite( + path, arr, photometric='minisblack', + planarconfig='contig', extratags=extratags, + ) + return m + + +def test_open_geotiff_rotated_emits_rotated_affine(tmp_path): + arr = np.arange(20, dtype='``. + + GDAL ``geo_transform`` ordering is + ``(origin_x, pixel_width, rot_x, origin_y, rot_y, pixel_height)``; + non-zero ``rot_x`` / ``rot_y`` is what the reader treats as a + rotated VRT (see ``_vrt_is_rotated`` in ``_backends/vrt.py``). + """ + vrt_xml = ( + f'\n' + f' {gt_str}\n' + f' \n' + f' \n' + f' {src_path}\n' + f' 1\n' + f' \n' + f' \n' + f' \n' + f' \n' + f'\n' + ) + vrt = tmp_path / "tmp_2129_rotated.vrt" + vrt.write_text(vrt_xml) + return str(vrt) + + +def test_open_geotiff_rotated_vrt_emits_rotated_affine(tmp_path): + """A VRT with a rotated ```` opened with + ``allow_rotated=True`` lands in ``georef_status='rotated_dropped'`` + and must surface ``rotated_affine`` so callers can recover the + mapping. Mirrors the non-VRT ``ModelTransformationTag`` path.""" + tifffile = pytest.importorskip("tifffile") + arr = np.arange(16, dtype='