From d22e0d48d8500ee8e8bfcfbd21c455fa74647d6d Mon Sep 17 00:00:00 2001 From: Peter Killick Date: Fri, 31 Oct 2014 11:55:46 +0000 Subject: [PATCH 1/4] Nearly-working bones of adding WFS to Cartopy --- lib/cartopy/io/ogc_clients.py | 202 +++++++++++++++++++++++++++++++++- lib/cartopy/mpl/geoaxes.py | 27 +++++ 2 files changed, 228 insertions(+), 1 deletion(-) diff --git a/lib/cartopy/io/ogc_clients.py b/lib/cartopy/io/ogc_clients.py index e9a6cf9e15..d683c0355 100644 --- a/lib/cartopy/io/ogc_clients.py +++ b/lib/cartopy/io/ogc_clients.py @@ -32,24 +32,29 @@ import io import math +import os 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,7 +71,9 @@ } _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.Mercator(), 'urn:ogc:def:crs:OGC:1.3:CRS84': ccrs.PlateCarree(), } @@ -407,3 +414,196 @@ 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 WFSRasterSource(RasterSource): + """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 :class:`owslib.wfs.WebFeatureService` + instance. + * features: + The name(s) of the features from the WFS to use. + + 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': ['*']} + + if len(features) == 0: + 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 + + # 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}" + + self._crs_urn_for_projection_id = {} + + def _crs_urn(self, projection): + """ + Confirm that the SRS of the features in the WFS response is recognised + and that it lines up with the requested `projection`. + + """ + key = id(projection) + crs_urn = self._crs_urn_for_projection_id.get(key) + if crs_urn is None: + # Confirm the CRS URNs are recognised and match `projection`. + for feature in self.features: + crs_options = self.service.contents[feature].crsOptions + for features_crs_urn in crs_options: + if features_crs_urn in _URN_TO_CRS: + this_crs = _URN_TO_CRS[features_crs_urn] + if this_crs == projection: + crs_urn = features_crs_urn + break + # Fail informatively if not. + if crs_urn is None: + raise ValueError('The projection {!r} was not convertible to ' + 'a suitable WFS SRS.'.format(projection)) + self._crs_urn_for_projection_id[key] = crs_urn + return crs_urn + + def validate_projection(self, projection): + self._crs_urn(projection) + + def fetch_raster(self, projection, extent): + """ + Return Shapely geometries of any Point, Linestring or LinearRing + geometries found in the WFS request. + + .. note:: + Kwargs passed to the get feature request are not guaranteed to + have an impact on the response data as the WFS server will not + necessarily be able to supply data that precisely matches the + request. In such cases the response will contain the default + data provided by the WFS server. + + """ + service = self.service + min_x, max_x, min_y, max_y = extent + response = service.getfeature(typename=self.features, + bbox=(min_x, min_y, max_x, max_y), + **self.getfeature_extra_kwargs) + wfs_features = self._to_shapely_geoms(response) + return wfs_features, extent + + def _to_shapely_geoms(self, response): + """ + Convert polygon coordinate strings in WFS response XML to Shapely + geometries. + + Args: + + * response: + WFS response XML data. + + """ + 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 _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 _node_has_child(node, find_str): + data = self._find_polygon_coords(node, find_str) + points_data.extend(data) + + linear_rings = {} + for k, x, y in linear_rings_data: + linear_rings.setdefault(k, []).append(LinearRing(zip(x, y))) + linestrings = {} + for k, x, y in linestrings_data: + linestrings.setdefault(k, []).append(LineString(zip(x, y))) + points = {} + for k, x, y in points_data: + points.setdefault(k, []).append(Point(zip(x, y))) + return linear_rings, linestrings, points + + def _find_polygon_coords(self, node, find_str): + """Find all coordinates data for a given Polygon `node`.""" + + data = [] + for polygon in node.findall(find_str): + feature_srs = polygon.attrib.get('srsName') + # Assume we can either have nodes called `coordinates` or `coords`. + find_str = '{}coordinates'.format(self.gml) + if _node_has_child(polygon, find_str): + points = polygon.findtext(find_str) + x, y = self._coordinates_to_points(points) + else: + x, y = [], [] + for coord in polygon.findall('{}coord'.format(self.gml)): + x.append(float(coord.findtext('{}X'.format(self.gml)))) + y.append(float(coord.findtext('{}Y'.format(self.gml)))) + + data.append([feature_srs, x, y]) + return data + + def _node_has_child(self, node, find_str): + """ + Determine whether `node` contains (at any sub-level), a node with name + equal to `find_str`. + + """ + found = list(node.iterfind(find_str)) + return len(found) > 0 + + def _coordinates_to_points(self, points): + """ + If our XML has a coordinates node then we'll get one string of + whitespace-separated points: + + "x,y x,y ... x,y" + + Parses this string to return lists of x and y points. + + """ + x = [] + y = [] + coords = points.strip().split(' ') + for coord in coords: + x_val, y_val = coord.split(',') + x.append(float(x_val)) + y.append(float(y_val)) + return x, y 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 From d664e67a4e6d0f076be169dc2cc65cfc54312a69 Mon Sep 17 00:00:00 2001 From: Ed Campbell Date: Mon, 3 Nov 2014 17:24:54 +0000 Subject: [PATCH 2/4] Initial pass --- lib/cartopy/feature.py | 27 ++++++ lib/cartopy/io/ogc_clients.py | 171 ++++++++++++++++++++++++---------- 2 files changed, 149 insertions(+), 49 deletions(-) diff --git a/lib/cartopy/feature.py b/lib/cartopy/feature.py index d91afea26..bdaf79f7a 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,32 @@ def intersecting_geometries(self, extent): yield geom +class WFSFeature(Feature): + """ + A class capable of drawing a collection of geometries + obtained from a Web Ferature Service (WFS) + + """ + def __init__(self, wfs, features, **kwargs): + 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 d683c0355..11f133377 100644 --- a/lib/cartopy/io/ogc_clients.py +++ b/lib/cartopy/io/ogc_clients.py @@ -75,6 +75,8 @@ 'urn:ogc:def:crs:EPSG::900913': ccrs.GOOGLE_MERCATOR, 'urn:ogc:def:crs:EPSG::32661': ccrs.Mercator(), '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) } @@ -416,7 +418,7 @@ def _wmts_images(self, wmts, layer, matrix_set_name, extent, return big_img, img_extent -class WFSRasterSource(RasterSource): +class WFSGeometrySource(RasterSource): """Web Feature Service (WFS) retrieval for Cartopy.""" def __init__(self, service, features, getfeature_extra_kwargs=None): @@ -427,7 +429,7 @@ def __init__(self, service, features, getfeature_extra_kwargs=None): The URL of a WFS, or an :class:`owslib.wfs.WebFeatureService` instance. * features: - The name(s) of the features from the WFS to use. + The typename(s) of the features from the WFS to use. Kwargs: @@ -448,7 +450,7 @@ def __init__(self, service, features, getfeature_extra_kwargs=None): if getfeature_extra_kwargs is None: getfeature_extra_kwargs = {'propertyname': ['*']} - if len(features) == 0: + if not features: raise ValueError('One or more features must be specified.') for feature in features: if feature not in service.contents: @@ -463,56 +465,126 @@ def __init__(self, service, features, getfeature_extra_kwargs=None): self.ms = "{http://mapserver.gis.umn.edu/mapserver}" self.gml = "{http://www.opengis.net/gml}" - self._crs_urn_for_projection_id = {} + #self._crs_urn_for_projection_id = {} + self._srs_for_projection_id = {} - def _crs_urn(self, projection): + #def _crs_urn(self, projection): + # """ + # Confirm that the SRS of the features in the WFS response is recognised + # and that it lines up with the requested `projection`. +# +# """ +# key = id(projection) +# crs_urn = self._crs_urn_for_projection_id.get(key) +# if crs_urn is None: +# # Confirm the CRS URNs are recognised and match `projection`. +# for feature in self.features: +# crs_options = self.service.contents[feature].crsOptions +# for features_crs_urn in crs_options: +# if features_crs_urn in _URN_TO_CRS: +# this_crs = _URN_TO_CRS[features_crs_urn] +# if this_crs == projection: +# crs_urn = features_crs_urn +# break +# # Fail informatively if not. +# if crs_urn is None: +# raise ValueError('The projection {!r} was not convertible to ' +# 'a suitable WFS SRS.'.format(projection)) +# self._crs_urn_for_projection_id[key] = crs_urn +# return crs_urn + + #def validate_projection(self, projection): + # self._crs_urn(projection) + + #def _srs(self, projection): + # key = id(projection) + # srs = self._srs_for_projection_id.get(key) + # if srs is None: + # srs = _CRS_TO_OGC_SRS.get(projection) + # if srs is None: + # raise ValueError('The projection {!r} was not convertible to ' + # 'a suitable OGC SRS.'.format(projection)) + # for feature in self.features: + # if srs not in self.service.contents[feature].crsOptions: + # raise ValueError('The SRS {} is not a valid SRS for the ' + # '{!r} WFS feature.'.format(srs, feature)) + # self._srs_for_projection_id[key] = srs + # return srs + + def default_projection(self): """ - Confirm that the SRS of the features in the WFS response is recognised - and that it lines up with the requested `projection`. + Return a :class:`cartopy.crs.Projection` and corresponding + in which the WFS service can supply the requested features. """ - key = id(projection) - crs_urn = self._crs_urn_for_projection_id.get(key) - if crs_urn is None: - # Confirm the CRS URNs are recognised and match `projection`. - for feature in self.features: - crs_options = self.service.contents[feature].crsOptions - for features_crs_urn in crs_options: - if features_crs_urn in _URN_TO_CRS: - this_crs = _URN_TO_CRS[features_crs_urn] - if this_crs == projection: - crs_urn = features_crs_urn - break - # Fail informatively if not. - if crs_urn is None: - raise ValueError('The projection {!r} was not convertible to ' - 'a suitable WFS SRS.'.format(projection)) - self._crs_urn_for_projection_id[key] = crs_urn - return crs_urn + # Using first element in crsOptions (default) + #import pdb;pdb.set_trace() + default_srs = set(self.service.contents[feature].crsOptions[0] for + feature in self.features) + if len(default_srs) != 1: + ValueError('Failed to find a single common default SRS ' + 'across all features (typenames).') + else: + default_srs = unicode(list(default_srs)[0]) + + if default_srs not in _URN_TO_CRS: + raise ValueError('Unknown mapping from SRS/CRS_URN {!r} to ' + 'cartopy projection.'.format(default_srs)) + + return _URN_TO_CRS[default_srs] + + #for proj, srs in _CRS_TO_OGC_SRS.items(): + # missing = False + # for feature in self.features: + # if srs not in self.service.contents[feature].crsOptions: + # missing = True + # if not missing: + # break + #if missing: + # raise ValueError('The requested features are not available in ' + # 'any of the supported SRSs.') + #return proj - def validate_projection(self, projection): - self._crs_urn(projection) - def fetch_raster(self, projection, extent): + def fetch_geometries(self, projection, extent): """ Return Shapely geometries of any Point, Linestring or LinearRing geometries found in the WFS request. .. note:: + Kwargs passed to the get feature request are not guaranteed to have an impact on the response data as the WFS server will not necessarily be able to supply data that precisely matches the request. In such cases the response will contain the default data provided by the WFS server. + Returns: + ... + """ - service = self.service - min_x, max_x, min_y, max_y = extent - response = service.getfeature(typename=self.features, - bbox=(min_x, min_y, max_x, max_y), - **self.getfeature_extra_kwargs) - wfs_features = self._to_shapely_geoms(response) - return wfs_features, extent + #import pdb;pdb.set_trace() + min_x, max_x, min_y, max_y = extent #TODO - handle non-native proj. ie. convert extent to query wfs - may depend on + # on default srs. + + response = self.service.getfeature(typename=self.features, + #srsname=self._srs(projection), - doesn't nec. work - could try query params + bbox=(min_x, min_y, max_x, max_y), # - what proj????? + **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('Geoms in multiple SRSs') # DEBUG code + else: + srs, geoms = geoms_by_srs.items()[0] + if srs is not None: + # Verify SRS. + geom_proj = _URN_TO_CRS[srs] + if geom_proj != projection: + raise ValueError('Geometries are not in requested proj - need to project...') + + return geoms def _to_shapely_geoms(self, response): """ @@ -528,6 +600,7 @@ def _to_shapely_geoms(self, response): linear_rings_data = [] linestrings_data = [] points_data = [] + #import pdb;pdb.set_trace() tree = ElementTree.parse(response) for node in tree.findall('.//{ms}msGeometry'.format(ms=self.ms)): @@ -539,26 +612,24 @@ def _to_shapely_geoms(self, response): # Find LineString geometries in our Polygon node. find_str = './/{gml}LineString'.format(gml=self.gml) - if _node_has_child(node, find_str): + 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 _node_has_child(node, find_str): + if self._node_has_child(node, find_str): data = self._find_polygon_coords(node, find_str) points_data.extend(data) - linear_rings = {} - for k, x, y in linear_rings_data: - linear_rings.setdefault(k, []).append(LinearRing(zip(x, y))) - linestrings = {} - for k, x, y in linestrings_data: - linestrings.setdefault(k, []).append(LineString(zip(x, y))) - points = {} - for k, x, y in points_data: - points.setdefault(k, []).append(Point(zip(x, y))) - return linear_rings, linestrings, points + 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): """Find all coordinates data for a given Polygon `node`.""" @@ -568,7 +639,7 @@ def _find_polygon_coords(self, node, find_str): feature_srs = polygon.attrib.get('srsName') # Assume we can either have nodes called `coordinates` or `coords`. find_str = '{}coordinates'.format(self.gml) - if _node_has_child(polygon, find_str): + if self._node_has_child(polygon, find_str): points = polygon.findtext(find_str) x, y = self._coordinates_to_points(points) else: @@ -586,8 +657,9 @@ def _node_has_child(self, node, find_str): equal to `find_str`. """ + #TODO Can this be done with a generator?? found = list(node.iterfind(find_str)) - return len(found) > 0 + return bool(found) def _coordinates_to_points(self, points): """ @@ -599,6 +671,7 @@ def _coordinates_to_points(self, points): Parses this string to return lists of x and y points. """ + #TODO optimise this using itertools. x = [] y = [] coords = points.strip().split(' ') From a9789c04fa13137aaac49c0a6d427269b324be4e Mon Sep 17 00:00:00 2001 From: Ed Campbell Date: Thu, 6 Nov 2014 16:31:30 +0000 Subject: [PATCH 3/4] further progress --- lib/cartopy/feature.py | 18 ++- lib/cartopy/io/ogc_clients.py | 242 +++++++++++++++------------------- 2 files changed, 124 insertions(+), 136 deletions(-) diff --git a/lib/cartopy/feature.py b/lib/cartopy/feature.py index bdaf79f7a..79757e96c 100644 --- a/lib/cartopy/feature.py +++ b/lib/cartopy/feature.py @@ -296,10 +296,26 @@ def intersecting_geometries(self, extent): class WFSFeature(Feature): """ A class capable of drawing a collection of geometries - obtained from a Web Ferature Service (WFS) + 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) diff --git a/lib/cartopy/io/ogc_clients.py b/lib/cartopy/io/ogc_clients.py index 11f133377..f53b1c297 100644 --- a/lib/cartopy/io/ogc_clients.py +++ b/lib/cartopy/io/ogc_clients.py @@ -33,6 +33,7 @@ import io import math import os +import warnings import weakref from xml.etree import ElementTree @@ -418,7 +419,7 @@ def _wmts_images(self, wmts, layer, matrix_set_name, extent, return big_img, img_extent -class WFSGeometrySource(RasterSource): +class WFSGeometrySource(object): """Web Feature Service (WFS) retrieval for Cartopy.""" def __init__(self, service, features, getfeature_extra_kwargs=None): @@ -426,10 +427,11 @@ def __init__(self, service, features, getfeature_extra_kwargs=None): Args: * service: - The URL of a WFS, or an :class:`owslib.wfs.WebFeatureService` - instance. + The URL of a WFS, or an instance of + :class:`owslib.wfs.WebFeatureService`. * features: - The typename(s) of the features from the WFS to use. + The typename(s) of the features from the WFS that + will be retrieved and made available as geometries. Kwargs: @@ -448,7 +450,8 @@ def __init__(self, service, features, getfeature_extra_kwargs=None): # Populate an empty kwargs dict with something harmless. if getfeature_extra_kwargs is None: - getfeature_extra_kwargs = {'propertyname': ['*']} + #getfeature_extra_kwargs = {'propertyname': ['*']} + getfeature_extra_kwargs = {} if not features: raise ValueError('One or more features must be specified.') @@ -461,129 +464,88 @@ def __init__(self, service, features, getfeature_extra_kwargs=None): 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}" - #self._crs_urn_for_projection_id = {} - self._srs_for_projection_id = {} - - #def _crs_urn(self, projection): - # """ - # Confirm that the SRS of the features in the WFS response is recognised - # and that it lines up with the requested `projection`. -# -# """ -# key = id(projection) -# crs_urn = self._crs_urn_for_projection_id.get(key) -# if crs_urn is None: -# # Confirm the CRS URNs are recognised and match `projection`. -# for feature in self.features: -# crs_options = self.service.contents[feature].crsOptions -# for features_crs_urn in crs_options: -# if features_crs_urn in _URN_TO_CRS: -# this_crs = _URN_TO_CRS[features_crs_urn] -# if this_crs == projection: -# crs_urn = features_crs_urn -# break -# # Fail informatively if not. -# if crs_urn is None: -# raise ValueError('The projection {!r} was not convertible to ' -# 'a suitable WFS SRS.'.format(projection)) -# self._crs_urn_for_projection_id[key] = crs_urn -# return crs_urn - - #def validate_projection(self, projection): - # self._crs_urn(projection) - - #def _srs(self, projection): - # key = id(projection) - # srs = self._srs_for_projection_id.get(key) - # if srs is None: - # srs = _CRS_TO_OGC_SRS.get(projection) - # if srs is None: - # raise ValueError('The projection {!r} was not convertible to ' - # 'a suitable OGC SRS.'.format(projection)) - # for feature in self.features: - # if srs not in self.service.contents[feature].crsOptions: - # raise ValueError('The SRS {} is not a valid SRS for the ' - # '{!r} WFS feature.'.format(srs, feature)) - # self._srs_for_projection_id[key] = srs - # return srs - def default_projection(self): """ - Return a :class:`cartopy.crs.Projection` and corresponding - in which the WFS service can supply the requested features. + Return a :class:`cartopy.crs.Projection` in which the WFS + service can supply the requested features. """ - # Using first element in crsOptions (default) - #import pdb;pdb.set_trace() - default_srs = set(self.service.contents[feature].crsOptions[0] for - feature in self.features) - if len(default_srs) != 1: - ValueError('Failed to find a single common default SRS ' - 'across all features (typenames).') - else: - default_srs = unicode(list(default_srs)[0]) - - if default_srs not in _URN_TO_CRS: - raise ValueError('Unknown mapping from SRS/CRS_URN {!r} to ' - 'cartopy projection.'.format(default_srs)) + # 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 - return _URN_TO_CRS[default_srs] + if unicode(default_urn) not in _URN_TO_CRS: + raise ValueError('Unknown mapping from SRS/CRS_URN {!r} to ' + 'cartopy projection.'.format(default_srs)) - #for proj, srs in _CRS_TO_OGC_SRS.items(): - # missing = False - # for feature in self.features: - # if srs not in self.service.contents[feature].crsOptions: - # missing = True - # if not missing: - # break - #if missing: - # raise ValueError('The requested features are not available in ' - # 'any of the supported SRSs.') - #return proj + self._default_urn = default_urn + return _URN_TO_CRS[unicode(self._default_urn)] def fetch_geometries(self, projection, extent): """ - Return Shapely geometries of any Point, Linestring or LinearRing - geometries found in the WFS request. + Return any Point, Linestring or LinearRing geometries available + from the WFS that lie within the specified extent. - .. note:: + Args: - Kwargs passed to the get feature request are not guaranteed to - have an impact on the response data as the WFS server will not - necessarily be able to supply data that precisely matches the - request. In such cases the response will contain the default - data provided by the WFS server. + * 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. """ - #import pdb;pdb.set_trace() - min_x, max_x, min_y, max_y = extent #TODO - handle non-native proj. ie. convert extent to query wfs - may depend on - # on default srs. + 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, - #srsname=self._srs(projection), - doesn't nec. work - could try query params - bbox=(min_x, min_y, max_x, max_y), # - what proj????? + 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('Geoms in multiple SRSs') # DEBUG code + 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: - # Verify SRS. - geom_proj = _URN_TO_CRS[srs] - if geom_proj != projection: - raise ValueError('Geometries are not in requested proj - need to project...') - + 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): @@ -593,14 +555,17 @@ def _to_shapely_geoms(self, response): Args: - * response: + * 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 = [] - #import pdb;pdb.set_trace() tree = ElementTree.parse(response) for node in tree.findall('.//{ms}msGeometry'.format(ms=self.ms)): @@ -632,51 +597,58 @@ def _to_shapely_geoms(self, response): return geoms_by_srs def _find_polygon_coords(self, node, find_str): - """Find all coordinates data for a given Polygon `node`.""" + """ + 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') - # Assume we can either have nodes called `coordinates` or `coords`. - find_str = '{}coordinates'.format(self.gml) - if self._node_has_child(polygon, find_str): - points = polygon.findtext(find_str) - x, y = self._coordinates_to_points(points) - else: - x, y = [], [] - for coord in polygon.findall('{}coord'.format(self.gml)): + 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]) + data.append((feature_srs, x, y)) return data - def _node_has_child(self, node, find_str): + @staticmethod + def _node_has_child(node, find_str): """ - Determine whether `node` contains (at any sub-level), a node with name + Return whether `node` contains (at any sub-level), a node with name equal to `find_str`. """ - #TODO Can this be done with a generator?? - found = list(node.iterfind(find_str)) - return bool(found) - - def _coordinates_to_points(self, points): - """ - If our XML has a coordinates node then we'll get one string of - whitespace-separated points: - - "x,y x,y ... x,y" - - Parses this string to return lists of x and y points. - - """ - #TODO optimise this using itertools. - x = [] - y = [] - coords = points.strip().split(' ') - for coord in coords: - x_val, y_val = coord.split(',') - x.append(float(x_val)) - y.append(float(y_val)) - return x, y + element = node.find(find_str) + return element is not None From 3e11599f4e2177824b2026def4aab55f07439f39 Mon Sep 17 00:00:00 2001 From: Laura Dreyer Date: Fri, 7 Nov 2014 11:13:41 +0000 Subject: [PATCH 4/4] fixed _URN_TO_CRS for stereographic and included scaling_factor parameter --- lib/cartopy/crs.py | 4 +++- lib/cartopy/io/ogc_clients.py | 6 +++++- 2 files changed, 8 insertions(+), 2 deletions(-) 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/io/ogc_clients.py b/lib/cartopy/io/ogc_clients.py index f53b1c297..cca64d3af 100644 --- a/lib/cartopy/io/ogc_clients.py +++ b/lib/cartopy/io/ogc_clients.py @@ -74,7 +74,11 @@ _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.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)