Change stereographic projection

lib/cartopy/crs.py

@@ -1098,12 +1098,14 @@ def y_limits(self):
 class Stereographic(Projection):
     def __init__(self, central_latitude=0.0, central_longitude=0.0,
                  false_easting=0.0, false_northing=0.0,
-                 true_scale_latitude=None, globe=None):
+                 true_scale_latitude=None, scaling_factor=None, globe=None):
         proj4_params = [('proj', 'stere'), ('lat_0', central_latitude),
                         ('lon_0', central_longitude),
                         ('x_0', false_easting), ('y_0', false_northing)]
         if true_scale_latitude:
             proj4_params.append(('lat_ts', true_scale_latitude))
+        if scaling_factor:
+            proj4_params.append(('k_0', scaling_factor))
         super(Stereographic, self).__init__(proj4_params, globe=globe)
 
         # TODO: Factor this out, particularly if there are other places using

lib/cartopy/feature.py

@@ -29,6 +29,7 @@
 import shapely.geometry
 import six
 
+from cartopy.io.ogc_clients import WFSGeometrySource
 import cartopy.io.shapereader as shapereader
 import cartopy.crs
 
@@ -292,6 +293,48 @@ def intersecting_geometries(self, extent):
                     yield geom
 
 
+class WFSFeature(Feature):
+    """
+    A class capable of drawing a collection of geometries
+    obtained from a OGC Web Feature Service (WFS).
+
+    """
+    def __init__(self, wfs, features, **kwargs):
+        """
+        Args:
+
+        * wfs: string or :class:`owslib.wfs.WebFeatureService` instance
+            The WebFeatureService instance, or URL of a WFS service, from which
+            to retrieve the geometries.
+
+        * features: string or list of strings
+            The typename(s) of features available from the web service that
+            will be retrieved. Somewhat analogous to layers in WMS/WMTS.
+
+        Kwargs:
+            Keyword arguments to be used when drawing this feature.
+
+        """
+
+        self.source = WFSGeometrySource(wfs, features)
+        crs = self.source.default_projection()
+        super(WFSFeature, self).__init__(crs, **kwargs)
+        # Default kwargs
+        self._kwargs.setdefault('edgecolor', 'red')
+        self._kwargs.setdefault('facecolor', 'none')
+
+    def geometries(self):
+        min_x, min_y, max_x, max_y = self.crs.boundary.bounds
+        geoms = self.source.fetch_geometries(self.crs,
+                                             extent=(min_x, max_x,
+                                                     min_y, max_y))
+        return iter(geoms)
+
+    def intersecting_geometries(self, extent):
+        geoms = self.source.fetch_geometries(self.crs, extent)
+        return iter(geoms)
+
+
 BORDERS = NaturalEarthFeature('cultural', 'admin_0_boundary_lines_land',
                               '110m', edgecolor='black', facecolor='none')
 """Small scale (1:110m) country boundaries."""

lib/cartopy/io/ogc_clients.py

@@ -32,24 +32,30 @@
 
 import io
 import math
+import os
+import warnings
 import weakref
+from xml.etree import ElementTree
 
 from PIL import Image
+from shapely.geometry import LinearRing, LineString, Point
 
 try:
     from owslib.wms import WebMapService
+    from owslib.wfs import WebFeatureService
     import owslib.util
     import owslib.wmts
     _OWSLIB_AVAILABLE = True
 except ImportError:
     WebMapService = None
+    WebFeatureService = None
     _OWSLIB_AVAILABLE = False
 
 from cartopy.io import RasterSource
 import cartopy.crs as ccrs
 
 
-_OWSLIB_REQUIRED = 'OWSLib is required to use the WMS or WMTS source.'
+_OWSLIB_REQUIRED = 'OWSLib is required to use OGC web services.'
 
 # Hardcode some known EPSG codes for now.
 _CRS_TO_OGC_SRS = {ccrs.PlateCarree(): 'EPSG:4326'
@@ -66,8 +72,16 @@
 }
 
 _URN_TO_CRS = {
+    'urn:ogc:def:crs:EPSG::4326': ccrs.PlateCarree(),
     'urn:ogc:def:crs:EPSG::900913': ccrs.GOOGLE_MERCATOR,
+    'urn:ogc:def:crs:EPSG::32661': ccrs.Stereographic(central_latitude=90,
+                                                      false_easting=2000000,
+                                                      false_northing=2000000,
+                                                      true_scale_latitude=90,
+                                                      scaling_factor=0.994),
     'urn:ogc:def:crs:OGC:1.3:CRS84': ccrs.PlateCarree(),
+    'urn:ogc:def:crs:EPSG::3031': ccrs.Stereographic(central_latitude=-90,
+                                                     true_scale_latitude=-71)
 }
 
 
