Merge pull request SciTools#2191 from mcuntz/epsg_to_ogc

Enhanced OGC reader

lib/cartopy/io/ogc_clients.py

@@ -19,6 +19,7 @@
 import collections
 import io
 import math
+from urllib.parse import urlparse
 import warnings
 import weakref
 from xml.etree import ElementTree
@@ -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,25 @@ 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))
+            # catch UserWarning from .to_dict(), because the resulting
+            # dictionary does not contain all information for all projections;
+            # need only 'units' here
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore")
+                crs_dict = tms_crs.to_dict()
+                crs_units = crs_dict.get('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 +698,15 @@ def __init__(self, service, features, getfeature_extra_kwargs=None):
             raise ImportError(_OWSLIB_REQUIRED)
 
         if isinstance(service, str):
+            # host name such as mapserver.gis.umn.edu or
+            # agroenvgeo.data.inra.fr from full address
+            # http://mapserver.gis.umn.edu/mapserver
+            # or https://agroenvgeo.data.inra.fr:443/geoserver/wfs
+            self.url = urlparse(service).hostname
+            # WebFeatureService of owslib
             service = WebFeatureService(service)
+        else:
+            self.url = ''
 
         if isinstance(features, str):
             features = [features]
@@ -703,14 +743,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 +805,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
 
@@ -788,24 +844,33 @@ def _to_shapely_geoms(self, response):
         points_data = []
         tree = ElementTree.parse(response)
 
-        for node in tree.findall(f'.//{_MAP_SERVER_NS}msGeometry'):
-            # 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)
+        # Get geometries from http://mapserver.gis.umn.edu/mapserver
+        # and other servers
+        for node in tree.iter():
+            snode = str(node)
+            if ((_MAP_SERVER_NS in snode) or
+                (self.url and (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:

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