Enhanced OGC reader

Added support for general EPSG projections in OGC clients WMSRasterSource, WMTSRasterSource, and WFSGeometrySource.

Try to determine geometries from other map servers than the hard-coded mapserver.gis.umn.edu.

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(

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