Merge pull request #1255 from snowman2/to_cf_bug

BUG: Add horizontal_datum_name for geographic CRS in CRS.to_cf
pyproj4 · Mar 24, 2023 · 04ce29a · 04ce29a
2 parents 79a8602 + 9405668
commit 04ce29a
Show file tree

Hide file tree

Showing 4 changed files with 38 additions and 28 deletions.
diff --git a/docs/history.rst b/docs/history.rst
@@ -14,6 +14,8 @@ Latest
 - REF: Raise error when :meth:`.CRS.to_wkt`, :meth:`.CRS.to_json`, or :meth:`.CRS.to_proj4` returns None (issue #1036)
 - CLN: Remove `AzumuthalEquidistantConversion` & :class:`LambertAzumuthalEqualAreaConversion`. :class:`AzimuthalEquidistantConversion` & :class:`LambertAzimuthalEqualAreaConversion` should be used instead (pull #1219)
 - BUG: Fix Derived Projected CRS support (issue #1182)
+- BUG: Add horizontal_datum_name for geographic CRS in :meth:`.CRS.to_cf` (issue #1251)
+- BUG: Add datum ensemble support to :class:`.GeographicCRS` (pull #1255)
 
 3.4.1
 -----

diff --git a/pyproj/crs/_cf1x8.py b/pyproj/crs/_cf1x8.py
@@ -34,6 +34,12 @@
 
 
 def _horizontal_datum_from_params(cf_params):
+    datum_name = cf_params.get("horizontal_datum_name")
+    if datum_name and datum_name not in ("undefined", "unknown"):
+        try:
+            return Datum.from_name(datum_name)
+        except CRSError:
+            pass
     # step 1: build ellipsoid
     ellipsoid = None
     ellipsoid_name = cf_params.get("reference_ellipsoid_name")
@@ -46,7 +52,7 @@ def _horizontal_datum_from_params(cf_params):
             radius=cf_params.get("earth_radius"),
         )
     except CRSError:
-        if ellipsoid_name:
+        if ellipsoid_name and ellipsoid_name not in ("undefined", "unknown"):
             ellipsoid = Ellipsoid.from_name(ellipsoid_name)
 
     # step 2: build prime meridian
@@ -58,19 +64,16 @@ def _horizontal_datum_from_params(cf_params):
             longitude=cf_params["longitude_of_prime_meridian"],
         )
     except KeyError:
-        if prime_meridian_name:
+        if prime_meridian_name and prime_meridian_name not in ("undefined", "unknown"):
             prime_meridian = PrimeMeridian.from_name(prime_meridian_name)
 
     # step 3: build datum
-    datum_name = cf_params.get("horizontal_datum_name")
     if ellipsoid or prime_meridian:
         return CustomDatum(
             name=datum_name or "undefined",
             ellipsoid=ellipsoid or "WGS 84",
             prime_meridian=prime_meridian or "Greenwich",
         )
-    if datum_name:
-        return Datum.from_name(datum_name)
     return None
 
 

diff --git a/pyproj/crs/crs.py b/pyproj/crs/crs.py
@@ -693,6 +693,8 @@ def to_cf(
         # handle geographic CRS
         if self.geodetic_crs:
             cf_dict["geographic_crs_name"] = self.geodetic_crs.name
+            if self.geodetic_crs.datum:
+                cf_dict["horizontal_datum_name"] = self.geodetic_crs.datum.name
 
         if self.is_geographic:
             if self.coordinate_operation:
@@ -711,15 +713,11 @@ def to_cf(
                         self.coordinate_operation.method_name.lower()
                     ](self.coordinate_operation)
                 )
-                if self.datum:
-                    cf_dict["horizontal_datum_name"] = self.datum.name
             else:
                 cf_dict["grid_mapping_name"] = "latitude_longitude"
             return cf_dict
 
         # handle projected CRS
-        if self.is_projected and self.datum:
-            cf_dict["horizontal_datum_name"] = self.datum.name
         coordinate_operation = None
         if not self.is_bound and self.is_projected:
             coordinate_operation = self.coordinate_operation
@@ -1758,31 +1756,35 @@ class GeographicCRS(CustomConstructorCRS):
     def __init__(
         self,
         name: str = "undefined",
-        datum: Any = "urn:ogc:def:datum:EPSG::6326",
+        datum: Any = "urn:ogc:def:ensemble:EPSG::6326",
         ellipsoidal_cs: Optional[Any] = None,
     ) -> None:
         """
         Parameters
         ----------
         name: str, default="undefined"
             Name of the CRS.
-        datum: Any, default="urn:ogc:def:datum:EPSG::6326"
+        datum: Any, default="urn:ogc:def:ensemble:EPSG::6326"
             Anything accepted by :meth:`pyproj.crs.Datum.from_user_input` or
             a :class:`pyproj.crs.datum.CustomDatum`.
         ellipsoidal_cs: Any, optional
             Input to create an Ellipsoidal Coordinate System.
             Anything accepted by :meth:`pyproj.crs.CoordinateSystem.from_user_input`
             or an Ellipsoidal Coordinate System created from :ref:`coordinate_system`.
         """
+        datum = Datum.from_user_input(datum).to_json_dict()
         geographic_crs_json = {
             "$schema": "https://proj.org/schemas/v0.2/projjson.schema.json",
             "type": "GeographicCRS",
             "name": name,
-            "datum": Datum.from_user_input(datum).to_json_dict(),
             "coordinate_system": CoordinateSystem.from_user_input(
                 ellipsoidal_cs or Ellipsoidal2DCS()
             ).to_json_dict(),
         }
+        if datum["type"] == "DatumEnsemble":
+            geographic_crs_json["datum_ensemble"] = datum
+        else:
+            geographic_crs_json["datum"] = datum
         super().__init__(geographic_crs_json)
 
 

diff --git a/test/crs/test_crs_cf.py b/test/crs/test_crs_cf.py
@@ -221,6 +221,7 @@ def test_cf_from_latlon():
         "grid_mapping_name": "latitude_longitude",
         "geographic_crs_name": "undefined",
         "reference_ellipsoid_name": "undefined",
+        "horizontal_datum_name": "undefined",
     }
     cf_dict = crs.to_cf()
     assert cf_dict.pop("crs_wkt").startswith("GEOGCRS[")
@@ -253,6 +254,7 @@ def test_cf_from_latlon__named():
         "reference_ellipsoid_name": "WGS 84",
         "longitude_of_prime_meridian": 0.0,
         "prime_meridian_name": "Greenwich",
+        "horizontal_datum_name": "World Geodetic System 1984 ensemble",
         "geographic_crs_name": "WGS 84",
         "grid_mapping_name": "latitude_longitude",
     }
@@ -360,7 +362,7 @@ def test_cf_rotated_latlon():
         "reference_ellipsoid_name": "WGS 84",
         "longitude_of_prime_meridian": 0.0,
         "prime_meridian_name": "Greenwich",
-        "horizontal_datum_name": "World Geodetic System 1984",
+        "horizontal_datum_name": "World Geodetic System 1984 ensemble",
         "grid_mapping_name": "rotated_latitude_longitude",
         "grid_north_pole_latitude": 32.5,
         "grid_north_pole_longitude": 170.0,
@@ -532,7 +534,7 @@ def test_cf_lambert_conformal_conic_1sp():
         "reference_ellipsoid_name": "WGS 84",
         "longitude_of_prime_meridian": 0.0,
         "prime_meridian_name": "Greenwich",
-        "horizontal_datum_name": "World Geodetic System 1984",
+        "horizontal_datum_name": "World Geodetic System 1984 ensemble",
         "grid_mapping_name": "lambert_conformal_conic",
         "longitude_of_central_meridian": 265.0,
         "false_easting": 0.0,
@@ -596,7 +598,7 @@ def test_cf_lambert_conformal_conic_2sp(standard_parallel):
         "reference_ellipsoid_name": "WGS 84",
         "longitude_of_prime_meridian": 0.0,
         "prime_meridian_name": "Greenwich",
-        "horizontal_datum_name": "World Geodetic System 1984",
+        "horizontal_datum_name": "World Geodetic System 1984 ensemble",
         "grid_mapping_name": "lambert_conformal_conic",
         "standard_parallel": (25.0, 30.0),
         "latitude_of_projection_origin": 25.0,
@@ -751,7 +753,7 @@ def test_geos_crs_sweep():
         "reference_ellipsoid_name": "WGS 84",
         "longitude_of_prime_meridian": 0.0,
         "prime_meridian_name": "Greenwich",
-        "horizontal_datum_name": "World Geodetic System 1984",
+        "horizontal_datum_name": "World Geodetic System 1984 ensemble",
         "grid_mapping_name": "geostationary",
         "sweep_angle_axis": "x",
         "perspective_point_height": 1.0,
@@ -799,7 +801,7 @@ def test_geos_crs_fixed_angle_axis():
         "reference_ellipsoid_name": "WGS 84",
         "longitude_of_prime_meridian": 0.0,
         "prime_meridian_name": "Greenwich",
-        "horizontal_datum_name": "World Geodetic System 1984",
+        "horizontal_datum_name": "World Geodetic System 1984 ensemble",
         "grid_mapping_name": "geostationary",
         "sweep_angle_axis": "x",
         "perspective_point_height": 1.0,
@@ -899,7 +901,7 @@ def test_mercator_b():
         "reference_ellipsoid_name": "WGS 84",
         "longitude_of_prime_meridian": 0.0,
         "prime_meridian_name": "Greenwich",
-        "horizontal_datum_name": "World Geodetic System 1984",
+        "horizontal_datum_name": "World Geodetic System 1984 ensemble",
         "grid_mapping_name": "mercator",
         "standard_parallel": 21.354,
         "longitude_of_projection_origin": 10.0,
@@ -1203,6 +1205,7 @@ def test_azimuthal_equidistant():
     assert cf_dict.pop("crs_wkt").startswith("PROJCRS[")
     assert cf_dict == expected_cf
     # test roundtrip
+    expected_cf["horizontal_datum_name"] = "World Geodetic System 1984 ensemble"
     _test_roundtrip(expected_cf, "PROJCRS[")
     # test coordinate system
     assert crs.cs_to_cf() == [
@@ -1230,7 +1233,7 @@ def test_lambert_azimuthal_equal_area():
         "reference_ellipsoid_name": "WGS 84",
         "longitude_of_prime_meridian": 0.0,
         "prime_meridian_name": "Greenwich",
-        "horizontal_datum_name": "World Geodetic System 1984",
+        "horizontal_datum_name": "World Geodetic System 1984 ensemble",
         "grid_mapping_name": "lambert_azimuthal_equal_area",
         "latitude_of_projection_origin": 1.0,
         "longitude_of_projection_origin": 2.0,
@@ -1270,7 +1273,7 @@ def test_lambert_cylindrical_equal_area():
         "reference_ellipsoid_name": "WGS 84",
         "longitude_of_prime_meridian": 0.0,
         "prime_meridian_name": "Greenwich",
-        "horizontal_datum_name": "World Geodetic System 1984",
+        "horizontal_datum_name": "World Geodetic System 1984 ensemble",
         "grid_mapping_name": "lambert_cylindrical_equal_area",
         "standard_parallel": 1.0,
         "longitude_of_central_meridian": 2.0,
@@ -1310,7 +1313,7 @@ def test_mercator_a():
         "reference_ellipsoid_name": "WGS 84",
         "longitude_of_prime_meridian": 0.0,
         "prime_meridian_name": "Greenwich",
-        "horizontal_datum_name": "World Geodetic System 1984",
+        "horizontal_datum_name": "World Geodetic System 1984 ensemble",
         "grid_mapping_name": "mercator",
         "standard_parallel": 0.0,
         "longitude_of_projection_origin": 2.0,
@@ -1351,7 +1354,7 @@ def test_orthographic():
         "reference_ellipsoid_name": "WGS 84",
         "longitude_of_prime_meridian": 0.0,
         "prime_meridian_name": "Greenwich",
-        "horizontal_datum_name": "World Geodetic System 1984",
+        "horizontal_datum_name": "World Geodetic System 1984 ensemble",
         "grid_mapping_name": "orthographic",
         "latitude_of_projection_origin": 1.0,
         "longitude_of_projection_origin": 2.0,
@@ -1391,7 +1394,7 @@ def test_polar_stereographic_a():
         "reference_ellipsoid_name": "WGS 84",
         "longitude_of_prime_meridian": 0.0,
         "prime_meridian_name": "Greenwich",
-        "horizontal_datum_name": "World Geodetic System 1984",
+        "horizontal_datum_name": "World Geodetic System 1984 ensemble",
         "grid_mapping_name": "polar_stereographic",
         "latitude_of_projection_origin": 90.0,
         "straight_vertical_longitude_from_pole": 1.0,
@@ -1432,7 +1435,7 @@ def test_polar_stereographic_b():
         "reference_ellipsoid_name": "WGS 84",
         "longitude_of_prime_meridian": 0.0,
         "prime_meridian_name": "Greenwich",
-        "horizontal_datum_name": "World Geodetic System 1984",
+        "horizontal_datum_name": "World Geodetic System 1984 ensemble",
         "grid_mapping_name": "polar_stereographic",
         "standard_parallel": 0.0,
         "straight_vertical_longitude_from_pole": 1.0,
@@ -1472,7 +1475,7 @@ def test_stereographic():
         "reference_ellipsoid_name": "WGS 84",
         "longitude_of_prime_meridian": 0.0,
         "prime_meridian_name": "Greenwich",
-        "horizontal_datum_name": "World Geodetic System 1984",
+        "horizontal_datum_name": "World Geodetic System 1984 ensemble",
         "grid_mapping_name": "stereographic",
         "latitude_of_projection_origin": 0.0,
         "longitude_of_projection_origin": 1.0,
@@ -1513,7 +1516,7 @@ def test_sinusoidal():
         "reference_ellipsoid_name": "WGS 84",
         "longitude_of_prime_meridian": 0.0,
         "prime_meridian_name": "Greenwich",
-        "horizontal_datum_name": "World Geodetic System 1984",
+        "horizontal_datum_name": "World Geodetic System 1984 ensemble",
         "grid_mapping_name": "sinusoidal",
         "longitude_of_projection_origin": 0.0,
         "false_easting": 1.0,
@@ -1552,7 +1555,7 @@ def test_vertical_perspective():
         "reference_ellipsoid_name": "WGS 84",
         "longitude_of_prime_meridian": 0.0,
         "prime_meridian_name": "Greenwich",
-        "horizontal_datum_name": "World Geodetic System 1984",
+        "horizontal_datum_name": "World Geodetic System 1984 ensemble",
         "grid_mapping_name": "vertical_perspective",
         "perspective_point_height": 50.0,
         "latitude_of_projection_origin": 0.0,
@@ -1889,7 +1892,7 @@ def test_3d_ellipsoidal_cs_depth():
             "name": "WGS 84 (geographic 3D)",
             "datum": {
                 "type": "GeodeticReferenceFrame",
-                "name": "World Geodetic System 1984",
+                "name": "World Geodetic System 1984 ensemble",
                 "ellipsoid": {
                     "name": "WGS 84",
                     "semi_major_axis": 6378137,