diff --git a/lib/cartopy/io/ogc_clients.py b/lib/cartopy/io/ogc_clients.py index 39fdbd6ee..038d2589b 100644 --- a/lib/cartopy/io/ogc_clients.py +++ b/lib/cartopy/io/ogc_clients.py @@ -22,6 +22,7 @@ import warnings import weakref from xml.etree import ElementTree +from urllib.parse import urlparse from PIL import Image import numpy as np @@ -412,7 +413,12 @@ def find_projection(match_projection): matrix_sets = self.wmts.tilematrixsets tile_matrix_set = matrix_sets[tile_matrix_set_name] crs_urn = tile_matrix_set.crs - tms_crs = _URN_TO_CRS.get(crs_urn) + tms_crs = None + if crs_urn in _URN_TO_CRS: + tms_crs = _URN_TO_CRS.get(crs_urn) + elif ':EPSG:' in crs_urn: + epsg_num = crs_urn.split(':')[-1] + tms_crs = ccrs.epsg(int(epsg_num)) if tms_crs == match_projection: result = tile_matrix_set_name break @@ -422,6 +428,7 @@ def find_projection(match_projection): matrix_set_name = find_projection(target_projection) if matrix_set_name is None: # Search instead for a set in _any_ projection we can use. + # Nothing to do for EPSG for possible_projection in _URN_TO_CRS.values(): # Look for supported projections (in a preferred order). matrix_set_name = find_projection(possible_projection) @@ -445,8 +452,15 @@ def validate_projection(self, projection): def fetch_raster(self, projection, extent, target_resolution): matrix_set_name = self._matrix_set_name(projection) - wmts_projection = _URN_TO_CRS[ - self.wmts.tilematrixsets[matrix_set_name].crs] + crs_urn = self.wmts.tilematrixsets[matrix_set_name].crs + if crs_urn in _URN_TO_CRS: + wmts_projection = _URN_TO_CRS[crs_urn] + elif ':EPSG:' in crs_urn: + epsg_num = crs_urn.split(':')[-1] + wmts_projection = ccrs.epsg(int(epsg_num)) + else: + raise ValueError(f'Unknown coordinate reference system string:' + f' {crs_urn}') if wmts_projection == projection: wmts_extents = [extent] @@ -574,7 +588,24 @@ def _wmts_images(self, wmts, layer, matrix_set_name, extent, # Find which tile matrix has the appropriate resolution. tile_matrix_set = wmts.tilematrixsets[matrix_set_name] tile_matrices = tile_matrix_set.tilematrix.values() - meters_per_unit = METERS_PER_UNIT[tile_matrix_set.crs] + if tile_matrix_set.crs in METERS_PER_UNIT: + meters_per_unit = METERS_PER_UNIT[tile_matrix_set.crs] + elif ':EPSG:' in tile_matrix_set.crs: + epsg_num = tile_matrix_set.crs.split(':')[-1] + tms_crs = ccrs.epsg(int(epsg_num)) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + crs_dict = tms_crs.to_dict() + crs_units = '' + if 'units' in crs_dict: + crs_units = crs_dict['units'] + if crs_units != 'm': + raise ValueError('Only unit "m" implemented for' + ' EPSG projections.') + meters_per_unit = 1 + else: + raise ValueError(f'Unknown coordinate reference system string:' + f' {tile_matrix_set.crs}') tile_matrix = self._choose_matrix(tile_matrices, meters_per_unit, max_pixel_span) @@ -666,7 +697,22 @@ def __init__(self, service, features, getfeature_extra_kwargs=None): raise ImportError(_OWSLIB_REQUIRED) if isinstance(service, str): + # network location without port + # service is, for example, http://mapserver.gis.umn.edu/mapserver + # or https://agroenvgeo.data.inra.fr:443/geoserver/wfs + url = service + # network location is then: mapserver.gis.umn.edu + # or agroenvgeo.data.inra.fr:443 + url = urlparse(service).netloc + # network location without port: mapserver.gis.umn.edu + # or agroenvgeo.data.inra.fr + if ':' in url: # remove port if present + url = url[:url.rfind(':')] + self.url = url + # WebFeatureService of owslib service = WebFeatureService(service) + else: + self.url = '' if isinstance(features, str): features = [features] @@ -703,14 +749,22 @@ def default_projection(self): else: default_urn = default_urn.pop() - if str(default_urn) not in _URN_TO_CRS: + if (str(default_urn) not in _URN_TO_CRS) and ( + ":EPSG:" not in str(default_urn) + ): raise ValueError( - f'Unknown mapping from SRS/CRS_URN {default_urn!r} to ' - 'cartopy projection.') - + f"Unknown mapping from SRS/CRS_URN {default_urn!r} to " + "cartopy projection.") self._default_urn = default_urn - return _URN_TO_CRS[str(self._default_urn)] + if str(self._default_urn) in _URN_TO_CRS: + return _URN_TO_CRS[str(self._default_urn)] + elif ':EPSG:' in str(self._default_urn): + epsg_num = str(self._default_urn).split(':')[-1] + return ccrs.epsg(int(epsg_num)) + else: + raise ValueError(f'Unknown coordinate reference system:' + f' {str(self._default_urn)}') def fetch_geometries(self, projection, extent): """ @@ -757,12 +811,20 @@ def fetch_geometries(self, projection, extent): geom_proj = _URN_TO_CRS[srs] if geom_proj != projection: raise ValueError( - 'The geometries are not in expected projection. ' + f'The geometries are not in expected projection. ' f'Expected {projection!r}, got {geom_proj!r}.') + elif ':EPSG:' in srs: + epsg_num = srs.split(':')[-1] + geom_proj = ccrs.epsg(int(epsg_num)) + if geom_proj != projection: + raise ValueError( + f'The EPSG geometries are not in expected ' + f' projection. Expected {projection!r}, ' + f' got {geom_proj!r}.') else: warnings.warn( - 'Unable to verify matching projections due to ' - 'incomplete mappings from SRS identifiers to cartopy ' + f'Unable to verify matching projections due to ' + f'incomplete mappings from SRS identifiers to cartopy ' f'projections. The geometries have an SRS of {srs!r}.') return geoms @@ -807,6 +869,36 @@ def _to_shapely_geoms(self, response): data = self._find_polygon_coords(node, find_str) points_data.extend(data) + # Try to get other geometries from other servers than + # http://mapserver.gis.umn.edu/mapserver + if ((self.url and + (len(linear_rings_data) == 0) and + (len(linestrings_data) == 0) and + (len(points_data) == 0))): + for node in tree.iter(): + snode = str(node) + if self.url in snode: + s1 = snode.split()[1] + tag = s1[s1.find('}') + 1:-1] + if ('geom' in tag) or ('Geom' in tag): + # Find LinearRing geometries in our msGeometry node. + find_str = f'.//{_GML_NS}LinearRing' + 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 msGeometry node. + find_str = f'.//{_GML_NS}LineString' + 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 msGeometry node. + find_str = f'.//{_GML_NS}Point' + 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( diff --git a/lib/cartopy/tests/mpl/baseline_images/mpl/test_features/wfs_france.png b/lib/cartopy/tests/mpl/baseline_images/mpl/test_features/wfs_france.png new file mode 100644 index 000000000..cbfc8fc03 Binary files /dev/null and b/lib/cartopy/tests/mpl/baseline_images/mpl/test_features/wfs_france.png differ diff --git a/lib/cartopy/tests/mpl/test_features.py b/lib/cartopy/tests/mpl/test_features.py index 4a2e2b0cd..07a66f8a9 100644 --- a/lib/cartopy/tests/mpl/test_features.py +++ b/lib/cartopy/tests/mpl/test_features.py @@ -70,3 +70,17 @@ def test_wfs(): edgecolor='red') ax.add_feature(feature) return ax.figure + + +@pytest.mark.network +@pytest.mark.skipif(not _OWSLIB_AVAILABLE, reason='OWSLib is unavailable.') +@pytest.mark.xfail(raises=ParseError, + reason="Bad XML returned from the URL") +@pytest.mark.mpl_image_compare(filename='wfs_france.png') +def test_wfs_france(): + ax = plt.axes(projection=ccrs.epsg(2154)) + url = 'https://agroenvgeo.data.inra.fr:443/geoserver/wfs' + typename = 'collectif:t_ser_l93' + feature = cfeature.WFSFeature(url, typename, edgecolor='red') + ax.add_feature(feature) + return ax.figure