diff --git a/lib/cartopy/_epsg.py b/lib/cartopy/_epsg.py index 65c03de05..ff3265fd5 100644 --- a/lib/cartopy/_epsg.py +++ b/lib/cartopy/_epsg.py @@ -21,23 +21,13 @@ from __future__ import (absolute_import, division, print_function) +from cartopy._proj4 import _PROJ4Projection import numpy as np -import shapely.geometry as sgeom import cartopy.crs as ccrs -_GLOBE_PARAMS = {'datum': 'datum', - 'ellps': 'ellipse', - 'a': 'semimajor_axis', - 'b': 'semiminor_axis', - 'f': 'flattening', - 'rf': 'inverse_flattening', - 'towgs84': 'towgs84', - 'nadgrids': 'nadgrids'} - - -class _EPSGProjection(ccrs.Projection): +class _EPSGProjection(_PROJ4Projection): def __init__(self, code): import pyepsg projection = pyepsg.get(code) @@ -45,20 +35,8 @@ def __init__(self, code): raise ValueError('EPSG code does not define a projection') self.epsg_code = code - - proj4_str = projection.as_proj4().strip() - terms = [term.strip('+').split('=') for term in proj4_str.split(' ')] - globe_terms = filter(lambda term: term[0] in _GLOBE_PARAMS, terms) - globe = ccrs.Globe(**{_GLOBE_PARAMS[name]: value for name, value in - globe_terms}) - other_terms = [] - for term in terms: - if term[0] not in _GLOBE_PARAMS: - if len(term) == 1: - other_terms.append([term[0], None]) - else: - other_terms.append(term) - super(_EPSGProjection, self).__init__(other_terms, globe) + proj4_str = str(projection.as_proj4().strip()) + super(_EPSGProjection, self).__init__(proj4_str) # Convert lat/lon bounds to projected bounds. # GML defines gmd:EX_GeographicBoundingBox as: @@ -78,23 +56,3 @@ def __init__(self, code): def __repr__(self): return '_EPSGProjection({})'.format(self.epsg_code) - @property - def boundary(self): - x0, x1, y0, y1 = self.bounds - return sgeom.LineString([(x0, y0), (x0, y1), (x1, y1), (x1, y0), - (x0, y0)]) - - @property - def x_limits(self): - x0, x1, y0, y1 = self.bounds - return (x0, x1) - - @property - def y_limits(self): - x0, x1, y0, y1 = self.bounds - return (y0, y1) - - @property - def threshold(self): - x0, x1, y0, y1 = self.bounds - return min(x1 - x0, y1 - y0) / 100. diff --git a/lib/cartopy/_proj4.py b/lib/cartopy/_proj4.py new file mode 100644 index 000000000..1aa6dfd41 --- /dev/null +++ b/lib/cartopy/_proj4.py @@ -0,0 +1,156 @@ +# (C) British Crown Copyright 2014 - 2018, Met Office +# +# This file is part of cartopy. +# +# cartopy is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# cartopy is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with cartopy. If not, see . +""" +Provide support for converting PROJ.4 strings to Projection instances. + +""" + +from __future__ import (absolute_import, division, print_function) + +import cartopy.crs as ccrs +import shapely.geometry as sgeom + +_GLOBE_PARAMS = {'datum': 'datum', + 'ellps': 'ellipse', + 'a': 'semimajor_axis', + 'b': 'semiminor_axis', + 'f': 'flattening', + 'rf': 'inverse_flattening', + 'towgs84': 'towgs84', + 'nadgrids': 'nadgrids'} +# Map PROJ.4 'proj' parameter to CRS class +PROJ_TO_CRS = {} + + +def get_proj4_dict(proj4_terms): + """Convert a PROJ.4 string to a dictionary. + + Parameters + ---------- + proj4_terms: (str, dict, or iterable of key-value pairs) + + Returns + ------- + get_proj4_dict + All PROJ.4 parameters in a dictionary. Any keys with no value + are set to `None`. + + """ + if isinstance(proj4_terms, dict): + return proj4_terms + elif isinstance(proj4_terms, str): + terms = [] + for term in proj4_terms.split(' '): + parts = term.strip('+').split('=') + if len(parts) == 1: + terms.append((parts[0], None)) + else: + terms.append(tuple(parts[:2])) + else: + # assume list of key value pairs + terms = proj4_terms + + return dict(terms) + + +def _split_globe_parameters(proj4_dict): + projection_terms = {} + globe_terms = {} + for name, value in proj4_dict.items(): + if name in _GLOBE_PARAMS: + globe_terms[name] = value + else: + projection_terms[name] = value + return projection_terms, globe_terms + + +def _globe_from_proj4(globe_terms): + """Create a `Globe` object from PROJ.4 parameters.""" + globe = ccrs.Globe(**{_GLOBE_PARAMS[name]: value for name, value in + globe_terms.items()}) + return globe + + +class _PROJ4Projection(ccrs.Projection): + def __init__(self, proj4_terms, globe=None, bounds=None): + terms = get_proj4_dict(proj4_terms) + projection_terms, globe_terms = _split_globe_parameters(terms) + if globe is None: + globe = _globe_from_proj4(globe_terms) + + super(_PROJ4Projection, self).__init__(projection_terms, globe) + + # FIXME: Can we guess at the bounds if not provided? Maybe transform + # an array of points and take the min/max of the result? + # It looks like that's what LambertConformal does. + self.bounds = bounds + + def __repr__(self): + return '_PROJ4Projection({})'.format(self.proj4_init) + + @property + def boundary(self): + x0, x1, y0, y1 = self.bounds + return sgeom.LineString([(x0, y0), (x0, y1), (x1, y1), (x1, y0), + (x0, y0)]) + + @property + def x_limits(self): + x0, x1, y0, y1 = self.bounds + return (x0, x1) + + @property + def y_limits(self): + x0, x1, y0, y1 = self.bounds + return (y0, y1) + + @property + def threshold(self): + x0, x1, y0, y1 = self.bounds + return min(abs(x1 - x0), abs(y1 - y0)) / 100. + + +def _all_subclasses(cls): + return cls.__subclasses__() + [g for s in cls.__subclasses__() + for g in _all_subclasses(s)] + + +def from_proj4(proj4_terms): + proj4_dict = get_proj4_dict(proj4_terms) + + if not PROJ_TO_CRS: + # initialize this here instead of at import + for crs_class in _all_subclasses(ccrs.CRS): + cls_proj = getattr(crs_class, '_proj4_proj', None) + if cls_proj is not None and cls_proj not in PROJ_TO_CRS: + PROJ_TO_CRS[cls_proj] = crs_class + + if 'proj' not in proj4_dict: + raise ValueError("Missing PROJ.4 parameter: proj") + + proj = proj4_dict['proj'] + crs_class = PROJ_TO_CRS.get(proj) + + # couldn't find a known CRS class + if crs_class is None: + # we don't want to allow non-CRS/generic Projection classes + raise ValueError("Projection '{}' is not implemented yet.".format( + proj)) + + projection_dict, globe_dict = _split_globe_parameters(proj4_dict) + globe = _globe_from_proj4(globe_dict) + return crs_class.from_proj4(projection_dict, globe=globe) diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py index 57fb8c061..8a00f737c 100644 --- a/lib/cartopy/crs.py +++ b/lib/cartopy/crs.py @@ -102,6 +102,11 @@ class Projection(six.with_metaclass(ABCMeta, CRS)): 'MultiPolygon': '_project_multipolygon', } + @classmethod + def from_proj4(cls, proj4_dict, **kwargs): + raise NotImplementedError("'{}' can not be created from a " + "PROJ.4 description.".format(cls.__name__)) + @abstractproperty def boundary(self): pass @@ -1039,6 +1044,8 @@ class LambertConformal(Projection): """ + _proj4_proj = 'lcc' + def __init__(self, central_longitude=-96.0, central_latitude=39.0, false_easting=0.0, false_northing=0.0, secant_latitudes=None, standard_parallels=None, @@ -1134,6 +1141,36 @@ def __init__(self, central_longitude=-96.0, central_latitude=39.0, self._x_limits = bounds[0], bounds[2] self._y_limits = bounds[1], bounds[3] + @classmethod + def from_proj4(cls, proj4_dict, **kwargs): + assert proj4_dict.get('proj') == 'lcc' + p_kwargs = {} + + if 'no_defs' in proj4_dict: + lat_1 = proj4_dict.get('lat_1') + lat_2 = proj4_dict.get('lat_2') + else: + lat_1 = proj4_dict.get('lat_1', 33.) + lat_2 = proj4_dict.get('lat_2', 45.) + if lat_1 is not None and lat_2 is not None: + p_kwargs['standard_parallels'] = (float(lat_1), float(lat_2)) + elif lat_1 is not None: + p_kwargs['standard_parallels'] = (float(lat_1),) + elif lat_2 is not None: + raise ValueError("'lat_2' specified without 'lat_1'") + + if 'lon_0' in proj4_dict: + p_kwargs['central_longitude'] = float(proj4_dict['lon_0']) + if 'lat_0' in proj4_dict: + p_kwargs['central_latitude'] = float(proj4_dict['lat_0']) + if 'x_0' in proj4_dict: + p_kwargs['false_easting'] = float(proj4_dict['x_0']) + if 'y_0' in proj4_dict: + p_kwargs['false_northing'] = float(proj4_dict['y_0']) + + kwargs.update(p_kwargs) + return cls(**kwargs) + def __eq__(self, other): res = super(LambertConformal, self).__eq__(other) if hasattr(other, "cutoff"): @@ -1309,6 +1346,9 @@ def y_limits(self): class Stereographic(Projection): + + _proj4_proj = 'stere' + def __init__(self, central_latitude=0.0, central_longitude=0.0, false_easting=0.0, false_northing=0.0, true_scale_latitude=None, @@ -1356,6 +1396,35 @@ def __init__(self, central_latitude=0.0, central_longitude=0.0, self._boundary = sgeom.LinearRing(coords.T) self._threshold = np.diff(self._x_limits)[0] * 1e-3 + @classmethod + def from_proj4(cls, proj4_dict, **kwargs): + assert proj4_dict.get('proj') == 'stere' + p_kwargs = {} + + if 'lon_0' in proj4_dict: + p_kwargs['central_longitude'] = float(proj4_dict['lon_0']) + if 'lat_ts' in proj4_dict: + p_kwargs['true_scale_latitude'] = float(proj4_dict['lat_ts']) + if 'k_0' in proj4_dict: + p_kwargs['scale_factor'] = float(proj4_dict['k_0']) + if 'x_0' in proj4_dict: + p_kwargs['false_easting'] = float(proj4_dict['x_0']) + if 'y_0' in proj4_dict: + p_kwargs['false_northing'] = float(proj4_dict['y_0']) + if 'lat_0' in proj4_dict: + # forced by North and South specific classes + lat_0 = float(proj4_dict['lat_0']) + if lat_0 == -90.: + cls = SouthPolarStereo + elif lat_0 == 90.: + cls = NorthPolarStereo + else: + p_kwargs['central_latitude'] = lat_0 + + kwargs.update(p_kwargs) + return cls(**kwargs) + + @property def boundary(self): return self._boundary @@ -2006,3 +2075,8 @@ def epsg(code): """ import cartopy._epsg return cartopy._epsg._EPSGProjection(code) + + +def from_proj4(proj4_terms): + from cartopy._proj4 import from_proj4 + return from_proj4(proj4_terms) diff --git a/lib/cartopy/tests/test_crs.py b/lib/cartopy/tests/test_crs.py index b54463ce0..431f01d80 100644 --- a/lib/cartopy/tests/test_crs.py +++ b/lib/cartopy/tests/test_crs.py @@ -226,6 +226,46 @@ def test_utm(self): decimal=1) +class TestFromPROJ4(object): + def test_lcc(self): + p_str1 = "+proj=lcc +datum=WGS84 +lat_0=25 +lat_1=25 +lon_0=-95" + crs1 = ccrs.from_proj4(p_str1) + p_str2 = "+proj=lcc +datum=WGS84 +lat_0=25 +lat_1=25 +lon_0=-95 " \ + "+lat_2=45" + crs2 = ccrs.from_proj4(p_str2) + p_str3 = "+proj=lcc +datum=WGS84 +lat_0=25 +lat_1=25 +lon_0=-95 " \ + "+no_defs" + crs3 = ccrs.from_proj4(p_str3) + + assert isinstance(crs1, ccrs.LambertConformal) + assert isinstance(crs2, ccrs.LambertConformal) + assert isinstance(crs3, ccrs.LambertConformal) + x, y = crs1.transform_point(-135.0, 42., ccrs.PlateCarree()) + assert_arr_almost_eq([x, y], + [-3199404.236416136, 2517302.7077927846]) + # crs2 should have the same parameters as crs1 through defaults + x, y = crs2.transform_point(-135.0, 42., ccrs.PlateCarree()) + assert_arr_almost_eq([x, y], + [-3199404.236416136, 2517302.7077927846]) + + assert crs1 == crs2 + + def test_stereographic(self): + p_str1 = "+proj=stere +datum=WGS84 +lat_0=90 +lat_ts=45 +lon_0=-150" + crs1 = ccrs.from_proj4(p_str1) + assert isinstance(crs1, ccrs.NorthPolarStereo) + p_str2 = "+proj=stere +datum=WGS84 +lat_0=-90 +lat_ts=45 +lon_0=-150" + crs2 = ccrs.from_proj4(p_str2) + assert isinstance(crs2, ccrs.SouthPolarStereo) + p_str3 = "+proj=stere +datum=WGS84 +lat_0=80 +lon_0=-150" + crs3 = ccrs.from_proj4(p_str3) + assert isinstance(crs3, ccrs.Stereographic) + + x, y = crs1.transform_point(-145., 74., ccrs.PlateCarree()) + assert_arr_almost_eq([x, y], + [133820.7681726163, -1529578.3794087067]) + + def test_pickle(): # check that we can pickle a simple CRS fh = BytesIO()