ENH: Add CRS.to_2d() (#1267)

docs/advanced_examples.rst

@@ -155,6 +155,28 @@ In PROJ 6+ you need to explicitly change your CRS to 3D if you have
     >>> transformer_3d.transform(8.37909, 47.01987, 1000)
     (2671499.8913080636, 1208075.1135782297, 951.4265527743846)
 
+
+Demote CRS to 2D
+----------------
+
+.. versionadded:: 3.6
+
+
+With the need for explicit 3D CRS since PROJ 6+, one might need to retrieve their 2D version,
+for example to create another 3D CRS compound between a 2D CRS and a vertical CRS.
+
+.. code-block:: python
+
+    >>> from pyproj import CRS, Transformer
+    >>> from pyproj.crs import CompoundCRS
+    >>> src_crs = CRS("EPSG:4979") # Any 3D CRS, here the 3D WGS 84
+    >>> vert_crs = CRS("EPSG:5773") # Any vertical CRS, here the EGM96 geoid
+    >>> dst_crs = CompoundCRS(src_crs.name + vert_crs.name, components=[src_crs.to_2d(), vert_crs])
+    >>> transformer_3d = Transformer.from_crs(src_crs, dst_crs, always_xy=True)
+    >>> transformer_3d.transform(8.37909, 47.01987, 1000)
+    (8.37909, 47.01987, 951.7851086745321)
+
+
 Projected CRS Bounds
 ----------------------
 

docs/history.rst

@@ -4,6 +4,8 @@ Change Log
 Latest
 ------
 
+- ENH: Added :meth:`CRS.to_2d` to demote 3D CRS to 2D (issue #1266)
+
 3.5.0
 ------
 - DEP: Minimum PROJ version 9.0 (issue #1223)

pyproj/_crs.pxd

@@ -9,6 +9,9 @@ cdef extern from "proj_experimental.h":
                                const char* crs_3D_name,
                                const PJ* crs_2D)
 
+    PJ *proj_crs_demote_to_2D(PJ_CONTEXT *ctx,
+                              const char *crs_2D_name,
+                              const PJ *crs_3D)
 
 cdef tuple _get_concatenated_operations(PJ_CONTEXT*context, PJ*concatenated_operation)
 cdef _to_proj4(

pyproj/_crs.pyi

@@ -227,6 +227,7 @@ class _CRS(Base):
         self, auth_name: Optional[str] = None, min_confidence: int = 70
     ) -> List[AuthorityMatchInfo]: ...
     def to_3d(self, name: Optional[str] = None) -> "_CRS": ...
+    def to_2d(self, name: Optional[str] = None) -> "_CRS": ...
     @property
     def is_geographic(self) -> bool: ...
     @property

pyproj/_crs.pyx

@@ -2965,6 +2965,44 @@ cdef class _CRS(Base):
             proj_destroy(projobj)
         return crs_3d
 
+    def to_2d(self, str name=None):
+        """
+        .. versionadded:: 3.6.0
+
+        Convert the current CRS to the 2D version if it makes sense.
+
+        Parameters
+        ----------
+        name: str, optional
+            CRS name. If None, it will use the name of the original CRS.
+
+        Returns
+        -------
+        _CRS
+        """
+        cdef char* c_name = NULL
+        cdef bytes b_name
+        if name is not None:
+            b_name = cstrencode(name)
+            c_name = b_name
+
+        cdef PJ * projobj = proj_crs_demote_to_2D(
+            self.context, c_name, self.projobj
+        )
+        _clear_proj_error()
+        if projobj == NULL:
+            return self
+        try:
+            crs_2d = _CRS(_to_wkt(
+                self.context,
+                projobj,
+                version=WktVersion.WKT2_2019,
+                pretty=False,
+            ))
+        finally:
+            proj_destroy(projobj)
+        return crs_2d
+
     def _is_crs_property(
         self, str property_name, tuple property_types, int sub_crs_index=0
     ):

pyproj/crs/crs.py

@@ -1433,6 +1433,23 @@ def to_3d(self, name: Optional[str] = None) -> "CRS":
         """
         return self.__class__(self._crs.to_3d(name=name))
 
+    def to_2d(self, name: Optional[str] = None) -> "CRS":
+        """
+        .. versionadded:: 3.6.0
+
+        Convert the current CRS to the 2D version if it makes sense.
+
+        Parameters
+        ----------
+        name: str, optional
+            CRS name. Defaults to use the name of the original CRS.
+
+        Returns
+        -------
+        CRS
+        """
+        return self.__class__(self._crs.to_2d(name=name))
+
     @property
     def is_geographic(self) -> bool:
         """

test/crs/test_crs.py

@@ -1608,6 +1608,36 @@ def test_to_3d__name():
     assert crs_3d.name == "TEST"
 
 
+@pytest.mark.parametrize(
+    "crs_input",
+    [
+        CRS("EPSG:4979"),  # native 3D
+        CRS("EPSG:2056").to_3d(),  # a 2D CRS converted to 3D
+        CRS("EPSG:4326+5773"),  # a 3D CRS based on a compound
+    ],
+)
+def test_to_2d(crs_input):
+    assert len(crs_input.axis_info) == 3
+    horizon_axis_crs_3d = crs_input.axis_info[:-1]
+    crs_2d = crs_input.to_2d()
+    horizon_axis_crs_2d = crs_input.axis_info
+    assert len(crs_2d.axis_info) == 2
+    assert horizon_axis_crs_2d[0] == horizon_axis_crs_3d[0]
+    assert horizon_axis_crs_2d[1] == horizon_axis_crs_3d[1]
+    assert crs_2d.to_2d() == crs_2d
+    # For CompoundCRS, the 3D name is initialized different from 2D
+    if crs_input.name == "WGS 84 + EGM96 height":
+        assert crs_2d.name == "WGS 84"
+    # Otherwise, no change
+    else:
+        assert crs_2d.name == crs_input.name
+
+
+def test_to_2d__name():
+    crs_2d = CRS("EPSG:2056").to_3d().to_2d(name="TEST")
+    assert crs_2d.name == "TEST"
+
+
 def test_crs__pickle(tmp_path):
     assert_can_pickle(CRS("epsg:4326"), tmp_path)