From 522db5a1f046172984c9385c5550fe68106caa4a Mon Sep 17 00:00:00 2001 From: snowman2 Date: Fri, 12 Jun 2020 20:48:29 -0500 Subject: [PATCH] ENH: Add support for coordinate systems with CRS using CF conventions --- docs/build_crs_cf.rst | 193 ++++++++ docs/history.rst | 1 + docs/index.rst | 7 +- pyproj/_crs.pyi | 1 + pyproj/_crs.pyx | 79 ++++ pyproj/crs/crs.py | 93 +++- test/crs/test_crs_cf.py | 625 ++++++++++++++++++++++++- test/crs/test_crs_coordinate_system.py | 101 ++++ 8 files changed, 1086 insertions(+), 14 deletions(-) create mode 100644 docs/build_crs_cf.rst diff --git a/docs/build_crs_cf.rst b/docs/build_crs_cf.rst new file mode 100644 index 000000000..4ee511fea --- /dev/null +++ b/docs/build_crs_cf.rst @@ -0,0 +1,193 @@ +.. _build_crs_cf: + +Managing CRS to and from CF +============================ + +http://cfconventions.org/cf-conventions/cf-conventions.html + +Exporting CRS to CF +-------------------- +When exporting a CRS to the Climate and Forecast (CF) conventions, +you need both the grid mapping as well as the coordinate system. +If you don't use the coordinate system, then you will lose the units +of your projection. + + +In this example, this is the CRS we will use: + +.. code-block:: python + + from pyproj import CRS + + crs = CRS("EPSG:4326") + + +To get the grid mapping you use :meth:`pyproj.crs.CRS.to_cf`: + +.. versionadded:: 2.2.0 + +.. code-block:: python + + cf_grid_mapping = crs.to_cf() + + +Contents of `cf_grid_mapping`:: + + + {'crs_wkt': 'GEOGCRS["WGS 84",DATUM["World Geodetic System ' + ....,ID["EPSG",4326]]', + 'geographic_crs_name': 'WGS 84', + 'grid_mapping_name': 'latitude_longitude', + 'inverse_flattening': 298.257223563, + 'longitude_of_prime_meridian': 0.0, + 'prime_meridian_name': 'Greenwich', + 'reference_ellipsoid_name': 'WGS 84', + 'semi_major_axis': 6378137.0, + 'semi_minor_axis': 6356752.314245179} + + +To get the coordinate system, you use :meth:`pyproj.crs.CRS.cs_to_cf`: + +.. versionadded:: 3.0.0 + +.. code-block:: python + + cf_coordinate_system = crs.cs_to_cf() + + +Contents of `cf_coordinate_system`:: + + [{'long_name': 'geodetic latitude coordinate', + 'standard_name': 'geodetic latitude', + 'unit': 'degrees_north'}, + {'long_name': 'geodetic longitude coordinate', + 'standard_name': 'geodetic longitude', + 'unit': 'degrees_east'}] + + +Importing CRS from CF +---------------------- + +When importing a CRS from the Climate and Forecast (CF) conventions, +you need both the grid mapping as well as the coordinate system. +If you don't use the coordinate system, then you will lose the units +of your projection. + +.. note:: If the CF `crs_wkt` attribute is available, the coordinate system is + inside of the WKT and can be used to create the CRS in a single step. + +.. warning:: If building from grid mapping, be mindful of the axis order. https://github.com/cf-convention/cf-conventions/pull/224 + + +Build the CRS from CF grid mapping: +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +In this example, this is the grid mapping and coordinate system we will use:: + + variables: + double x(x) ; + x:standard_name = "projection_x_coordinate" ; + x:long_name = "Easting" ; + x:units = "m" ; + double y(y) ; + y:standard_name = "projection_y_coordinate" ; + y:long_name = "Northing" ; + y:units = "m" ; + int crsOSGB ; + crsOSGB:grid_mapping_name = "transverse_mercator"; + crsOSGB:semi_major_axis = 6377563.396 ; + crsOSGB:inverse_flattening = 299.3249646 ; + crsOSGB:longitude_of_prime_meridian = 0.0 ; + crsOSGB:latitude_of_projection_origin = 49.0 ; + crsOSGB:longitude_of_central_meridian = -2.0 ; + crsOSGB:scale_factor_at_central_meridian = 0.9996012717 ; + crsOSGB:false_easting = 400000.0 ; + crsOSGB:false_northing = -100000.0 ; + +.. note:: If the units are meters as in this example, + then no further changes are necessary. + +.. code-block:: python + + from pyproj import CRS + + crs = CRS.from_cf( + { + "grid_mapping_name": "transverse_mercator", + "semi_major_axis": 6377563.396, + "inverse_flattening": 299.3249646, + "longitude_of_prime_meridian": 0.0, + "latitude_of_projection_origin": 49.0, + "longitude_of_central_meridian": -2.0, + "scale_factor_at_central_meridian": 0.9996012717, + "false_easting": 400000.0, + "false_northing": -100000.0, + } + ) + + +Modify the CRS with coordinate system: +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. versionadded:: 3.0.0 + +.. note:: If the CF `crs_wkt` attribute is available, the coordinate system is + inside of the WKT and can be used to create the CRS in a single step. + +.. warning:: Be mindful of the axis order. https://github.com/cf-convention/cf-conventions/pull/224 + + +In this example, assume everything is the same as above. +However, the units are instead `US_Survey_Foot`:: + + variables: + double x(x) ; + x:standard_name = "projection_x_coordinate" ; + x:long_name = "Easting" ; + x:units = "US_Survey_Foot" ; + double y(y) ; + y:standard_name = "projection_y_coordinate" ; + y:long_name = "Northing" ; + y:units = "US_Survey_Foot" ; + ... + + +In this case, you will need to get the unit conversion factor: + +https://github.com/SciTools/cf-units + + +.. code-block:: python + + from cf_units import Unit + from pyproj import CRS + + cf_unit = Unit("US_Survey_Foot") + unit = { + "type": "LinearUnit", + "name": "US Survey Foot", + "conversion_factor": cf_unit.convert(1, "m"), + } + cartesian_cs = { + "$schema": "https://proj.org/schemas/v0.2/projjson.schema.json", + "type": "CoordinateSystem", + "subtype": "Cartesian", + "axis": [ + {"name": "Easting", "abbreviation": "E", "direction": "east", "unit": unit}, + {"name": "Northing", "abbreviation": "N", "direction": "north", "unit": unit}, + ], + } + crs = CRS.from_cf( + { + "grid_mapping_name": "transverse_mercator", + "semi_major_axis": 6377563.396, + "inverse_flattening": 299.3249646, + "longitude_of_prime_meridian": 0.0, + "latitude_of_projection_origin": 49.0, + "longitude_of_central_meridian": -2.0, + "scale_factor_at_central_meridian": 0.9996012717, + "false_easting": 400000.0, + "false_northing": -100000.0, + }, + cartesian_cs=cartesian_cs, + ) diff --git a/docs/history.rst b/docs/history.rst index 8af43f302..a9742221d 100644 --- a/docs/history.rst +++ b/docs/history.rst @@ -15,6 +15,7 @@ Change Log * ENH: Use from_user_input in __eq__ when comparing CRS sub-classes (i.e. PrimeMeridian, Datum, Ellipsoid, etc.) (issue #606) * ENH: Added :meth:`pyproj.transformer.TransformerGroup.download_grids` (pull #643) * ENH: Added transformation grid sync API/CLI (issue #572) +* ENH: Add support for coordinate systems with CRS using CF conventions (issue #536) 2.6.1 ~~~~~ diff --git a/docs/index.rst b/docs/index.rst index 3ab5b196e..256bdcb01 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -18,12 +18,13 @@ GitHub Repository: https://github.com/pyproj4/pyproj examples transformation_grids gotchas - advanced_examples - build_crs api/index cli - optimize_transformations + advanced_examples + build_crs + build_crs_cf crs_compatibility + optimize_transformations history past_versions diff --git a/pyproj/_crs.pyi b/pyproj/_crs.pyi index 1aa7e680d..1d2825f08 100644 --- a/pyproj/_crs.pyi +++ b/pyproj/_crs.pyi @@ -115,6 +115,7 @@ class CoordinateSystem(_CRSParts): def from_json_dict(coordinate_system_dict: dict) -> "CoordinateSystem": ... @staticmethod def from_json(coordinate_system_json_str: str) -> "CoordinateSystem": ... + def to_cf(self, rotated_pole: bool = False) -> List[dict]: ... class Param: name: str diff --git a/pyproj/_crs.pyx b/pyproj/_crs.pyx index 40be4f930..e9972d7f5 100644 --- a/pyproj/_crs.pyx +++ b/pyproj/_crs.pyx @@ -710,6 +710,85 @@ cdef class CoordinateSystem(_CRSParts): _load_proj_json(coordinate_system_json_str) ) + def to_cf(self, rotated_pole=False): + """ + .. versionadded:: 3.0.0 + + This converts a :obj:`pyproj.crs.CoordinateSystem` axis + to a list of Climate and Forecast (CF) Version 1.8 dicts. + + Parameters + ---------- + rotated_pole: bool, optional + If True, the geographic coordinates are on a rotated pole grid. + This corresponds to the rotated_latitude_longitude grid_mapping_name. + Default is False. + + Returns + ------- + List[dict]: + CF-1.8 version of the CoordinateSystem. + """ + axis_list = self.to_json_dict()["axis"] + cf_params = [] + def get_linear_unit(axis): + try: + return f'{axis["unit"]["conversion_factor"]} metre' + except TypeError: + return axis["unit"] + + if self.name == "cartesian": + for axis in axis_list: + if axis["name"].lower() == "easting": + cf_axis = "X" + else: + cf_axis = "Y" + cf_params.append(dict( + axis=cf_axis, + long_name=axis["name"], + standard_name=f"projection_{cf_axis.lower()}_coordinate", + unit=get_linear_unit(axis), + )) + elif self.name == "ellipsoidal": + for axis in axis_list: + if axis["abbreviation"].upper() == "H": + cf_params.append(dict( + standard_name="height_above_reference_ellipsoid", + long_name=axis["name"], + unit=axis["unit"], + positive=axis["direction"], + axis="Z", + )) + else: + name = axis["name"].lower() + if rotated_pole: + cf_params.append(dict( + standard_name=f"grid_{name}", + long_name=f"{name} in rotated pole grid", + unit="degrees", + )) + else: + cf_params.append(dict( + standard_name=name, + long_name=f"{name} coordinate", + unit=f'degrees_{axis["direction"]}', + )) + elif self.name == "vertical": + for axis in axis_list: + if axis["abbreviation"].upper() == "H": + standard_name = "height" + else: + standard_name = "depth" + cf_params.append(dict( + standard_name=standard_name, + long_name=axis["name"], + unit=get_linear_unit(axis), + positive=axis["direction"], + axis="Z", + )) + + return cf_params + cdef class Ellipsoid(_CRSParts): """ diff --git a/pyproj/crs/crs.py b/pyproj/crs/crs.py index 87283bf44..28927731a 100644 --- a/pyproj/crs/crs.py +++ b/pyproj/crs/crs.py @@ -562,9 +562,7 @@ def to_cf( This converts a :obj:`pyproj.crs.CRS` object to a Climate and Forecast (CF) Grid Mapping Version 1.8 dict. - .. warning:: The full projection will be stored in the - crs_wkt attribute. However, other parameters may be lost - if a mapping to the CF parameter is not found. + :ref:`build_crs_cf` Parameters ---------- @@ -672,21 +670,39 @@ def to_cf( return cf_dict @staticmethod - def from_cf(in_cf: dict, errcheck=False) -> "CRS": + def from_cf( + in_cf: dict, + ellipsoidal_cs: Any = None, + cartesian_cs: Any = None, + vertical_cs: Any = None, + errcheck=False, + ) -> "CRS": """ .. versionadded:: 2.2.0 + .. versionadded:: 3.0.0 ellipsoidal_cs, cartesian_cs, vertical_cs + This converts a Climate and Forecast (CF) Grid Mapping Version 1.8 dict to a :obj:`pyproj.crs.CRS` object. - .. warning:: Parameters may be lost if a mapping - from the CF parameter is not found. For best results - store the WKT of the projection in the crs_wkt attribute. + :ref:`build_crs_cf` Parameters ---------- in_cf: dict CF version of the projection. + 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`. + cartesian_cs: Any, optional + Input to create a Cartesian Coordinate System. + Anything accepted by :meth:`pyproj.crs.CoordinateSystem.from_user_input` + or :class:`pyproj.crs.coordinate_system.Cartesian2DCS`. + vertical_cs: Any, optional + Input to create a Vertical Coordinate System accepted by + :meth:`pyproj.crs.CoordinateSystem.from_user_input` + or :class:`pyproj.crs.coordinate_system.VerticalCS` errcheck: bool, optional This parameter is for backwards compatibility with the old version. It currently does nothing when True or False. @@ -719,17 +735,27 @@ def from_cf(in_cf: dict, errcheck=False) -> "CRS": geographic_crs_name = in_cf.get("geographic_crs_name") if datum: geographic_crs = GeographicCRS( - name=geographic_crs_name or "undefined", datum=datum, + name=geographic_crs_name or "undefined", + datum=datum, + ellipsoidal_cs=ellipsoidal_cs, ) # type: CRS elif geographic_crs_name and geographic_crs_name not in unknown_names: geographic_crs = CRS(geographic_crs_name) + if ellipsoidal_cs is not None: + geographic_crs_json = geographic_crs.to_json_dict() + geographic_crs_json[ + "coordinate_system" + ] = CoordinateSystem.from_user_input(ellipsoidal_cs).to_json_dict() + geographic_crs = CRS(geographic_crs_json) else: - geographic_crs = GeographicCRS() + geographic_crs = GeographicCRS(ellipsoidal_cs=ellipsoidal_cs) if grid_mapping_name == "latitude_longitude": return geographic_crs if geographic_conversion_method is not None: return DerivedGeographicCRS( - base_crs=geographic_crs, conversion=geographic_conversion_method(in_cf), + base_crs=geographic_crs, + conversion=geographic_conversion_method(in_cf), + ellipsoidal_cs=ellipsoidal_cs, ) # build projected CRS @@ -741,6 +767,7 @@ def from_cf(in_cf: dict, errcheck=False) -> "CRS": name=in_cf.get("projected_crs_name", "undefined"), conversion=conversion_method(in_cf), geodetic_crs=geographic_crs, + cartesian_cs=cartesian_cs, ) # build bound CRS if exists @@ -761,6 +788,7 @@ def from_cf(in_cf: dict, errcheck=False) -> "CRS": name="undefined", datum=in_cf["geopotential_datum_name"], geoid_model=in_cf.get("geoid_name"), + vertical_cs=vertical_cs, ) # build compound CRS @@ -768,6 +796,51 @@ def from_cf(in_cf: dict, errcheck=False) -> "CRS": name="undefined", components=[bound_crs or projected_crs, vertical_crs] ) + def cs_to_cf(self) -> List[dict]: + """ + .. versionadded:: 3.0.0 + + This converts all coordinate systems (cs) in the CRS + to a list of Climate and Forecast (CF) Version 1.8 dicts. + + :ref:`build_crs_cf` + + Returns + ------- + List[dict]: + CF-1.8 version of the coordinate systems. + """ + cf_axis_list = [] + + def rotated_pole(crs): + try: + return ( + crs.coordinate_operation + and crs.coordinate_operation.method_name.lower() + in _INVERSE_GEOGRAPHIC_GRID_MAPPING_NAME_MAP + ) + except KeyError: + return False + + if self.coordinate_system: + cf_axis_list.extend( + self.coordinate_system.to_cf(rotated_pole=rotated_pole(self)) + ) + elif self.is_bound and self.source_crs and self.source_crs.coordinate_system: + cf_axis_list.extend( + self.source_crs.coordinate_system.to_cf( + rotated_pole=rotated_pole(self.source_crs) + ) + ) + else: + for sub_crs in self.sub_crs_list: + if not sub_crs.coordinate_system: + continue + cf_axis_list.extend( + sub_crs.coordinate_system.to_cf(rotated_pole=rotated_pole(sub_crs)) + ) + return cf_axis_list + def is_exact_same(self, other: Any, ignore_axis_order: bool = False) -> bool: """ Check if the CRS objects are the exact same. diff --git a/test/crs/test_crs_cf.py b/test/crs/test_crs_cf.py index e461d17b6..7679ce37a 100644 --- a/test/crs/test_crs_cf.py +++ b/test/crs/test_crs_cf.py @@ -79,6 +79,21 @@ def test_to_cf_transverse_mercator(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "BOUNDCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] with pytest.warns(UserWarning): assert crs.to_dict() == { @@ -140,6 +155,21 @@ def test_from_cf_transverse_mercator(towgs84_test): # test roundtrip expected_cf["towgs84"] = _try_list_if_string(towgs84_test) _test_roundtrip(expected_cf, "BOUNDCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_cf_from_latlon(): @@ -165,6 +195,19 @@ def test_cf_from_latlon(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "GEOGCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "standard_name": "longitude", + "long_name": "longitude coordinate", + "unit": "degrees_east", + }, + { + "standard_name": "latitude", + "long_name": "latitude coordinate", + "unit": "degrees_north", + }, + ] def test_cf_from_latlon__named(): @@ -210,6 +253,21 @@ def test_cf_from_utm(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_cf_from_utm__nad83(): @@ -236,6 +294,21 @@ def test_cf_from_utm__nad83(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_cf_rotated_latlon(): @@ -265,6 +338,19 @@ def test_cf_rotated_latlon(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "GEOGCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "standard_name": "grid_longitude", + "long_name": "longitude in rotated pole grid", + "unit": "degrees", + }, + { + "standard_name": "grid_latitude", + "long_name": "latitude in rotated pole grid", + "unit": "degrees", + }, + ] with pytest.warns(UserWarning): proj_dict = crs.to_dict() assert proj_dict == { @@ -332,6 +418,21 @@ def test_cf_lambert_conformal_conic_1sp(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] with pytest.warns(UserWarning): proj_dict = crs.to_dict() @@ -382,7 +483,21 @@ def test_cf_lambert_conformal_conic_2sp(standard_parallel): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") - + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] with pytest.warns(UserWarning): proj_dict = crs.to_dict() assert proj_dict == { @@ -435,6 +550,21 @@ def test_oblique_mercator(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] with pytest.warns(UserWarning): assert crs.to_dict() == { "proj": "omerc", @@ -509,6 +639,21 @@ def test_geos_crs_sweep(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_geos_crs_fixed_angle_axis(): @@ -542,6 +687,21 @@ def test_geos_crs_fixed_angle_axis(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_geos_proj_string(): @@ -569,6 +729,21 @@ def test_geos_proj_string(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_ob_tran_not_rotated_latlon(): @@ -622,6 +797,21 @@ def test_mercator_b(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_osgb_1936(): @@ -651,6 +841,21 @@ def test_osgb_1936(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_export_compound_crs(): @@ -678,6 +883,28 @@ def test_export_compound_crs(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "COMPOUNDCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "standard_name": "height", + "long_name": "Gravity-related height", + "unit": "metre", + "positive": "up", + "axis": "Z", + }, + ] def test_geoid_model_name(): @@ -756,6 +983,27 @@ def test_geoid_model_name(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "COMPOUNDCRS[") + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + { + "standard_name": "height", + "long_name": "Gravity-related height", + "unit": "metre", + "positive": "up", + "axis": "Z", + }, + ] def test_albers_conical_equal_area(): @@ -782,6 +1030,21 @@ def test_albers_conical_equal_area(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_azimuthal_equidistant(): @@ -807,6 +1070,21 @@ def test_azimuthal_equidistant(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_lambert_azimuthal_equal_area(): @@ -832,6 +1110,21 @@ def test_lambert_azimuthal_equal_area(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_lambert_cylindrical_equal_area(): @@ -857,6 +1150,21 @@ def test_lambert_cylindrical_equal_area(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_mercator_a(): @@ -883,6 +1191,21 @@ def test_mercator_a(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_orthographic(): @@ -908,6 +1231,21 @@ def test_orthographic(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_polar_stereographic_a(): @@ -934,6 +1272,21 @@ def test_polar_stereographic_a(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_polar_stereographic_b(): @@ -959,6 +1312,21 @@ def test_polar_stereographic_b(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_stereographic(): @@ -985,6 +1353,21 @@ def test_stereographic(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_sinusoidal(): @@ -1009,6 +1392,21 @@ def test_sinusoidal(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_vertical_perspective(): @@ -1035,6 +1433,21 @@ def test_vertical_perspective(): assert cf_dict == expected_cf # test roundtrip _test_roundtrip(expected_cf, "PROJCRS[") + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + ] def test_build_custom_datum(): @@ -1095,3 +1508,213 @@ def test_build_custom_datum__default_ellipsoid(): assert crs.ellipsoid.name == "WGS 84" assert crs.prime_meridian.name == "Paris" assert str(crs.prime_meridian.longitude).startswith("2.") + + +def test_cartesian_cs(): + unit = {"type": "LinearUnit", "name": "US Survey Foot", "conversion_factor": 0.3048} + cartesian_cs = { + "$schema": "https://proj.org/schemas/v0.2/projjson.schema.json", + "type": "CoordinateSystem", + "subtype": "Cartesian", + "axis": [ + {"name": "Easting", "abbreviation": "E", "direction": "east", "unit": unit}, + { + "name": "Northing", + "abbreviation": "N", + "direction": "north", + "unit": unit, + }, + ], + } + crs = CRS.from_cf( + { + "grid_mapping_name": "transverse_mercator", + "semi_major_axis": 6377563.396, + "inverse_flattening": 299.3249646, + "longitude_of_prime_meridian": 0.0, + "latitude_of_projection_origin": 49.0, + "longitude_of_central_meridian": -2.0, + "scale_factor_at_central_meridian": 0.9996012717, + "false_easting": 400000.0, + "false_northing": -100000.0, + }, + cartesian_cs=cartesian_cs, + ) + assert crs.coordinate_system.to_json_dict() == cartesian_cs + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "0.3048 metre", + }, + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "0.3048 metre", + }, + ] + + +def test_ellipsoidal_cs(): + ellipsoidal_cs = { + "$schema": "https://proj.org/schemas/v0.2/projjson.schema.json", + "type": "CoordinateSystem", + "subtype": "ellipsoidal", + "axis": [ + { + "name": "Latitude", + "abbreviation": "lat", + "direction": "north", + "unit": "degree", + }, + { + "name": "Longitude", + "abbreviation": "lon", + "direction": "east", + "unit": "degree", + }, + ], + } + crs = CRS.from_cf( + dict( + grid_mapping_name="latitude_longitude", + semi_major_axis=6378137.0, + inverse_flattening=298.257223, + ), + ellipsoidal_cs=ellipsoidal_cs, + ) + assert crs.coordinate_system.to_json_dict() == ellipsoidal_cs + # test coordinate system + assert crs.cs_to_cf() == [ + { + "standard_name": "latitude", + "long_name": "latitude coordinate", + "unit": "degrees_north", + }, + { + "standard_name": "longitude", + "long_name": "longitude coordinate", + "unit": "degrees_east", + }, + ] + + +def test_ellipsoidal_cs__from_name(): + ellipsoidal_cs = { + "$schema": "https://proj.org/schemas/v0.2/projjson.schema.json", + "type": "CoordinateSystem", + "subtype": "ellipsoidal", + "axis": [ + { + "name": "Longitude", + "abbreviation": "lon", + "direction": "east", + "unit": "degree", + }, + { + "name": "Latitude", + "abbreviation": "lat", + "direction": "north", + "unit": "degree", + }, + ], + } + crs = CRS.from_cf( + dict(grid_mapping_name="latitude_longitude", geographic_crs_name="WGS 84",), + ellipsoidal_cs=ellipsoidal_cs, + ) + assert crs.coordinate_system.to_json_dict() == ellipsoidal_cs + # test coordinate system + assert crs.cs_to_cf() == [ + { + "standard_name": "longitude", + "long_name": "longitude coordinate", + "unit": "degrees_east", + }, + { + "standard_name": "latitude", + "long_name": "latitude coordinate", + "unit": "degrees_north", + }, + ] + + +def test_export_compound_crs_cs(): + unit = {"type": "LinearUnit", "name": "US Survey Foot", "conversion_factor": 0.3048} + cartesian_cs = { + "$schema": "https://proj.org/schemas/v0.2/projjson.schema.json", + "type": "CoordinateSystem", + "subtype": "Cartesian", + "axis": [ + { + "name": "Northing", + "abbreviation": "N", + "direction": "north", + "unit": unit, + }, + {"name": "Easting", "abbreviation": "E", "direction": "east", "unit": unit}, + ], + } + vertical_cs = { + "$schema": "https://proj.org/schemas/v0.2/projjson.schema.json", + "type": "CoordinateSystem", + "subtype": "vertical", + "axis": [ + { + "name": "Gravity-related height", + "abbreviation": "H", + "direction": "up", + "unit": unit, + } + ], + } + + crs = CRS.from_cf( + { + "semi_major_axis": 6378388.0, + "semi_minor_axis": 6356911.9461279465, + "inverse_flattening": 297.0, + "reference_ellipsoid_name": "International 1924", + "longitude_of_prime_meridian": 0.0, + "prime_meridian_name": "Greenwich", + "geographic_crs_name": "KKJ", + "horizontal_datum_name": "Kartastokoordinaattijarjestelma (1966)", + "projected_crs_name": "KKJ / Finland Uniform Coordinate System", + "grid_mapping_name": "transverse_mercator", + "latitude_of_projection_origin": 0.0, + "longitude_of_central_meridian": 27.0, + "false_easting": 3500000.0, + "false_northing": 0.0, + "scale_factor_at_central_meridian": 1.0, + "geopotential_datum_name": "Helsinki 1960", + }, + cartesian_cs=cartesian_cs, + vertical_cs=vertical_cs, + ) + assert crs.sub_crs_list[0].coordinate_system.to_json_dict() == cartesian_cs + assert crs.sub_crs_list[1].coordinate_system.to_json_dict() == vertical_cs + # test coordinate system + assert crs.cs_to_cf() == [ + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "0.3048 metre", + }, + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "0.3048 metre", + }, + { + "standard_name": "height", + "long_name": "Gravity-related height", + "unit": "0.3048 metre", + "positive": "up", + "axis": "Z", + }, + ] diff --git a/test/crs/test_crs_coordinate_system.py b/test/crs/test_crs_coordinate_system.py index 9ef2ceb4c..96ee2f415 100644 --- a/test/crs/test_crs_coordinate_system.py +++ b/test/crs/test_crs_coordinate_system.py @@ -163,3 +163,104 @@ def test_ellipsoidal_3d_cs(axis, name_0, direction_0, name_1, direction_1): assert vcs.axis_list[2].direction == "up" assert vcs.axis_list[2].name == "Ellipsoidal height" assert vcs.axis_list[2].unit_name == "metre" + + +def test_ellipsoidal2dcs_to_cf(): + ecs = Ellipsoidal2DCS(axis=Ellipsoidal2DCSAxis.LATITUDE_LONGITUDE) + assert ecs.to_cf() == [ + { + "standard_name": "latitude", + "long_name": "latitude coordinate", + "unit": "degrees_north", + }, + { + "standard_name": "longitude", + "long_name": "longitude coordinate", + "unit": "degrees_east", + }, + ] + + +def test_ellipsoidal3dcs_to_cf(): + ecs = Ellipsoidal3DCS(axis=Ellipsoidal3DCSAxis.LONGITUDE_LATITUDE_HEIGHT) + assert ecs.to_cf() == [ + { + "standard_name": "longitude", + "long_name": "longitude coordinate", + "unit": "degrees_east", + }, + { + "standard_name": "latitude", + "long_name": "latitude coordinate", + "unit": "degrees_north", + }, + { + "standard_name": "height_above_reference_ellipsoid", + "long_name": "Ellipsoidal height", + "unit": "metre", + "positive": "up", + "axis": "Z", + }, + ] + + +def test_cartesian2dcs_ft_to_cf(): + csft = Cartesian2DCS(axis=Cartesian2DCSAxis.NORTHING_EASTING_FT) + csft.to_cf() == [ + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "0.3048 metre", + }, + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "0.3048 metre", + }, + ] + + +def test_cartesian2dcs_to_cf(): + csm = Cartesian2DCS(axis=Cartesian2DCSAxis.EASTING_NORTHING_FT) + csm.to_cf() == [ + { + "axis": "Y", + "long_name": "Northing", + "standard_name": "projection_y_coordinate", + "unit": "metre", + }, + { + "axis": "X", + "long_name": "Easting", + "standard_name": "projection_x_coordinate", + "unit": "metre", + }, + ] + + +def test_verticalcs_depth_to_cf(): + vcs = VerticalCS(axis=VerticalCSAxis.DEPTH) + vcs.to_cf() == [ + { + "standard_name": "depth", + "long_name": "Depth", + "unit": "metre", + "positive": "down", + "axis": "Z", + } + ] + + +def test_verticalcs_height_to_cf(): + vcs = VerticalCS(axis=VerticalCSAxis.GRAVITY_HEIGHT_US_FT) + vcs.to_cf() == [ + { + "standard_name": "height", + "long_name": "Gravity-related height", + "unit": "0.304800609601219 metre", + "positive": "up", + "axis": "Z", + } + ]