diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py index 8a6dd4ccf..c5fc6f3e2 100644 --- a/lib/cartopy/crs.py +++ b/lib/cartopy/crs.py @@ -66,6 +66,12 @@ def _safe_pj_transform(src_crs, tgt_crs, x, y, z=None, trap=True): return transformer.transform(x, y, z, errcheck=trap) +def _ensure_tuple(value): + if np.isscalar(value): + return (value,) + return tuple(value) + + class Globe: """ Define an ellipsoid and, optionally, how to relate it to the real world. @@ -124,6 +130,23 @@ def to_proj4_params(self): ) return OrderedDict((k, v) for k, v in proj4_params if v is not None) + @classmethod + def from_cf(cls, **cf_params): + if cf_params.get('inverse_flattening') == 0: + cf_params['reference_ellipsoid_name'] = 'sphere' + cf_params.pop('inverse_flattening', None) + return cls( + datum=cf_params.get('horizontal_datum_name'), + ellipse=cf_params.get('reference_ellipsoid_name', 'WGS84'), + semimajor_axis=( + cf_params.get('earth_radius') + or cf_params.get('semi_major_axis') + ), + semiminor_axis=cf_params.get('semi_minor_axis'), + inverse_flattening=cf_params.get('inverse_flattening'), + towgs84=cf_params.get('towgs84') + ) + class CRS(_CRS): """ @@ -565,6 +588,13 @@ def transform_vectors(self, src_proj, x, y, u, v): projected_v = vector_magnitudes * np.sin(projected_angles) return projected_u, projected_v + @classmethod + def from_cf(self): + raise NotImplementedError( + "'from_cf' constructor only implemented on specific projections. " + "Try using `cartopy.crs.from_cf` instead." + ) + class Geodetic(CRS): """ @@ -1360,6 +1390,11 @@ def __init__(self, central_longitude=0.0, globe=None): self.threshold = x_max / 360 super().__init__(proj4_params, x_max, y_max, globe=globe) + @classmethod + def from_cf(cls, **cf_params): + globe = Globe.from_cf(**cf_params) + return cls(globe=globe) + def _bbox_and_offset(self, other_plate_carree): """ Return a pair of (xmin, xmax) pairs and an offset which can be used @@ -1500,6 +1535,18 @@ def x_limits(self): def y_limits(self): return (-1e7, 1e7) + @classmethod + def from_cf(cls, **cf_params): + globe = Globe.from_cf(**cf_params) + return cls( + globe=globe, + scale_factor=cf_params['scale_factor_at_central_meridian'], + central_longitude=cf_params['longitude_of_central_meridian'], + central_latitude=cf_params['latitude_of_projection_origin'], + false_easting=cf_params.get('false_easting', 0), + false_northing=cf_params.get('false_northing', 0) + ) + class OSGB(TransverseMercator): def __init__(self, approx=False): @@ -1689,6 +1736,18 @@ def __init__(self, central_longitude=0.0, self.threshold = min(np.diff(self.x_limits)[0] / 720, np.diff(self.y_limits)[0] / 360) + @classmethod + def from_cf(cls, **cf_params): + globe = Globe.from_cf(**cf_params) + return cls( + central_longitude=cf_params['longitude_of_projection_origin'], + latitude_true_scale=cf_params.get('standard_parallel'), + scale_factor=cf_params.get('scale_factor_at_projection_origin'), + false_easting=cf_params.get('false_easting', 0), + false_northing=cf_params.get('false_northing', 0), + globe=globe + ) + def __eq__(self, other): res = super().__eq__(other) if hasattr(other, "_y_limits") and hasattr(other, "_x_limits"): @@ -1738,6 +1797,32 @@ def __init__(self, central_longitude=0.0, globe=None): globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS))] super().__init__(proj4_params, 180, math.degrees(1), globe=globe) + @classmethod + def from_cf(cls, **cf_params): + globe = Globe.from_cf(**cf_params) + if ( + cf_params.get('false_northing', 0) != 0 + or cf_params.get('false_easting', 0) != 0 + ): + raise ValueError( + 'Lambert cylindrical equal area projections with non-zero ' + 'false_northing or false_easting are not supported in cartopy.' + ) + if cf_params.get('standard_parallel', 0) != 0: + raise ValueError( + 'Lambert cylindrical equal area projections with non-zero standard ' + 'parallel are not supported in cartopy.' + ) + if cf_params.get('scale_factor_at_projection_origin', 1) != 1: + raise ValueError( + 'Lambert cylindrical equal area projections with scale ' + 'factor different than 1 are not supported in cartopy.' + ) + return cls( + central_longitude=cf_params['longitude_of_central_meridian'], + globe=globe + ) + class LambertConformal(Projection): """ @@ -1828,6 +1913,18 @@ def __init__(self, central_longitude=-96.0, central_latitude=39.0, self.threshold = 1e5 + @classmethod + def from_cf(cls, **cf_params): + globe = Globe.from_cf(**cf_params) + return cls( + standard_parallels=_ensure_tuple(cf_params['standard_parallel']), + central_longitude=cf_params['longitude_of_central_meridian'], + central_latitude=cf_params['latitude_of_projection_origin'], + false_easting=cf_params.get('false_easting', 0), + false_northing=cf_params.get('false_northing', 0), + globe=globe + ) + def __eq__(self, other): res = super().__eq__(other) if hasattr(other, "cutoff"): @@ -1919,6 +2016,17 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0, self._y_limits = mins[1], maxs[1] self.threshold = np.diff(self._x_limits)[0] * 1e-3 + @classmethod + def from_cf(cls, **cf_params): + globe = Globe.from_cf(**cf_params) + return cls( + central_longitude=cf_params['longitude_of_projection_origin'], + central_latitude=cf_params['latitude_of_projection_origin'], + false_easting=cf_params.get('false_easting', 0), + false_northing=cf_params.get('false_northing', 0), + globe=globe + ) + @property def boundary(self): return self._boundary @@ -1990,6 +2098,16 @@ def __init__(self, pole_longitude=0.0, pole_latitude=90.0, globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS))] super().__init__(proj4_params, 180, 90, globe=globe) + @classmethod + def from_cf(cls, **cf_params): + globe = Globe.from_cf(**cf_params) + return cls( + pole_longitude=cf_params['grid_north_pole_longitude'], + pole_latitude=cf_params['grid_north_pole_latitude'], + central_rotated_longitude=cf_params.get('north_pole_grid_longitude', 0), + globe=globe + ) + class Gnomonic(Projection): _handles_ellipses = False @@ -2074,6 +2192,44 @@ def x_limits(self): def y_limits(self): return self._y_limits + @classmethod + def from_cf(cls, **cf_params): + globe = Globe.from_cf(**cf_params) + central_longitude = cf_params.get( + 'straight_vertical_longitude_from_pole', + cf_params.get('longitude_of_projection_origin') + ) + if central_longitude is None: + raise KeyError('longitude_of_projection_origin') + scale_factor = cf_params.get('scale_factor_at_projection_origin') + false_easting = cf_params.get('false_easting', 0) + false_northing = cf_params.get('false_northing', 0) + central_latitude = cf_params['latitude_of_projection_origin'] + if scale_factor is None and false_northing == 0 and false_easting == 0: + # We can use our PolarStereo classes + if central_latitude == 90: + return NorthPolarStereo( + central_longitude=central_longitude, + true_scale_latitude=cf_params.get('standard_parallel'), + globe=globe + ) + if central_latitude == -90: + return SouthPolarStereo( + central_longitude=central_longitude, + true_scale_latitude=cf_params.get('standard_parallel'), + globe=globe + ) + # Else, we must use a more generic Stereographic + return Stereographic( + central_longitude=central_longitude, + central_latitude=central_latitude, + scale_factor=scale_factor, + false_easting=false_easting, + false_northing=false_northing, + true_scale_latitude=cf_params.get('standard_parallel'), + globe=globe + ) + class NorthPolarStereo(Stereographic): def __init__(self, central_longitude=0.0, true_scale_latitude=None, @@ -2085,6 +2241,7 @@ def __init__(self, central_longitude=0.0, true_scale_latitude=None, globe=globe) + class SouthPolarStereo(Stereographic): def __init__(self, central_longitude=0.0, true_scale_latitude=None, globe=None): @@ -2122,6 +2279,22 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0, self._y_limits = mins[1], maxs[1] self.threshold = np.diff(self._x_limits)[0] * 0.02 + @classmethod + def from_cf(cls, **cf_params): + globe = Globe.from_cf(**cf_params) + if (cf_params.get('false_northing', 0) != 0 + and cf_params.get('false_easting', 0) != 0 + ): + raise ValueError( + 'Orthographic projections with non-zero false northing or ' + 'false easting are not yet supported in cartopy.' + ) + return cls( + central_longitude=cf_params['longitude_of_projection_origin'], + central_latitude=cf_params['latitude_of_projection_origin'], + globe=globe + ) + @property def boundary(self): return self._boundary @@ -2797,6 +2970,27 @@ def __init__(self, central_longitude=0.0, satellite_height=35785831, coords += np.array([[false_easting], [false_northing]]) self._set_boundary(coords) + @classmethod + def from_cf(cls, **cf_params): + globe = Globe.from_cf(**cf_params) + if 'sweep_angle_axis' in cf_params: + sweep_axis = cf_params['sweep_angle_axis'] + else: + sweep_axis = {'x': 'y', 'y': 'x'}[cf_params['fixed_angle_axis']] + if cf_params.get('latitude_of_projection_origin', 0) != 0: + raise ValueError( + "Cartopy doesn't support Geostationary projections with " + "non-zero latitude of projection origin." + ) + return cls( + central_longitude=cf_params['longitude_of_projection_origin'], + satellite_height=cf_params['perspective_point_height'], + false_easting=cf_params.get('false_easting', 0), + false_northing=cf_params.get('false_northing', 0), + sweep_axis=sweep_axis, + globe=globe + ) + class NearsidePerspective(_Satellite): """ @@ -2852,6 +3046,18 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0, false_easting, false_northing, 61) self._set_boundary(coords) + @classmethod + def from_cf(cls, **cf_params): + globe = Globe.from_cf(**cf_params) + return cls( + central_longitude=cf_params['longitude_of_projection_origin'], + central_latitude=cf_params['latitude_of_projection_origin'], + satellite_height=cf_params['perspective_point_height'], + false_easting=cf_params.get('false_easting', 0), + false_northing=cf_params.get('false_northing', 0), + globe=globe + ) + class AlbersEqualArea(Projection): """ @@ -2923,6 +3129,18 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0, self.threshold = 1e5 + @classmethod + def from_cf(cls, **cf_params): + globe = Globe.from_cf(**cf_params) + return cls( + standard_parallels=_ensure_tuple(cf_params['standard_parallel']), + central_longitude=cf_params['longitude_of_central_meridian'], + central_latitude=cf_params['latitude_of_projection_origin'], + false_easting=cf_params.get('false_easting', 0), + false_northing=cf_params.get('false_northing', 0), + globe=globe + ) + @property def boundary(self): return self._boundary @@ -2985,6 +3203,17 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0, self.threshold = 1e5 + @classmethod + def from_cf(cls, **cf_params): + globe = Globe.from_cf(**cf_params) + return cls( + central_longitude=cf_params['longitude_of_projection_origin'], + central_latitude=cf_params['latitude_of_projection_origin'], + false_easting=cf_params.get('false_easting', 0), + false_northing=cf_params.get('false_northing', 0), + globe=globe + ) + @property def boundary(self): return self._boundary @@ -3049,6 +3278,16 @@ def __init__(self, central_longitude=0.0, false_easting=0.0, self._y_limits = mins[1], maxs[1] self.threshold = max(np.abs(self.x_limits + self.y_limits)) * 1e-5 + @classmethod + def from_cf(cls, **cf_params): + globe = Globe.from_cf(**cf_params) + return cls( + central_longitude=cf_params['longitude_of_projection_origin'], + false_easting=cf_params.get('false_easting', 0), + false_northing=cf_params.get('false_northing', 0), + globe=globe + ) + @property def boundary(self): return self._boundary @@ -3213,6 +3452,19 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0, self._y_limits = mercator.y_limits self.threshold = mercator.threshold + @classmethod + def from_cf(cls, **cf_params): + globe = Globe.from_cf(**cf_params) + return cls( + central_longitude=cf_params['longitude_of_projection_origin'], + central_latitude=cf_params['latitude_of_projection_origin'], + scale_factor=cf_params['scale_factor_at_projection_origin'], + azimuth=cf_params['azimuth_of_central_line'], + false_easting=cf_params.get('false_easting', 0), + false_northing=cf_params.get('false_northing', 0), + globe=globe + ) + @property def boundary(self): x0, x1 = self.x_limits @@ -3317,3 +3569,33 @@ def epsg(code): """ import cartopy._epsg return cartopy._epsg._EPSGProjection(code) + + +_cf_mapping = { + 'albers_conical_equal_area': AlbersEqualArea, + 'azimuthal_equidistant': AzimuthalEquidistant, + 'geostationary': Geostationary, + 'lambert_azimuthal_equal_area': LambertAzimuthalEqualArea, + 'lambert_conformal_conic': LambertConformal, + 'lambert_cylindrical_equal_area': LambertCylindrical, + 'latitude_longitude': PlateCarree, + 'mercator': Mercator, + 'oblique_mercator': ObliqueMercator, + 'orthographic': Orthographic, + 'polar_stereographic': Stereographic, + 'rotated_latitude_longitude': RotatedPole, + 'sinusoidal': Sinusoidal, + 'stereographic': Stereographic, + 'transverse_mercator': TransverseMercator, + 'vertical_perspective': NearsidePerspective +} + +def from_cf(**cf_attrs): + """ + Return the projection defined by the grid mapping attributes + as defined by the CF conventions. + """ + name = cf_attrs.pop('grid_mapping_name') + if name not in _cf_mapping: + raise ValueError(f'Grid mapping name "{name}" not known to cartopy.') + return _cf_mapping[name].from_cf(**cf_attrs) diff --git a/lib/cartopy/tests/crs/test_cf.py b/lib/cartopy/tests/crs/test_cf.py new file mode 100644 index 000000000..5b30dbc69 --- /dev/null +++ b/lib/cartopy/tests/crs/test_cf.py @@ -0,0 +1,88 @@ +# Copyright Crown and Cartopy Contributors +# +# This file is part of Cartopy and is released under the BSD 3-clause license. +# See LICENSE in the root of the repository for full licensing details. +""" +Tests for creating CRS from CF attributes +""" +import numpy as np +import pytest + +import cartopy.crs as ccrs + + +# All taken from examples of the CF-Convenvtions, section 5.6 +CRS = [ + ( + {'grid_mapping_name': 'latitude_longitude'}, + ccrs.PlateCarree() + ), + ( + { + 'grid_mapping_name': 'latitude_longitude', + 'semi_major_axis': 6371000.0, + 'inverse_flattening': 0, + }, + ccrs.PlateCarree(globe=ccrs.Globe( + semimajor_axis=6371000.0, + ellipse='sphere' + )) + ), + ( + { + 'grid_mapping_name': 'rotated_latitude_longitude', + 'grid_north_pole_latitude': 36.0, + 'grid_north_pole_longitude': 74.0 + }, + ccrs.RotatedPole(pole_latitude=36, pole_longitude=74) + ), + ( + { + 'grid_mapping_name': 'lambert_conformal_conic', + 'standard_parallel': 25.0, + 'longitude_of_central_meridian': 265.0, + 'latitude_of_projection_origin': 25.0 + }, + ccrs.LambertConformal( + central_longitude=265.0, + central_latitude=25.0, + standard_parallels=(25,) + ) + ), +] + + +@pytest.mark.parametrize("cfattrs,exp", CRS) +def test_from_cf(cfattrs, exp): + crs = ccrs.from_cf(**cfattrs) + assert crs == exp + + +def test_from_cf_OSGB(): + # These attrs are from the cf-conventions examples + # We can't compare the TransverseMercator created with those directly with OSGB + # So we compare values from a transform + crs = ccrs.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, + "unit": "metre", + } + ) + exp = ccrs.OSGB() + + # Arbitrary points, all within the british isles + ptsLon = np.array([-4.256495470173718, -0.23201298137905724, -4.5156834308392035]) + ptsLat = np.array([58.3199544202207960, 52.59337104346588, 50.3827602907713]) + + cfpts = crs.transform_points(ccrs.PlateCarree(), ptsLon, ptsLat) + exppts = exp.transform_points(ccrs.PlateCarree(), ptsLon, ptsLat) + + np.testing.assert_array_equal(cfpts, exppts)