BUG: Add datum ensemble support to GeographicCRS

Author: snowman2

docs/history.rst

@@ -13,7 +13,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 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
 -----

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
 
 

pyproj/crs/crs.py

@@ -718,8 +718,6 @@ def to_cf(
             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)
 
 

test/crs/test_crs_cf.py

@@ -362,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,
@@ -534,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,
@@ -598,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,
@@ -753,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,
@@ -801,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,
@@ -901,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,
@@ -1205,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() == [
@@ -1232,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,
@@ -1272,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,
@@ -1312,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,
@@ -1353,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,
@@ -1393,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,
@@ -1434,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,
@@ -1474,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,
@@ -1515,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,
@@ -1554,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,
@@ -1891,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,