@@ -407,3 +421,238 @@ def _wmts_images(self, wmts, layer, matrix_set_name, extent,
             img_extent = (min_img_x, min_img_x + n_cols * tile_span_x,
                           max_img_y - n_rows * tile_span_y, max_img_y)
         return big_img, img_extent
+
+
+class WFSGeometrySource(object):
+    """Web Feature Service (WFS) retrieval for Cartopy."""
+
+    def __init__(self, service, features, getfeature_extra_kwargs=None):
+        """
+        Args:
+
+        * service:
+            The URL of a WFS, or an instance of
+            :class:`owslib.wfs.WebFeatureService`.
+        * features:
+            The typename(s) of the features from the WFS that
+            will be retrieved and made available as geometries.
+
+        Kwargs:
+
+        * getfeature_extra_kwargs:
+            Extra keyword args to pass to the service's `getfeature` call.
+
+        """
+        if WebFeatureService is None:
+            raise ImportError(_OWSLIB_REQUIRED)
+
+        if isinstance(service, basestring):
+            service = WebFeatureService(service)
+
+        if isinstance(features, basestring):
+            features = [features]
+
+        # Populate an empty kwargs dict with something harmless.
+        if getfeature_extra_kwargs is None:
+            #getfeature_extra_kwargs = {'propertyname': ['*']}
+            getfeature_extra_kwargs = {}
+
+        if not features:
+            raise ValueError('One or more features must be specified.')
+        for feature in features:
+            if feature not in service.contents:
+                raise ValueError('The {!r} feature does not exist in this '
+                                 'service.'.format(feature))
+
+        self.service = service
+        self.features = features
+        self.getfeature_extra_kwargs = getfeature_extra_kwargs
+
+        self._default_urn = None
+        # For parsing GML. The enclosing {} are required as part of this.
+        self.ms = "{http://mapserver.gis.umn.edu/mapserver}"
+        self.gml = "{http://www.opengis.net/gml}"
+
+    def default_projection(self):
+        """
+        Return a :class:`cartopy.crs.Projection` in which the WFS
+        service can supply the requested features.
+
+        """
+        # Using first element in crsOptions (default).
+        if self._default_urn is None:
+            default_urn = set(self.service.contents[feature].crsOptions[0] for
+                              feature in self.features)
+            if len(default_urn) != 1:
+                ValueError('Failed to find a single common default SRS '
+                           'across all features (typenames).')
+            else:
+                default_urn = default_urn.pop()
+                default_srs = default_urn.id
+
+            if unicode(default_urn) not in _URN_TO_CRS:
+                raise ValueError('Unknown mapping from SRS/CRS_URN {!r} to '
+                                 'cartopy projection.'.format(default_srs))
+
+            self._default_urn = default_urn
+
+        return _URN_TO_CRS[unicode(self._default_urn)]
+
+    def fetch_geometries(self, projection, extent):
+        """
+        Return any Point, Linestring or LinearRing geometries available
+        from the WFS that lie within the specified extent.
+
+        Args:
+
+        * projection: :class:`cartopy.crs.Projection`
+            The projection in which the extent is specified and in
+            which the geometries should be returned. Only the default
+            (native) projection is supported.
+
+        * extent: four element tuple
+            (min_x, max_x, min_y, max_y) tuple defining the geographic extent
+            of the geometries to obtain.
+
+        Returns:
+            A list of Shapely geometries.
+
+        """
+        if self.default_projection() != projection:
+            raise ValueError('Geometries are only available in projection '
+                             '{!r}.'.format(self.default_projection()))
+
+        min_x, max_x, min_y, max_y = extent
+        response = self.service.getfeature(typename=self.features,
+                                           bbox=(min_x, min_y, max_x, max_y),
+                                           **self.getfeature_extra_kwargs)
+        geoms_by_srs = self._to_shapely_geoms(response)
+        if not geoms_by_srs:
+            geoms = []
+        elif len(geoms_by_srs) > 1:
+            raise ValueError('Unexpected response from the WFS server. The '
+                             'geometries are in multiple SRSs, when only one '
+                             'was expected.')
+        else:
+            srs, geoms = geoms_by_srs.items()[0]
+            # Attempt to verify the SRS associated with the geometries (if any)
+            # matches the specified projection.
+            if srs is not None:
+                if srs in _URN_TO_CRS:
+                    geom_proj = _URN_TO_CRS[srs]
+                    if geom_proj != projection:
+                        raise ValueError('The geometries are not in expected '
+                                         'projection. Expected {!r}, got '
+                                         '{!r}.'.format(projection, geom_proj))
+                else:
+                    warnings.warn('Unable to verify matching projections due '
+                                  'to incomplete mappings from srs identifiers '
+                                  'to cartopy projections. The geometries have '
+                                  'an SRS of {!r}.'.format(srs))
+
+        return geoms
+
+    def _to_shapely_geoms(self, response):
+        """
+        Convert polygon coordinate strings in WFS response XML to Shapely
+        geometries.
+
+        Args:
+
+        * response: (file-like object)
+            WFS response XML data.
+
+        Returns:
+            A dictionary containing geometries, with key-value pairs of
+            the form {srsname: [geoms]}.
+
+        """
+        linear_rings_data = []
+        linestrings_data = []
+        points_data = []
+        tree = ElementTree.parse(response)
+
+        for node in tree.findall('.//{ms}msGeometry'.format(ms=self.ms)):
+            # Find LinearRing geometries in our Polygon node.
+            find_str = './/{gml}LinearRing'.format(gml=self.gml)
+            if self._node_has_child(node, find_str):
+                data = self._find_polygon_coords(node, find_str)
+                linear_rings_data.extend(data)
+
+            # Find LineString geometries in our Polygon node.
+            find_str = './/{gml}LineString'.format(gml=self.gml)
+            if self._node_has_child(node, find_str):
+                data = self._find_polygon_coords(node, find_str)
+                linestrings_data.extend(data)
+
+            # Find Point geometries in our Polygon node.
+            find_str = './/{gml}Point'.format(gml=self.gml)
+            if self._node_has_child(node, find_str):
+                data = self._find_polygon_coords(node, find_str)
+                points_data.extend(data)
+
+        geoms_by_srs = {}
+        for srs, x, y in linear_rings_data:
+            geoms_by_srs.setdefault(srs, []).append(LinearRing(zip(x, y)))
+        for srs, x, y in linestrings_data:
+            geoms_by_srs.setdefault(srs, []).append(LineString(zip(x, y)))
+        for srs, x, y in points_data:
+            geoms_by_srs.setdefault(srs, []).append(Point(zip(x, y)))
+        return geoms_by_srs
+
+    def _find_polygon_coords(self, node, find_str):
+        """
+        Return the x, y coordinate values for all the geometries in
+        a given`node`.
+
+        Args:
+
+        * node: :class:`xml.etree.ElementTree.Element`
+            Node of the parsed XML response. 
+
+        * find_str: string
+            A search string used to match subelements that contain
+            the coordinates of interest, for example:
+            './/{http://www.opengis.net/gml}LineString'
+
+        Returns:
+            A list of (srsName, x_vals, y_vals) tuples.
+
+        """
+
+        data = []
+        for polygon in node.findall(find_str):
+            feature_srs = polygon.attrib.get('srsName')
+            x, y = [], []
+
+            # We can have nodes called `coordinates` or `coord`.
+            coordinates_find_str = '{}coordinates'.format(self.gml)
+            coords_find_str = '{}coord'.format(self.gml)
+
+            if self._node_has_child(polygon, coordinates_find_str): 
+                points = polygon.findtext(coordinates_find_str)
+                coords = points.strip().split(' ')
+                for coord in coords:
+                    x_val, y_val = coord.split(',')
+                    x.append(float(x_val))
+                    y.append(float(y_val))
+            elif self._node_has_child(polygon, coords_find_str):
+                for coord in polygon.findall(coords_find_str):
+                    x.append(float(coord.findtext('{}X'.format(self.gml))))
+                    y.append(float(coord.findtext('{}Y'.format(self.gml))))
+            else:
+                raise ValueError('Unable to find or parse coordinate values '
+                                 'from the XML.')
+
+            data.append((feature_srs, x, y))
+        return data
+
+    @staticmethod
+    def _node_has_child(node, find_str):
+        """
+        Return whether `node` contains (at any sub-level), a node with name
+        equal to `find_str`.
+
+        """
+        element = node.find(find_str)
+        return element is not None

lib/cartopy/mpl/geoaxes.py

@@ -741,6 +741,33 @@ def add_raster(self, raster_source, **slippy_image_kwargs):
             self.add_image(img)
         return img
 
+    def add_web_feature(self, geoms_source, crs, **kwargs):
+        """
+        Add the given shapely geometries (in the given crs) from a
+        Web Feature Service (WFS) to the axes.
+
+        """
+        ax = self.get_axes()
+        # Fail early if the WFS source cannot provide features in the
+        # current projection.
+        #geoms_source.validate_projection(self.projection)
+
+        [x1, y1], [x2, y2] = ax.viewLim.get_points()
+        features, extent = geoms_source.fetch_raster(ax.projection,
+                                                     extent=[x1, x2, y1, y2])
+        if features is None or extent is None:
+            raise ValueError('Uh-oh.')
+        with ax.hold_limits():
+            self.set_extent(extent)
+
+        geoms = []
+        for source in features:
+            if source.values():
+                for polygon in source.values():
+                    geoms.extend(polygon)
+        feature = cartopy.feature.ShapelyFeature(geoms, crs, **kwargs)
+        return self.add_feature(feature)
+
     def _regrid_shape_aspect(self, regrid_shape, target_extent):
         """
         Helper for setting regridding shape which is used in several