geotiff: synthesize WKT for citation-only user-defined CRSes (#1930) (#2054)

* geotiff: synthesize WKT for citation-only user-defined CRSes (#1930)

Closes the citation-only CRS parity gap that every phase-3 backend
module flagged. The crs_citation_only fixture in the golden corpus
carries a user-defined geographic CRS (no EPSG, no WKT in the
citation) and the reader was decoding it only into the deprecated
geog_citation / datum_code / angular_units attrs. attrs['crs_wkt']
stayed None, so the oracle's _candidate_crs returned None on the
candidate side and compare_to_oracle failed.

This change adds _synthesize_user_defined_wkt to _geotags.py. When
the file declares a user-defined geographic CRS and exposes the
ellipsoid (semi_major plus semi_minor or inv_flattening) and the
angular-units GeoKeys, the helper feeds those parameters to
pyproj.CRS.from_dict and stamps the resulting WKT on
attrs['crs_wkt']. The synthesis fails closed (returns None) when
pyproj is missing, when there is no semi_major, or when the
ellipsoid shape is ambiguous, so the existing deprecated-attrs path
stays in place.

Projected user-defined CRSes (ModelTypeProjected with
ProjectedCSType == 32767) are not yet reconstructible from GeoKeys
alone -- they need the GeogPrime / Projection parameters that the
corpus does not exercise -- so the helper returns None for that
case.

Drops the crs_citation_only entry from _PARITY_GAPS in each of the
four phase-3 backend modules (eager numpy, dask numpy, GPU,
dask+GPU). xfail(strict=True) flips to a real failure when the test
starts passing, so the fix and the entry cleanups have to land
together. Each module's top-of-file docstring moves the citation
entry from "Real parity gaps" to "Resolved gaps", mirroring the
masked-nodata fix from PR #2046.

Adds unit tests for _synthesize_user_defined_wkt in
test_user_defined_crs_wkt_1632.py: the sphere case (the citation
fixture shape), the oblate ellipsoid case, the projected-CRS skip
path, and the missing-ellipsoid skip path. Adds end-to-end fixture
tests in test_oracle.py for attrs['crs_wkt'] stamping, PROJ-dict
equality with the rasterio reference, and compare_to_oracle
round-trip.

After the fix the four phase-3 backend modules report 115 passed,
8 skipped, 4 xfailed (the remaining JPEG axis-order gap). The full
golden_corpus / corpus suite reports 253 passed, 8 skipped, 4
xfailed with no XPASS.

* Address review nits: named constant, drop unused params (#1930)

* Adds ``GEOKEY_GEOG_SEMI_MINOR_AXIS = 2058`` next to the other
  GeoKey ID constants and uses it in ``extract_geo_info`` instead of
  the raw int literal. Matches the existing convention for the
  semi-major-axis and inv-flattening keys.

* Drops ``angular_units_code`` and ``geog_citation`` from
  ``_synthesize_user_defined_wkt``. The function never read them.
  PROJ's ``longlat`` always emits degrees, and the corpus has no
  radian-unit user-defined fixture; the docstring now calls that
  limit out so a future radian fixture forces a deliberate signature
  change rather than silently degrading the units.

* Expands the docstring's "not handled" section to call out
  ``MODEL_TYPE_GEOCENTRIC`` (3) and unknown / zero model_type. Adds
  a ``test_synthesize_user_defined_wkt_geocentric_returns_none``
  pin so a future change that promotes geocentric to a real
  proj_dict has to touch the test deliberately.

Full corpus suite still reports 287 passed, 8 skipped, 4 xfailed
(JPEG axis-order gap), 0 XPASS.

Author: brendancol

xrspatial/geotiff/_geotags.py

@@ -67,6 +67,7 @@
 GEOKEY_GEOG_LINEAR_UNITS = 2052
 GEOKEY_GEOG_ANGULAR_UNITS = 2054
 GEOKEY_GEOG_SEMI_MAJOR_AXIS = 2057
+GEOKEY_GEOG_SEMI_MINOR_AXIS = 2058
 GEOKEY_GEOG_INV_FLATTENING = 2059
 GEOKEY_PROJECTED_CS_TYPE = 3072
 GEOKEY_PROJ_CITATION = 3073
@@ -308,6 +309,92 @@ def _looks_like_wkt(text) -> bool:
     return any(head.startswith(p) for p in _WKT_PREFIXES)
 
 
+def _synthesize_user_defined_wkt(
+    model_type: int,
+    semi_major: float | None,
+    semi_minor: float | None,
+    inv_flattening: float | None,
+) -> str | None:
+    """Build a canonical WKT for a user-defined geographic CRS.
+
+    Triggered from :func:`extract_geo_info` when the file declares a
+    user-defined geographic CRS (``GeographicTypeGeoKey == 32767`` and
+    no projected EPSG) and the citation does not already carry a WKT
+    string. The classic example is a GeoTIFF whose CRS lives only in
+    ``GeogCitationGeoKey`` as a free-form name, with the ellipsoid
+    exposed via separate GeoKeys (``GeogSemiMajorAxisGeoKey``,
+    ``GeogSemiMinorAxisGeoKey``, ``GeogInvFlatteningGeoKey``).
+    Rasterio reads such files and reports a GEOGCS CRS; pre-#1930
+    xrspatial dropped the projection silently because nothing stamped
+    a canonical ``attrs['crs_wkt']``. See issue #1930 (Phase 3
+    ``crs_citation_only`` parity gap).
+
+    Returns a WKT string when pyproj is installed and the GeoKeys
+    expose enough of the ellipsoid to feed ``pyproj.CRS.from_dict``,
+    otherwise ``None``. The synthesized CRS is name-stripped (PROJ
+    drops the GEOGCS name on ``to_dict()`` anyway), so callers that
+    need the user-supplied name should still consult
+    ``attrs['geog_citation']`` -- this helper exists to close the
+    canonical-CRS parity gap, not to round-trip the citation field.
+
+    Angular units are not threaded through. PROJ's ``longlat`` always
+    emits degrees, and the corpus has no radian-unit user-defined
+    fixture; if one shows up, this helper needs an angular_units_code
+    parameter and the proj_dict needs a ``to_meter`` / units shim.
+
+    Non-geographic model types fall through to ``None``:
+
+    * ``MODEL_TYPE_PROJECTED`` (1) user-defined CRSes (``ProjectedCSType
+      == 32767``) need the GeogPrime / Projection parameters that this
+      helper does not yet read.
+    * ``MODEL_TYPE_GEOCENTRIC`` (3) and unknown / zero model types are
+      not exercised by the corpus and would need their own per-type
+      proj_dict construction.
+
+    In all three cases the caller falls through to the existing
+    deprecated-attrs path.
+    """
+    if model_type != MODEL_TYPE_GEOGRAPHIC:
+        return None
+    if semi_major is None:
+        return None
+    # No ellipsoid: a degenerate file. Refuse to fabricate one rather
+    # than emit a misleading WGS84 fallback.
+    if semi_minor is None and inv_flattening is None:
+        return None
+    try:
+        from pyproj import CRS
+    except ImportError:
+        return None
+
+    proj_dict: dict[str, object] = {'proj': 'longlat', 'no_defs': True}
+    proj_dict['a'] = float(semi_major)
+    if inv_flattening is not None and float(inv_flattening) != 0.0:
+        proj_dict['rf'] = float(inv_flattening)
+    elif semi_minor is not None:
+        proj_dict['b'] = float(semi_minor)
+    else:
+        # inv_flattening == 0 is the GeoTIFF convention for a sphere
+        # (b == a). Encode it as such for PROJ.
+        proj_dict['b'] = float(semi_major)
+
+    try:
+        crs = CRS.from_dict(proj_dict)
+        return crs.to_wkt()
+    except Exception as exc:
+        import warnings
+        from . import _geotiff_strict_mode, GeoTIFFFallbackWarning
+        if _geotiff_strict_mode():
+            raise
+        warnings.warn(
+            f"_synthesize_user_defined_wkt failed "
+            f"({type(exc).__name__}: {exc}); returning None.",
+            GeoTIFFFallbackWarning,
+            stacklevel=2,
+        )
+        return None
+
+
 def _parse_geokeys(ifd: IFD, data: bytes | memoryview,
                    byte_order: str) -> dict[int, int | float | str]:
     """Parse the GeoKeyDirectory and resolve values from param tags.
@@ -609,6 +696,12 @@ def extract_geo_info(ifd: IFD, data: bytes | memoryview,
     semi_major = geokeys.get(GEOKEY_GEOG_SEMI_MAJOR_AXIS)
     if not isinstance(semi_major, (int, float)):
         semi_major = None
+    # GeogSemiMinorAxisGeoKey. Kept local for WKT synthesis below; not
+    # surfaced as a separate attr because the canonical ellipsoid lives
+    # in attrs['crs_wkt'].
+    semi_minor = geokeys.get(GEOKEY_GEOG_SEMI_MINOR_AXIS)
+    if not isinstance(semi_minor, (int, float)):
+        semi_minor = None
     inv_flat = geokeys.get(GEOKEY_GEOG_INV_FLATTENING)
     if not isinstance(inv_flat, (int, float)):
         inv_flat = None
@@ -728,6 +821,23 @@ def extract_geo_info(ifd: IFD, data: bytes | memoryview,
         # to_geotiff only consults attrs['crs'] / attrs['crs_wkt']. See
         # issue #1632.
         crs_wkt = crs_name
+    else:
+        # Citation-only user-defined geographic CRS: no EPSG, no WKT
+        # in the citation, but enough ellipsoid / units GeoKeys to
+        # synthesize a canonical WKT (issue #1930 Phase 3). Without
+        # this branch ``attrs['crs_wkt']`` stays None and the
+        # golden-corpus oracle reports a CRS parity gap.
+        synth = _synthesize_user_defined_wkt(
+            model_type=(
+                int(model_type)
+                if isinstance(model_type, (int, float)) else 0
+            ),
+            semi_major=semi_major,
+            semi_minor=semi_minor,
+            inv_flattening=inv_flat,
+        )
+        if synth is not None:
+            crs_wkt = synth
 
     return GeoInfo(
         transform=transform,

xrspatial/geotiff/tests/golden_corpus/test_oracle.py

@@ -517,6 +517,64 @@ def test_crs_citation_only_fixture_rejects_unrelated_crs() -> None:
         compare_to_oracle(path, cand)
 
 
+def test_crs_citation_only_open_geotiff_stamps_canonical_wkt() -> None:
+    """``open_geotiff`` on the citation fixture stamps ``attrs['crs_wkt']``.
+
+    The user-defined geographic CRS in ``crs_citation_only.tif`` has no
+    EPSG code and no WKT in the citation; only the ellipsoid radius,
+    inverse flattening, and angular-units GeoKeys are populated. The
+    reader synthesizes a canonical WKT from those parameters via
+    :func:`xrspatial.geotiff._geotags._synthesize_user_defined_wkt`,
+    which closes the Phase 3 ``crs_citation_only`` parity gap.
+
+    Pinned here so any future refactor that drops the synthesis branch
+    re-opens the gap visibly. The companion oracle test
+    (``test_crs_citation_only_xrspatial_round_trips_through_oracle``)
+    drives the same fixture through ``compare_to_oracle`` to confirm the
+    stamp lands somewhere the parity check accepts.
+    """
+    from xrspatial.geotiff import open_geotiff
+
+    path = _CRS_FIXTURE_DIR / 'crs_citation_only.tif'
+    da = open_geotiff(str(path))
+
+    # The fixture has no EPSG, so the canonical EPSG attr stays absent.
+    assert da.attrs.get('crs') is None
+    wkt = da.attrs.get('crs_wkt')
+    assert isinstance(wkt, str) and wkt, (
+        f"open_geotiff must stamp a non-empty crs_wkt on the "
+        f"citation-only fixture; got {wkt!r}"
+    )
+    # PROJ structural sanity: the synthesized WKT must parse back
+    # through rasterio's CRS.from_wkt (which delegates to PROJ).
+    parsed = rasterio.crs.CRS.from_wkt(wkt)
+    assert parsed is not None
+    # And it must share PROJ-dict shape with the rasterio reference.
+    with rasterio.open(path) as src:
+        ref = src.crs
+    assert parsed.to_dict() == ref.to_dict(), (
+        f"synthesized WKT must round-trip to the same PROJ dict as "
+        f"the rasterio reference; got {parsed.to_dict()} vs "
+        f"{ref.to_dict()}"
+    )
+
+
+def test_crs_citation_only_xrspatial_round_trips_through_oracle() -> None:
+    """``compare_to_oracle`` accepts the xrspatial-stamped citation CRS.
+
+    Drives the citation fixture through ``open_geotiff`` (the exact
+    code path the Phase 3 backend parametrizations use) and runs the
+    result through ``compare_to_oracle``. This is the end-to-end
+    parity check that flips the corpus from xfail to pass once
+    ``_synthesize_user_defined_wkt`` is wired in.
+    """
+    from xrspatial.geotiff import open_geotiff
+
+    path = _CRS_FIXTURE_DIR / 'crs_citation_only.tif'
+    da = open_geotiff(str(path))
+    compare_to_oracle(path, da)
+
+
 def test_crs_wkt_utm10n_fixture_accepts_wkt_attr() -> None:
     """``crs_wkt_utm10n`` also accepts a candidate that carries crs_wkt.
 

xrspatial/geotiff/tests/test_golden_corpus_dask_gpu_1930.py

@@ -13,6 +13,14 @@
 the eager / dask / GPU modules. ``_DASK_GPU_SKIPS`` is reserved for
 gaps that surface only when chunked GPU reads stitch through nvCOMP
 plus dask; it starts empty.
+
+Resolved gaps (no longer xfail):
+
+* Integer nodata masking -- closed by the oracle's
+  ``_normalise_for_masked_nodata`` helper.
+* ``crs_citation_only`` -- closed by ``_synthesize_user_defined_wkt``
+  in ``_geotags.py``, which stamps a canonical
+  ``attrs['crs_wkt']`` on user-defined geographic CRSes.
 """
 from __future__ import annotations
 
@@ -53,15 +61,13 @@
 # Integer-nodata masking used to live here too; the oracle's
 # _normalise_for_masked_nodata helper (#2046) closes that gap so it is
 # no longer xfailed on any backend. The multi-band axis-order gap for
-# the JPEG-YCbCr fixture is also closed (see ``_normalise_axis_order``
-# in ``_oracle.py``).
-_PARITY_GAPS: dict[str, str] = {
-    "crs_citation_only": (
-        "citation-only CRS: xrspatial decodes the citation into deprecated "
-        "attrs['geog_citation'] but does not emit a canonical attrs['crs'] "
-        "or attrs['crs_wkt']. Real parity gap; needs a fix in _crs.py."
-    ),
-}
+# the JPEG-YCbCr fixture is also closed on the eager / dask paths
+# (see ``_normalise_axis_order`` in ``_oracle.py``); on the dask+GPU
+# pipeline the decode error wins first so the JPEG fixture lives in
+# ``_INTENTIONAL_SKIPS`` below. crs_citation_only moved off this list
+# once ``_synthesize_user_defined_wkt`` started stamping a canonical
+# WKT for user-defined geographic CRSes.
+_PARITY_GAPS: dict[str, str] = {}
 
 # Empty: failures unique to the dask+GPU combo (eager / dask / pure-
 # GPU all pass) would land here. Kept around as the documented home

xrspatial/geotiff/tests/test_golden_corpus_dask_numpy_1930.py

@@ -10,12 +10,20 @@
 shows up here.
 
 The skip / xfail taxonomy is intentionally identical to the eager
-module's. The parity gaps that motivated them (integer-nodata masking,
-citation CRS encoding, RGB band axis order, MinIsWhite inversion) live
-in the codec / attrs layer, which both backends share. If a future
-parity gap turns out to be dask-specific, add it to ``_DASK_SKIPS``
-rather than ``_PARITY_GAPS`` so the difference between "shared
-codec/attrs gap" and "dask plumbing gap" stays legible.
+module's. The parity gaps that motivated them (RGB band axis order,
+MinIsWhite inversion) live in the codec / attrs layer, which both
+backends share. If a future parity gap turns out to be dask-specific,
+add it to ``_DASK_SKIPS`` rather than ``_PARITY_GAPS`` so the
+difference between "shared codec/attrs gap" and "dask plumbing gap"
+stays legible.
+
+Resolved gaps (no longer xfail):
+
+* Integer nodata masking -- closed by the oracle's
+  ``_normalise_for_masked_nodata`` helper.
+* ``crs_citation_only`` -- closed by ``_synthesize_user_defined_wkt``
+  in ``_geotags.py``, which stamps a canonical
+  ``attrs['crs_wkt']`` on user-defined geographic CRSes.
 
 The chunk size is fixed at 32. Most corpus fixtures are 64x64 or
 smaller, so 32 produces either a 2x2 chunk grid or a single chunk;
@@ -52,14 +60,11 @@
 # oracle's _normalise_for_masked_nodata helper closes that gap so it
 # is no longer xfailed on any backend. The multi-band axis-order gap
 # for the JPEG-YCbCr fixture is also closed (see
-# ``_normalise_axis_order`` in ``_oracle.py``).
-_PARITY_GAPS: dict[str, str] = {
-    "crs_citation_only": (
-        "citation-only CRS: xrspatial decodes the citation into deprecated "
-        "attrs['geog_citation'] but does not emit a canonical attrs['crs'] "
-        "or attrs['crs_wkt']. Real parity gap; needs a fix in _crs.py."
-    ),
-}
+# ``_normalise_axis_order`` in ``_oracle.py``). The citation-only CRS
+# gap is closed by ``_synthesize_user_defined_wkt`` in ``_geotags.py``,
+# which stamps a canonical WKT for user-defined geographic CRSes onto
+# ``attrs['crs_wkt']``.
+_PARITY_GAPS: dict[str, str] = {}
 
 # Dask-only gaps go here. Empty in the first pass; add an entry only when
 # a fixture is dask-specific (i.e. eager passes, dask does not).

xrspatial/geotiff/tests/test_golden_corpus_eager_numpy_1930.py

@@ -24,10 +24,9 @@
 
 Real parity gaps (``xfail``):
 
-* ``crs_citation_only`` -- xrspatial decodes the citation into the
-  deprecated ``attrs['geog_citation']`` but does not emit a canonical
-  ``attrs['crs']`` or ``attrs['crs_wkt']``. Real parity gap; needs a fix
-  in ``_crs.py`` to round-trip citation WKT.
+None today. Every parity gap the corpus has surfaced so far has either
+been fixed (see Resolved gaps below) or moved to a per-backend skip
+table because it is not a shared decode/attrs gap.
 
 Resolved gaps (no longer xfail):
 
@@ -43,6 +42,12 @@
   The oracle's ``_normalise_axis_order`` helper now transposes the
   trailing-band candidate to leading-band before the shape and pixel
   checks run, so this fixture passes directly.
+* ``crs_citation_only`` -- the GeoTIFF reader now synthesizes a
+  canonical WKT for user-defined geographic CRSes from the ellipsoid
+  / units GeoKeys and stamps it on ``attrs['crs_wkt']`` (see
+  ``_synthesize_user_defined_wkt`` in ``_geotags.py``). The oracle's
+  ``_candidate_crs`` picks it up via the WKT branch, and ``_crs_equal``
+  compares the PROJ dicts because neither side has an EPSG code.
 
 Intentional skip (``skip``):
 
@@ -78,13 +83,7 @@
 )
 
 
-_PARITY_GAPS: dict[str, str] = {
-    "crs_citation_only": (
-        "citation-only CRS: xrspatial decodes the citation into deprecated "
-        "attrs['geog_citation'] but does not emit a canonical attrs['crs'] "
-        "or attrs['crs_wkt']. Real parity gap; needs a fix in _crs.py."
-    ),
-}
+_PARITY_GAPS: dict[str, str] = {}
 
 _INTENTIONAL_SKIPS: dict[str, str] = {
     "nodata_miniswhite_uint8": (

xrspatial/geotiff/tests/test_golden_corpus_gpu_1930.py

@@ -11,13 +11,22 @@
 this module; runs with a GPU exercise every fixture the eager and dask
 backends already do.
 
-``_PARITY_GAPS`` carries over the codec / attrs gaps that all three
-backends share (citation CRS, integer nodata masking). ``_GPU_SKIPS``
-holds GPU-only failures, currently the JPEG-YCbCr fixture: the GPU
-decoder does not handle it and ``on_gpu_failure='strict'`` raises
-rather than falling back, so the read fails before the oracle can
-compare. On eager / dask the same fixture exposes the RGB axis-order
-divergence; on GPU strict mode it never gets that far.
+``_PARITY_GAPS`` is empty for the GPU backend now that integer nodata
+masking and the citation-only CRS gap are closed at the codec / attrs
+layer. ``_GPU_SKIPS`` holds GPU-only failures, currently the
+JPEG-YCbCr fixture: the GPU decoder does not handle it and
+``on_gpu_failure='strict'`` raises rather than falling back, so the
+read fails before the oracle can compare. On eager / dask the same
+fixture exposes the RGB axis-order divergence; on GPU strict mode it
+never gets that far.
+
+Resolved gaps (no longer xfail):
+
+* Integer nodata masking -- closed by the oracle's
+  ``_normalise_for_masked_nodata`` helper.
+* ``crs_citation_only`` -- closed by ``_synthesize_user_defined_wkt``
+  in ``_geotags.py``, which stamps a canonical
+  ``attrs['crs_wkt']`` on user-defined geographic CRSes.
 
 The GPU read is configured with ``on_gpu_failure='strict'`` so a codec
 that would silently CPU-fall-back instead surfaces as an xfail / fail
@@ -56,16 +65,10 @@
 )
 
 
-# Integer-nodata masking used to live here too; the oracle's
-# _normalise_for_masked_nodata helper now closes that gap so it is no
-# longer xfailed on any backend.
-_PARITY_GAPS: dict[str, str] = {
-    "crs_citation_only": (
-        "citation-only CRS: xrspatial decodes the citation into deprecated "
-        "attrs['geog_citation'] but does not emit a canonical attrs['crs'] "
-        "or attrs['crs_wkt']. Real parity gap; needs a fix in _crs.py."
-    ),
-}
+# Integer-nodata masking used to live here (closed by the oracle's
+# _normalise_for_masked_nodata helper), and crs_citation_only used to
+# live here too (closed by _synthesize_user_defined_wkt in _geotags.py).
+_PARITY_GAPS: dict[str, str] = {}
 
 # GPU-only gaps. Empty in the current corpus: failures here would be
 # fixtures the eager and dask backends decode cleanly but the GPU