diff --git a/pyproj/cf1x8.py b/pyproj/cf1x8.py new file mode 100644 index 000000000..1843cb35a --- /dev/null +++ b/pyproj/cf1x8.py @@ -0,0 +1,83 @@ +""" +This module contains mappings necessary to convert from +a CRS to a CF-1.8 compliant projection. + +http://cfconventions.org/cf-conventions/cf-conventions.html#appendix-grid-mappings + + +It is not complete and does not support all projections. + +CF PARAMS (NOT SURE OF THE MAPPING): +------------------------------------ +geoid_name (OGC WKT VERT_DATUM - ex GEOID12B) +geopotential_datum_name (OGC WKT VERT_DATUM - ex NAVD88) + +""" + +GRID_MAPPING_NAME_MAP = { + "albers_conical_equal_area": "aea", + "azimuthal_equidistant": "aeqd", + "geostationary": "geos", + "lambert_azimuthal_equal_area": "laea", + "lambert_conformal_conic": "lcc", + "lambert_cylindrical_equal_area": "cea", + "mercator": "merc", + "oblique_mercator": "omerc", + "orthographic": "ortho", + "polar_stereographic": "stere", + "sinusoidal": "sinu", + "stereographic": "stere", + "transverse_mercator": "tmerc", + "vertical_perspective": "nsper", + "rotated_latitude_longitude": "ob_tran", + "latitude_longitude": "latlon", +} + + +INVERSE_GRID_MAPPING_NAME_MAP = { + value: key for key, value in GRID_MAPPING_NAME_MAP.items() +} + +PROJ_PARAM_MAP = { + "azimuth_of_central_line": "gamma", + "earth_radius": "R", + "fase_easting": "x_0", + "fase_northing": "y_0", + "latitude_of_projection_origin": "lat_0", + "north_pole_grid_longitude": "lon_0", + "straight_vertical_longitude_from_pole": "lon_0", + "longitude_of_central_meridian": "lon_0", + "longitude_of_projection_origin": "lon_0", + "horizontal_datum_name": "datum", + "reference_ellipsoid_name": "ellps", + "towgs84": "towgs84", + "prime_meridian_name": "pm", + "scale_factor_at_central_meridian": "k_0", + "scale_factor_at_projection_origin": "k_0", + "unit": "units", + "perspective_point_height": "h", + "grid_north_pole_longitude": "o_lon_p", + "grid_north_pole_latitude": "o_lat_p", + "semi_major_axis": "a", + "semi_minor_axis": "b", + "inverse_flattening": "rf", +} + + +INVERSE_PROJ_PARAM_MAP = {value: key for key, value in PROJ_PARAM_MAP.items()} + + +LON_0_MAP = { + "DEFAULT": "longitude_of_projection_origin", + "rotated_latitude_longitude": "north_pole_grid_longitude", + "polar_stereographic": "straight_vertical_longitude_from_pole", + "transverse_mercator": "longitude_of_central_meridian", + "lambert_cylindrical_equal_area": "longitude_of_central_meridian", + "lambert_conformal_conic": "longitude_of_central_meridian", + "albers_conical_equal_area": "longitude_of_central_meridian", +} + +K_0_MAP = { + "DEFAULT": "scale_factor_at_projection_origin", + "transverse_mercator": "scale_factor_at_central_meridian", +} diff --git a/pyproj/crs.py b/pyproj/crs.py index 8ace7b85c..8f8bb7daf 100644 --- a/pyproj/crs.py +++ b/pyproj/crs.py @@ -22,8 +22,17 @@ __all__ = ["CRS", "is_wkt"] import json +import warnings from pyproj._crs import _CRS, is_wkt +from pyproj.cf1x8 import ( + GRID_MAPPING_NAME_MAP, + INVERSE_GRID_MAPPING_NAME_MAP, + INVERSE_PROJ_PARAM_MAP, + K_0_MAP, + LON_0_MAP, + PROJ_PARAM_MAP, +) from pyproj.compat import string_types from pyproj.exceptions import CRSError from pyproj.geod import Geod @@ -33,6 +42,9 @@ def _dict2string(projparams): # convert a dict to a proj4 string. pjargs = [] for key, value in projparams.items(): + # the towgs84 as list + if isinstance(value, list): + value = ",".join([str(val) for val in value]) # issue 183 (+ no_rot) if value is None or value is True: pjargs.append("+" + key + " ") @@ -49,8 +61,8 @@ class CRS(_CRS): The functionality is based on other fantastic projects: - * https://github.com/opendatacube/datacube-core/blob/83bae20d2a2469a6417097168fd4ede37fd2abe5/datacube/utils/geometry/_base.py - * https://github.com/mapbox/rasterio/blob/c13f0943b95c0eaa36ff3f620bd91107aa67b381/rasterio/_crs.pyx + * `rasterio `_ + * `opendatacube `_ """ @@ -245,6 +257,179 @@ def get_geod(self): in_kwargs["b"] = self.ellipsoid.semi_minor_metre return Geod(**in_kwargs) + def to_proj4_dict(self): + """ + Converts the PROJ string to a dict. + + Returns + ------- + dict: PROJ params in dict format. + + """ + + def parse(val): + if val.lower() == "true": + return True + elif val.lower() == "false": + return False + try: + return int(val) + except ValueError: + pass + try: + return float(val) + except ValueError: + pass + val_split = val.split(",") + if len(val_split) > 1: + val = [float(sval.strip()) for sval in val_split] + return val + + items = map( + lambda kv: len(kv) == 2 and (kv[0], parse(kv[1])) or (kv[0], None), + ( + part.lstrip("+").split("=", 1) + for part in self.to_proj4().strip().split() + ), + ) + + return {key: value for key, value in items if value is not False} + + def to_cf(self, wkt_version="WKT2_2018", errcheck=False): + """ + This converts a :obj:`~pyproj.crs.CRS` object + to a CF-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. + + Parameters + ---------- + wkt_version: str + Version of WKT supported by ~CRS.to_wkt. + errcheck: bool, optional + If True, will warn when parameters are ignored. Defaults to False. + + Returns + ------- + dict: CF-1.8 version of the projection. + + """ + + cf_dict = {"crs_wkt": self.to_wkt(wkt_version)} + if self.is_geographic and self.name != "unknown": + cf_dict["geographic_crs_name"] = self.name + elif self.is_projected and self.name != "unknown": + cf_dict["projected_crs_name"] = self.name + + proj_dict = self.to_proj4_dict() + proj_name = proj_dict.pop("proj") + if proj_name in ("lonlat", "latlon", "longlat", "latlong"): + grid_mapping_name = "latitude_longitude" + else: + grid_mapping_name = INVERSE_GRID_MAPPING_NAME_MAP.get(proj_name, "unknown") + cf_dict["grid_mapping_name"] = grid_mapping_name + + # get best match for lon_0 value for projetion name + lon_0 = proj_dict.pop("lon_0", None) + if lon_0 is not None: + try: + cf_dict[LON_0_MAP[grid_mapping_name]] = lon_0 + except KeyError: + cf_dict[LON_0_MAP["DEFAULT"]] = lon_0 + + # get best match for k_0 value for projetion name + k_0 = proj_dict.pop("k_0", None) + if k_0 is not None: + try: + cf_dict[K_0_MAP[grid_mapping_name]] = k_0 + except KeyError: + cf_dict[K_0_MAP["DEFAULT"]] = k_0 + + # format the lat_1 and lat_2 for the standard parallel + if "lat_1" in proj_dict and "lat_2" in proj_dict: + cf_dict["standard_parallel"] = [ + proj_dict.pop("lat_1"), + proj_dict.pop("lat_2"), + ] + elif "lat_1" in proj_dict: + cf_dict["standard_parallel"] = proj_dict.pop("lat_1") + + skipped_params = [] + for proj_param, proj_val in proj_dict.items(): + try: + cf_dict[INVERSE_PROJ_PARAM_MAP[proj_param]] = proj_val + except KeyError: + skipped_params.append(proj_param) + + if errcheck and skipped_params: + warnings.warn( + "PROJ parameters not mapped to CF: {}".format(tuple(skipped_params)) + ) + return cf_dict + + @staticmethod + def from_cf(in_cf, errcheck=False): + """ + This converts a CF-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. + + Parameters + ---------- + in_cf: dict + CF version of the projection. + errcheck: bool, optional + If True, will warn when parameters are ignored. Defaults to False. + + Returns + ------- + ~pyproj.crs.CRS: CRS object. + """ + in_cf = in_cf.copy() # preserve user input + if "crs_wkt" in in_cf: + return CRS(in_cf["crs_wkt"]) + elif "spatial_ref" in in_cf: # for previous supported WKT key + return CRS(in_cf["spatial_ref"]) + + grid_mapping_name = in_cf.pop("grid_mapping_name", None) + if grid_mapping_name is None: + raise CRSError("CF projection parameters missing 'grid_mapping_name'") + proj_name = GRID_MAPPING_NAME_MAP.get(grid_mapping_name) + if proj_name is None: + raise CRSError( + "Unsupported grid mapping name: {}".format(grid_mapping_name) + ) + proj_dict = {"proj": proj_name} + if grid_mapping_name == "rotated_latitude_longitude": + proj_dict["o_proj"] = "latlon" + + if "standard_parallel" in in_cf: + standard_parallel = in_cf.pop("standard_parallel") + if isinstance(standard_parallel, list): + proj_dict["lat_1"] = standard_parallel[0] + proj_dict["lat_2"] = standard_parallel[1] + else: + proj_dict["lat_1"] = standard_parallel + + skipped_params = [] + for cf_param, proj_val in in_cf.items(): + try: + proj_dict[PROJ_PARAM_MAP[cf_param]] = proj_val + except KeyError: + skipped_params.append(cf_param) + + if errcheck and skipped_params: + warnings.warn( + "CF parameters not mapped to PROJ: {}".format(tuple(skipped_params)) + ) + + return CRS(proj_dict) + def __reduce__(self): """special method that allows CRS instance to be pickled""" return self.__class__, (self.srs,) diff --git a/test/test_crs_cf.py b/test/test_crs_cf.py new file mode 100644 index 000000000..a7c8de893 --- /dev/null +++ b/test/test_crs_cf.py @@ -0,0 +1,267 @@ +import pytest +from numpy.testing import assert_almost_equal + +from pyproj import CRS +from pyproj.exceptions import CRSError + + +def test_to_cf_transverse_mercator(): + crs = CRS( + init="epsg:3004", towgs84="-122.74,-34.27,-22.83,-1.884,-3.400,-3.030,-15.62" + ) + with pytest.warns(UserWarning): + cf_dict = crs.to_cf(errcheck=True) + towgs84_test = [-122.74, -34.27, -22.83, -1.884, -3.4, -3.03, -15.62] + assert cf_dict.pop("crs_wkt").startswith("BOUNDCRS[") + assert cf_dict == { + "grid_mapping_name": "transverse_mercator", + "latitude_of_projection_origin": 0, + "longitude_of_central_meridian": 15, + "fase_easting": 2520000, + "fase_northing": 0, + "reference_ellipsoid_name": "intl", + "towgs84": towgs84_test, + "unit": "m", + } + assert crs.to_proj4_dict() == { + "proj": "tmerc", + "lat_0": 0, + "lon_0": 15, + "k": 0.9996, + "x_0": 2520000, + "y_0": 0, + "ellps": "intl", + "towgs84": towgs84_test, + "units": "m", + "no_defs": None, + "type": "crs", + } + + +def test_from_cf_transverse_mercator(): + towgs84_test = [-122.74, -34.27, -22.83, -1.884, -3.4, -3.03, -15.62] + crs = CRS.from_cf( + { + "grid_mapping_name": "transverse_mercator", + "latitude_of_projection_origin": 0, + "longitude_of_central_meridian": 15, + "fase_easting": 2520000, + "fase_northing": 0, + "reference_ellipsoid_name": "intl", + "towgs84": towgs84_test, + "unit": "m", + } + ) + with pytest.warns(UserWarning): + cf_dict = crs.to_cf(errcheck=True) + assert len(cf_dict) == 9 + assert cf_dict["crs_wkt"].startswith("BOUNDCRS[") + assert cf_dict["grid_mapping_name"] == "transverse_mercator" + assert cf_dict["latitude_of_projection_origin"] == 0 + assert cf_dict["longitude_of_central_meridian"] == 15 + assert cf_dict["fase_easting"] == 2520000 + assert cf_dict["fase_northing"] == 0 + assert cf_dict["reference_ellipsoid_name"] == "intl" + assert cf_dict["unit"] == "m" + assert_almost_equal(cf_dict["towgs84"], towgs84_test) + + +def test_cf_from_latlon(): + crs = CRS.from_cf( + dict( + grid_mapping_name="latitude_longitude", + semi_major_axis=6378137.0, + inverse_flattening=298.257223, + ) + ) + with pytest.warns(UserWarning): + cf_dict = crs.to_cf(errcheck=True) + assert len(cf_dict) == 4 + assert cf_dict["crs_wkt"].startswith("GEOGCRS[") + assert cf_dict["grid_mapping_name"] == "latitude_longitude" + assert cf_dict["semi_major_axis"] == 6378137.0 + assert cf_dict["inverse_flattening"] == 298.257223 + + +def test_cf_from_latlon__named(): + crs = CRS.from_cf(dict(spatial_ref="epsg:4326")) + with pytest.warns(UserWarning): + cf_dict = crs.to_cf(errcheck=True) + assert cf_dict.pop("crs_wkt").startswith("GEOGCRS[") + assert cf_dict == { + "geographic_crs_name": "WGS 84", + "grid_mapping_name": "latitude_longitude", + "horizontal_datum_name": "WGS84", + } + + +def test_cf_from_utm(): + crs = CRS.from_cf(dict(crs_wkt="epsg:32615")) + with pytest.warns(UserWarning): + cf_dict = crs.to_cf(errcheck=True) + assert cf_dict.pop("crs_wkt").startswith("PROJCRS[") + assert cf_dict == { + "projected_crs_name": "WGS 84 / UTM zone 15N", + "grid_mapping_name": "unknown", + "horizontal_datum_name": "WGS84", + "unit": "m", + } + + +def test_cf_rotated_latlon(): + crs = CRS.from_cf( + dict( + grid_mapping_name="rotated_latitude_longitude", + grid_north_pole_latitude=32.5, + grid_north_pole_longitude=170.0, + ) + ) + assert crs.to_proj4_dict() == { + "proj": "ob_tran", + "o_proj": "latlon", + "o_lat_p": 32.5, + "o_lon_p": 170.0, + "type": "crs", + } + cf_dict = crs.to_cf() + assert cf_dict.pop("crs_wkt").startswith("PROJCRS[") + assert cf_dict == dict( + grid_mapping_name="rotated_latitude_longitude", + grid_north_pole_latitude=32.5, + grid_north_pole_longitude=170.0, + ) + + +def test_cf_rotated_latlon__grid(): + crs = CRS.from_cf( + dict( + grid_mapping_name="rotated_latitude_longitude", + grid_north_pole_latitude=32.5, + grid_north_pole_longitude=170.0, + north_pole_grid_longitude=0, + ) + ) + assert crs.to_proj4_dict() == { + "proj": "ob_tran", + "o_proj": "latlon", + "o_lat_p": 32.5, + "o_lon_p": 170.0, + "lon_0": 0, + "type": "crs", + } + + +def test_cf_lambert_conformal_conic(): + crs = CRS.from_cf( + dict( + grid_mapping_name="lambert_conformal_conic", + standard_parallel=25.0, + longitude_of_central_meridian=265.0, + latitude_of_projection_origin=25.0, + ) + ) + with pytest.warns(UserWarning): + cf_dict = crs.to_cf(errcheck=True) + assert cf_dict.pop("crs_wkt").startswith("PROJCRS[") + assert cf_dict == { + "grid_mapping_name": "lambert_conformal_conic", + "longitude_of_central_meridian": 265, + "scale_factor_at_projection_origin": 1, + "standard_parallel": 25, + "latitude_of_projection_origin": 25, + "fase_easting": 0, + "fase_northing": 0, + "horizontal_datum_name": "WGS84", + "unit": "m", + } + + +def test_cf_lambert_conformal_conic_1sp(): + crs = CRS.from_cf( + dict( + grid_mapping_name="lambert_conformal_conic", + standard_parallel=25.0, + longitude_of_central_meridian=265.0, + latitude_of_projection_origin=25.0, + ) + ) + with pytest.warns(UserWarning): + cf_dict = crs.to_cf(errcheck=True) + assert cf_dict.pop("crs_wkt").startswith("PROJCRS[") + assert cf_dict == { + "grid_mapping_name": "lambert_conformal_conic", + "longitude_of_central_meridian": 265, + "scale_factor_at_projection_origin": 1, + "standard_parallel": 25, + "latitude_of_projection_origin": 25, + "fase_easting": 0, + "fase_northing": 0, + "horizontal_datum_name": "WGS84", + "unit": "m", + } + proj_dict = crs.to_proj4_dict() + assert proj_dict == { + "proj": "lcc", + "lat_1": 25, + "lat_0": 25, + "lon_0": 265, + "k_0": 1, + "x_0": 0, + "y_0": 0, + "datum": "WGS84", + "units": "m", + "no_defs": None, + "type": "crs", + } + + +def test_cf_lambert_conformal_conic_2sp(): + crs = CRS.from_cf( + dict( + grid_mapping_name="lambert_conformal_conic", + standard_parallel=[25.0, 30.0], + longitude_of_central_meridian=265.0, + latitude_of_projection_origin=25.0, + ) + ) + with pytest.warns(UserWarning): + cf_dict = crs.to_cf(errcheck=True) + assert cf_dict.pop("crs_wkt").startswith("PROJCRS[") + assert cf_dict == { + "grid_mapping_name": "lambert_conformal_conic", + "longitude_of_central_meridian": 265, + "standard_parallel": [25, 30], + "latitude_of_projection_origin": 25, + "fase_easting": 0, + "fase_northing": 0, + "horizontal_datum_name": "WGS84", + "unit": "m", + } + proj_dict = crs.to_proj4_dict() + assert proj_dict == { + "proj": "lcc", + "lat_1": 25, + "lat_2": 30, + "lat_0": 25, + "lon_0": 265, + "x_0": 0, + "y_0": 0, + "datum": "WGS84", + "units": "m", + "no_defs": None, + "type": "crs", + } + + +def test_cf_from_invalid(): + with pytest.raises(CRSError): + CRS.from_cf( + dict( + longitude_of_central_meridian=265.0, latitude_of_projection_origin=25.0 + ) + ) + + with pytest.raises(CRSError): + CRS.from_cf( + dict(grid_mapping_name="invalid", latitude_of_projection_origin=25.0) + )