Reject compound EPSG codes on the stable-EPSG writer path (#2418)

docs/source/user_guide/geotiff_safe_io.rst

@@ -235,6 +235,13 @@ the ambiguous-metadata family at once.
      - The writer cannot classify an EPSG code as geographic or
        projected.
      - Pass a code pyproj can resolve, or install pyproj.
+   * - :class:`~xrspatial.geotiff.NonRepresentableEPSGCRSError`
+     - The integer EPSG code resolves to a compound (horizontal +
+       vertical) CRS, which the writer cannot represent in a single
+       ``GeographicTypeGeoKey`` or ``ProjectedCSTypeGeoKey`` slot.
+     - Pass the full compound CRS as WKT to take the user-defined CRS
+       fallback path, or pass the horizontal sub-CRS EPSG directly if
+       the vertical component is not needed.
    * - :class:`~xrspatial.geotiff.UnsupportedGeoTIFFFeatureError`
      - The input declares a feature the GeoTIFF module does not
        implement (warped / reprojection VRTs, pansharpened or derived

xrspatial/geotiff/__init__.py

@@ -72,9 +72,9 @@
     transform_tuple_from_pixel_geometry as _transform_tuple_from_pixel_geometry  # noqa: F401
 from ._crs import _resolve_crs_to_wkt, _wkt_to_epsg  # noqa: F401
 from ._errors import (ConflictingCRSError, ConflictingNodataError, GeoTIFFAmbiguousMetadataError,
-                      InvalidCRSCodeError, MixedBandMetadataError, NonUniformCoordsError,
-                      RotatedTransformError, UnknownCRSModelTypeError, UnparseableCRSError,
-                      UnsupportedGeoTIFFFeatureError)
+                      InvalidCRSCodeError, MixedBandMetadataError, NonRepresentableEPSGCRSError,
+                      NonUniformCoordsError, RotatedTransformError, UnknownCRSModelTypeError,
+                      UnparseableCRSError, UnsupportedGeoTIFFFeatureError)
 from ._geotags import RASTER_PIXEL_IS_AREA, RASTER_PIXEL_IS_POINT, GeoTransform  # noqa: F401
 from ._reader import _MAX_CLOUD_BYTES_SENTINEL, UnsafeURLError
 from ._reader import read_to_array as _read_to_array
@@ -110,6 +110,7 @@
     'GEOREF_STATUS_VALUES',
     'InvalidCRSCodeError',
     'MixedBandMetadataError',
+    'NonRepresentableEPSGCRSError',
     'NonUniformCoordsError',
     'RotatedTransformError',
     'SUPPORTED_FEATURES',

xrspatial/geotiff/_crs.py

@@ -15,7 +15,7 @@
 import numbers
 import warnings
 
-from ._errors import UnparseableCRSError
+from ._errors import NonRepresentableEPSGCRSError, UnparseableCRSError
 from ._runtime import GeoTIFFFallbackWarning, _geotiff_strict_mode
 
 #: WKT root keywords. A string that starts with one of these (after
@@ -45,6 +45,49 @@ def _looks_like_wkt(s: str) -> bool:
     return s.lstrip().upper().startswith(_WKT_ROOT_KEYWORDS)
 
 
+def _reject_non_representable_epsg(crs_int: int, crs_obj) -> None:
+    """Refuse compound EPSG codes that the stable EPSG writer cannot
+    round-trip through rasterio / GDAL (#2418).
+
+    The GeoTIFF writer's stable-path emits only two EPSG-carrying
+    GeoKeys: ``GeographicTypeGeoKey`` (2048) and
+    ``ProjectedCSTypeGeoKey`` (3072). Both are spec'd to hold a 2D
+    horizontal CRS code. When the caller passes a compound EPSG code
+    (horizontal + vertical, e.g. EPSG:6349 = "NAD83(2011) + NAVD88
+    height"), pyproj reports the compound CRS as ``is_geographic`` or
+    ``is_projected`` based on its horizontal sub-CRS, and the writer
+    stores the compound code in the horizontal-CRS slot. The resulting
+    GeoTIFF either reads back through GDAL as the wrong CRS (the
+    horizontal sub-CRS, with the vertical component silently dropped)
+    or as no CRS at all, because the GeoKey value claims to be a 2D
+    horizontal CRS but the EPSG database registers it as compound.
+
+    Reject the code at validation time so the failure surfaces at the
+    write call instead of after the file is on disk. Callers who need
+    to preserve the compound semantics can pass the full WKT (which
+    takes the user-defined / citation path -- see issue #1768) or
+    pass the horizontal sub-CRS code directly.
+
+    pyproj exposes ``is_compound`` on every CRS object since 2.x;
+    callers without pyproj never reach this function (the caller
+    site short-circuits on ``ImportError``).
+    """
+    if getattr(crs_obj, "is_compound", False):
+        raise NonRepresentableEPSGCRSError(
+            f"EPSG:{crs_int} is a compound CRS ({crs_obj.name!r}) and "
+            "cannot be written through the integer-EPSG path: the "
+            "GeoTIFF writer only emits GeographicTypeGeoKey / "
+            "ProjectedCSTypeGeoKey, which are specified to carry a 2D "
+            "horizontal CRS code. Storing a compound EPSG there "
+            "produces a file that rasterio / GDAL reads back as the "
+            "horizontal sub-CRS (vertical component lost) or as no "
+            "CRS at all. Pass the full compound WKT as crs= to take "
+            "the user-defined CRS path (see issue #1768), or pass the "
+            "horizontal sub-CRS EPSG directly if the vertical "
+            "component is not needed."
+        )
+
+
 def _validate_crs_arg(crs) -> None:
     """Reject malformed ``crs=`` arguments before they reach the writer.
 
@@ -90,7 +133,7 @@ def _validate_crs_arg(crs) -> None:
         except ImportError:
             return
         try:
-            CRS.from_epsg(crs_int)
+            crs_obj = CRS.from_epsg(crs_int)
         except Exception as e:
             if _geotiff_strict_mode():
                 raise
@@ -99,6 +142,7 @@ def _validate_crs_arg(crs) -> None:
                 f"(pyproj: {type(e).__name__}: {e}). Pass a valid "
                 f"EPSG integer, a WKT string, or None."
             ) from e
+        _reject_non_representable_epsg(crs_int, crs_obj)
         return
     if isinstance(crs, str):
         return
@@ -122,6 +166,15 @@ def _wkt_to_epsg(wkt_or_proj: str) -> int | None:
         from pyproj import CRS
         crs = CRS.from_user_input(wkt_or_proj)
         epsg = crs.to_epsg()
+        # Issue #2418: ``epsg`` is set when pyproj recognises the WKT,
+        # but if the resolved CRS is compound (horizontal + vertical),
+        # the integer-EPSG writer path cannot represent it -- the GeoKey
+        # slots only carry 2D horizontal codes. Return None so the
+        # caller routes through the WKT / citation fallback instead of
+        # writing a compound EPSG into a horizontal slot. The original
+        # WKT is preserved on the fallback path.
+        if epsg is not None and getattr(crs, "is_compound", False):
+            return None
         return epsg
     except Exception as e:
         if _geotiff_strict_mode():

xrspatial/geotiff/_errors.py

@@ -136,6 +136,28 @@ class UnknownCRSModelTypeError(GeoTIFFAmbiguousMetadataError):
     """
 
 
+class NonRepresentableEPSGCRSError(GeoTIFFAmbiguousMetadataError):
+    """EPSG code resolves, but the GeoTIFF writer cannot represent it (#2418).
+
+    Raised by the GeoTIFF writer when the caller supplies an integer
+    EPSG code that pyproj resolves to a compound CRS (horizontal +
+    vertical). The writer can only emit ``GeographicTypeGeoKey`` (2048)
+    or ``ProjectedCSTypeGeoKey`` (3072); storing a compound EPSG in
+    either slot produces a GeoTIFF that rasterio / GDAL reads back as
+    a different CRS (or no CRS at all), because the GeoKey value
+    declares "this is a 2D horizontal CRS" but the EPSG database
+    registers it as compound.
+
+    Callers who need to preserve a compound CRS should either pass the
+    horizontal sub-CRS as the integer EPSG (the vertical component is
+    not representable in a single GeoTIFF without additional GeoKeys
+    that this release does not emit), or pass the full compound WKT as
+    ``crs=`` so the writer takes the ``ProjectedCSType=32767`` /
+    ``GTCitationGeoKey`` fallback path. The WKT fallback path round-trips
+    through libgeotiff and GDAL but not all readers; see issue #1768.
+    """
+
+
 class UnsupportedGeoTIFFFeatureError(ValueError):
     """Caller asked for a feature this release does not implement (#2349).
 
@@ -163,5 +185,6 @@ class UnsupportedGeoTIFFFeatureError(ValueError):
     "ConflictingNodataError",
     "VRTUnsupportedError",
     "UnknownCRSModelTypeError",
+    "NonRepresentableEPSGCRSError",
     "UnsupportedGeoTIFFFeatureError",
 ]

xrspatial/geotiff/_geotags.py

@@ -4,7 +4,7 @@
 from dataclasses import dataclass, field
 
 from ._dtypes import resolve_bits_per_sample
-from ._errors import RotatedTransformError, UnknownCRSModelTypeError
+from ._errors import NonRepresentableEPSGCRSError, RotatedTransformError, UnknownCRSModelTypeError
 from ._header import (IFD, TAG_BITS_PER_SAMPLE, TAG_COMPRESSION, TAG_EXTRA_SAMPLES,
                       TAG_GDAL_METADATA, TAG_GDAL_NODATA, TAG_GEO_ASCII_PARAMS,
                       TAG_GEO_DOUBLE_PARAMS, TAG_GEO_KEY_DIRECTORY, TAG_IMAGE_LENGTH,
@@ -1452,6 +1452,24 @@ def _model_type_from_epsg(crs_epsg: int) -> int:
                 "known EPSG or installing a current pyproj database will "
                 "resolve this."
             ) from e
+        # Defense in depth (#2418): the stable-EPSG path can only emit
+        # GeographicTypeGeoKey or ProjectedCSTypeGeoKey, both of which
+        # are specified to hold a 2D horizontal CRS. A compound CRS
+        # (horizontal + vertical) reports ``is_geographic`` or
+        # ``is_projected`` for its horizontal sub-CRS but is not
+        # representable in either slot -- rasterio / GDAL reads the
+        # output back as the wrong CRS. ``_validate_crs_arg`` should
+        # have rejected this already; trip the same error here so a
+        # direct call into ``build_geo_tags`` cannot bypass the check.
+        if getattr(crs, "is_compound", False):
+            raise NonRepresentableEPSGCRSError(
+                f"EPSG:{crs_epsg} is a compound CRS ({crs.name!r}); "
+                "the GeoTIFF integer-EPSG writer path only emits 2D "
+                "horizontal GeoKeys. Pass the compound CRS as WKT "
+                "(``crs=<wkt>``) for the user-defined CRS fallback, "
+                "or pass the horizontal sub-CRS EPSG directly. "
+                "See issue #2418."
+            )
         if crs.is_geographic:
             return MODEL_TYPE_GEOGRAPHIC
         return MODEL_TYPE_PROJECTED