Proj4 projection

docs/source/contributors.rst

@@ -27,6 +27,7 @@ the package wouldn't be as rich or diverse as it is today:
  * Thomas Lecocq
  * Lion Krischer
  * Nat Wilson
+ * Henry Walshaw
 
 
 Thank you!

lib/cartopy/_epsg.py

@@ -19,22 +19,10 @@
 
 """
 import cartopy.crs as ccrs
-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(ccrs._Proj4Projection):
     def __init__(self, code):
         projection = pyepsg.get(code)
         if not isinstance(projection, pyepsg.ProjectedCRS):
@@ -43,54 +31,11 @@ def __init__(self, code):
         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)
-
-        # Convert lat/lon bounds to projected bounds.
-        # GML defines gmd:EX_GeographicBoundingBox as:
-        #   Geographic area of the entire dataset referenced to WGS 84
-        # NB. We can't use a polygon transform at this stage because
-        # that relies on the existence of the map boundary... the very
-        # thing we're trying to work out! ;-)
         x0, x1, y0, y1 = projection.domain_of_validity()
-        geodetic = ccrs.Geodetic()
-        lons = np.array([x0, x0, x1, x1])
-        lats = np.array([y0, y1, y1, y0])
-        points = self.transform_points(geodetic, lons, lats)
-        x = points[:, 0]
-        y = points[:, 1]
-        self.bounds = (x.min(), x.max(), y.min(), y.max())
-
-    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)])
+        super(_EPSGProjection, self).__init__(proj4_str, x0, x1, y0, y1)
 
-    @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)
+    def __repr__(self):
+        return '_EPSGProjection({})'.format(self.epsg_code)
 
-    @property
-    def threshold(self):
-        x0, x1, y0, y1 = self.bounds
-        return min(x1 - x0, y1 - y0) / 100.

lib/cartopy/crs.py

@@ -1421,6 +1421,73 @@ def _find_gt(a, x):
     return a[0]
 
 
+_GLOBE_PARAMS = {'datum': 'datum',
+                 'ellps': 'ellipse',
+                 'a': 'semimajor_axis',
+                 'b': 'semiminor_axis',
+                 'f': 'flattening',
+                 'rf': 'inverse_flattening',
+                 'towgs84': 'towgs84',
+                 'nadgrids': 'nadgrids'}
+
+
+class _Proj4Projection(Projection):
+    def __init__(self, proj4_str, x0=-180., x1=180., y0=-90., y1=90.):
+        self.proj4_str = proj4_str.strip() #just in case
+
+        terms = [term.strip('+').split('=') for term in self.proj4_str.split(' ')]
+        globe_terms = filter(lambda term: term[0] in _GLOBE_PARAMS, terms)
+        globe = 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(_Proj4Projection, self).__init__(other_terms, globe)
+
+        # Convert lat/lon bounds to projected bounds.
+        # GML defines gmd:EX_GeographicBoundingBox as:
+        #   Geographic area of the entire dataset referenced to WGS 84
+        # NB. We can't use a polygon transform at this stage because
+        # that relies on the existence of the map boundary... the very
+        # thing we're trying to work out! ;-)
+        geodetic = Geodetic()
+        lons = np.array([x0, x0, x1, x1])
+        lats = np.array([y0, y1, y1, y0])
+        points = self.transform_points(geodetic, lons, lats)
+        x = points[:, 0]
+        y = points[:, 1]
+        self.bounds = (x.min(), x.max(), y.min(), y.max())
+
+
+    def __repr__(self):
+        return '_Proj4Projection({})'.format(self.proj4_str)
+
+    @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 epsg(code):
     """
     Return the projection which corresponds to the given EPSG code.
@@ -1436,3 +1503,18 @@ def epsg(code):
     """
     import cartopy._epsg
     return cartopy._epsg._EPSGProjection(code)
+
+
+def proj4(proj4_str, xmin=-180., ymin=-90., xmax=180., ymax=90.):
+    """
+    Return the projection which corresponds to the given proj4 string.
+
+    The proj4 string must correspond to a "projected coordinate system"
+    so EPSG codes such as 4326 (WGS-84) which define a "geodetic coordinate
+    system" will not work.
+
+    .. note::
+        A proj4 string doesn't contain the area of validity, so you can either
+        specify the bounds, or have them default to the entire globe.
+    """
+    return _Proj4Projection(proj4_str, xmin, xmax, ymin, ymax)

lib/cartopy/tests/test_crs.py

@@ -84,6 +84,20 @@ def _check_osgb(self, osgb):
     def test_osgb(self):
         self._check_osgb(ccrs.OSGB())
 
+    def test_proj4(self):
+        #taken from http://epsg.io/27700 for consistency with pyepsg
+        proj4_str = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs"
+        xmin, ymin, xmax, ymax = -8.73, 49.81, 1.83, 60.89
+        uk = ccrs.proj4(proj4_str, xmin, ymin, xmax, ymax)
+        #Assume a bit of rounding error, but this defaults to 7 decimal places
+        #so it's pretty close
+        self.assertAlmostEqual(uk.x_limits, (-83948.465999040171,
+                                             675634.89881823619))
+        self.assertAlmostEqual(uk.y_limits, (-2994.0109472532495,
+                                             1241785.8617898584))
+        self.assertAlmostEqual(uk.threshold, 7595.8336481727638)
+        self._check_osgb(uk)
+
     @unittest.skipIf(pyepsg is None, 'requires pyepsg')
     def test_epsg(self):
         uk = ccrs.epsg(27700)