WIP: Add from_proj4 function to create CRS from PROJ.4 definitions

lib/cartopy/_epsg.py

@@ -22,42 +22,20 @@
 from __future__ import (absolute_import, division, print_function)
 
 import cartopy.crs as ccrs
+from cartopy._proj4 import _PROJ4Projection
 import numpy as np
 import pyepsg
-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'}
-
-
-class _EPSGProjection(ccrs.Projection):
+class _EPSGProjection(_PROJ4Projection):
     def __init__(self, code):
         projection = pyepsg.get(code)
         if not isinstance(projection, pyepsg.ProjectedCRS):
             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:
@@ -77,23 +55,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.

lib/cartopy/_proj4.py

@@ -0,0 +1,154 @@
+# (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 <https://www.gnu.org/licenses/>.
+"""
+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'}
+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(x1 - x0, 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, globe=None, bounds=None):
+    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
+
+    proj = proj4_dict['proj']
+    crs_class = PROJ_TO_CRS.get(proj)
+
+    # couldn't find a known CRS class
+    if crs_class is None:
+        # return _PROJ4Projection(proj4_dict, globe=globe, bounds=bounds)
+        # 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)
+    if globe is None:
+        globe = _globe_from_proj4(globe_dict)
+    return crs_class.from_proj4(projection_dict, globe=globe)

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,35 @@ 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):
+        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 +1345,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 +1395,34 @@ 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):
+        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 +2073,8 @@ def epsg(code):
     """
     import cartopy._epsg
     return cartopy._epsg._EPSGProjection(code)
+
+
+def from_proj4(proj4_terms, globe=None, bounds=None):
+    from cartopy._proj4 import from_proj4
+    return from_proj4(proj4_terms, globe, bounds)

lib/cartopy/tests/test_crs.py

@@ -1,4 +1,4 @@
-# (C) British Crown Copyright 2011 - 2017, Met Office
+# (C) British Crown Copyright 2011 - 2018, Met Office
 #
 # This file is part of cartopy.
 #
@@ -225,6 +225,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()