From c36fbb130c878ffc0be444e739e6451a7df96fc0 Mon Sep 17 00:00:00 2001 From: snowman2 Date: Sat, 27 Mar 2021 20:05:59 -0500 Subject: [PATCH] ENH: Added pyproj.crs.CRS.to_3d --- docs/advanced_examples.rst | 25 +++++++++++++++++++++++++ docs/history.rst | 2 +- pyproj/_crs.pxd | 5 +++++ pyproj/_crs.pyi | 1 + pyproj/_crs.pyx | 37 +++++++++++++++++++++++++++++++++++++ pyproj/crs/crs.py | 22 ++++++++++++++++++++++ test/crs/test_crs.py | 25 +++++++++++++++++++++++++ 7 files changed, 116 insertions(+), 1 deletion(-) diff --git a/docs/advanced_examples.rst b/docs/advanced_examples.rst index 4a6a1c10a..b4f3912ab 100644 --- a/docs/advanced_examples.rst +++ b/docs/advanced_examples.rst @@ -137,6 +137,31 @@ which transformation operation is selected in the transformation. - bounds: (-136.46, 49.0, -60.72, 83.17) +Promote CRS to 3D +------------------- + +.. versionadded:: 3.1 + + +In PROJ 6+ you neeed to explictly change your CRS to 3D if you have +2D CRS and you want the ellipsoidal height taken into account. + + +.. code-block:: python + + >>> from pyproj import CRS, Transformer + >>> transformer = Transformer.from_crs("EPSG:4326", "EPSG:2056", always_xy=True) + >>> transformer.transform(8.37909, 47.01987, 1000) + (2671499.8913080636, 1208075.1135782297, 1000.0) + >>> transformer_3d = Transformer.from_crs( + ... CRS("EPSG:4326").to_3d(), + ... CRS("EPSG:2056").to_3d(), + ... always_xy=True, + ...) + >>> transformer_3d.transform(8.37909, 47.01987, 1000) + (2671499.8913080636, 1208075.1135782297, 951.4265527743846) + + Multithreading -------------- diff --git a/docs/history.rst b/docs/history.rst index b1de6e6a0..3152d2609 100644 --- a/docs/history.rst +++ b/docs/history.rst @@ -10,7 +10,7 @@ Change Log * ENH: Added :meth:`pyproj.transformer.Transformer.to_proj4` (pull #798) * ENH: Added authority, accuracy, and allow_ballpark kwargs to :meth:`pyproj.transformer.Transformer.from_crs` (issue #754) * ENH: Added support for "AUTH:CODE" input to :meth:`pyproj.transformer.Transformer.from_pipeline` (issue #755) - +* ENH: Added :meth:`pyproj.crs.CRS.to_3d` (pull #808) 3.0.1 ----- diff --git a/pyproj/_crs.pxd b/pyproj/_crs.pxd index 9cedcb1ad..4cc78d6bc 100644 --- a/pyproj/_crs.pxd +++ b/pyproj/_crs.pxd @@ -1,5 +1,10 @@ include "proj.pxi" +cdef extern from "proj_experimental.h": + PJ *proj_crs_promote_to_3D(PJ_CONTEXT *ctx, + const char* crs_3D_name, + const PJ* crs_2D) + cdef _get_concatenated_operations(PJ_CONTEXT*context, PJ*concatenated_operation) cdef _to_proj4( diff --git a/pyproj/_crs.pyi b/pyproj/_crs.pyi index 1d2825f08..5a476eb9e 100644 --- a/pyproj/_crs.pyi +++ b/pyproj/_crs.pyi @@ -218,6 +218,7 @@ class _CRS(Base): def to_authority( self, auth_name: Optional[str] = None, min_confidence: int = 70 ): ... + def to_3d(self, name: Optional[str] = None) -> "_CRS": ... @property def is_geographic(self) -> bool: ... @property diff --git a/pyproj/_crs.pyx b/pyproj/_crs.pyx index 6b3300bbc..12b8c6414 100644 --- a/pyproj/_crs.pyx +++ b/pyproj/_crs.pyx @@ -2782,6 +2782,43 @@ cdef class _CRS(Base): return None + def to_3d(self, name=None): + """ + .. versionadded:: 3.1 + + Convert the current CRS to the 3D version if it makes sense. + + New vertical axis attributes: + - ellipsoidal height + - oriented upwards + - metre units + + Parameters + ---------- + name: str, optional + CRS name. Defaults to use the name of the original CRS. + + Returns + ------- + _CRS + """ + cdef char* name_3D = NULL + if name is not None: + b_name = cstrencode(name) + name_3D = b_name + + cdef PJ * projobj = proj_crs_promote_to_3D( + self.context, name_3D, self.projobj + ) + CRSError.clear() + if projobj == NULL: + return self + try: + crs_3d = _CRS(_to_wkt(self.context, projobj)) + finally: + proj_destroy(projobj) + return crs_3d + def _is_crs_property(self, property_name, property_types, sub_crs_index=0): """ .. versionadded:: 2.2.0 diff --git a/pyproj/crs/crs.py b/pyproj/crs/crs.py index edcc1112f..846aaa7cb 100644 --- a/pyproj/crs/crs.py +++ b/pyproj/crs/crs.py @@ -1294,6 +1294,28 @@ def to_authority(self, auth_name: Optional[str] = None, min_confidence: int = 70 auth_name=auth_name, min_confidence=min_confidence ) + def to_3d(self, name: Optional[str] = None) -> "CRS": + """ + .. versionadded:: 3.1 + + Convert the current CRS to the 3D version if it makes sense. + + New vertical axis attributes: + - ellipsoidal height + - oriented upwards + - metre units + + Parameters + ---------- + name: str, optional + CRS name. Defaults to use the name of the original CRS. + + Returns + ------- + CRS + """ + return CRS(self._crs.to_3d(name=name)) + @property def is_geographic(self) -> bool: """ diff --git a/test/crs/test_crs.py b/test/crs/test_crs.py index aed4e4697..79c27035d 100644 --- a/test/crs/test_crs.py +++ b/test/crs/test_crs.py @@ -1445,3 +1445,28 @@ def test_coordinate_operation__to_proj4__pretty(): proj_string = operation.to_proj4(pretty=True) assert "+proj=pipeline" in proj_string assert "\n" in proj_string + + +@pytest.mark.parametrize( + "crs_input", + [ + "EPSG:4326", + "EPSG:2056", + ], +) +def test_to_3d(crs_input): + crs = CRS(crs_input) + assert len(crs.axis_info) == 2 + crs_3d = crs.to_3d() + assert len(crs_3d.axis_info) == 3 + vert_axis = crs_3d.axis_info[-1] + assert vert_axis.name == "Ellipsoidal height" + assert vert_axis.unit_name == "metre" + assert vert_axis.direction == "up" + assert crs_3d.to_3d() == crs_3d + assert crs_3d.name == crs.name + + +def test_to_3d__name(): + crs_3d = CRS("EPSG:2056").to_3d(name="TEST") + assert crs_3d.name == "TEST"