diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py index f658640aa..9ffca4492 100644 --- a/lib/cartopy/crs.py +++ b/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 diff --git a/lib/cartopy/feature.py b/lib/cartopy/feature.py index d91afea26..79757e96c 100644 --- a/lib/cartopy/feature.py +++ b/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.""" diff --git a/lib/cartopy/io/ogc_clients.py b/lib/cartopy/io/ogc_clients.py index e9a6cf9e15..cca64d3af 100644 --- a/lib/cartopy/io/ogc_clients.py +++ b/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 diff --git a/lib/cartopy/mpl/geoaxes.py b/lib/cartopy/mpl/geoaxes.py index 630314ea4..0290cdf02 100644 --- a/lib/cartopy/mpl/geoaxes.py +++ b/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