Added tests

lib/cartopy/mpl/feature_artist.py

@@ -285,8 +285,6 @@ def draw(self, renderer, *args, **kwargs):
         stylised_paths = OrderedDict()
         key = ax.projection
 
-        skip_geoms = ["LineString", "LinearRing", "MultiLineString"]
-
         for geom in geoms:
             # As Shapely geometries cannot be relied upon to be
             # hashable, we have to use a WeakValueDictionary to manage
@@ -314,12 +312,7 @@ def draw(self, renderer, *args, **kwargs):
                                     self._styler(geom))
                 style = _freeze(style)
             if geom_paths is None:
-                # Shapely can't represent the difference between a line and a
-                # polygon. If all geoms are lines, bypass building the polygon
-                # containing points projecting to infinity and use different
-                # masking method.
-                if (ax.projection != feature_crs) & \
-                   (geom.geom_type not in skip_geoms):
+                if ax.projection != feature_crs:
                     if feature_crs not in ax.projection.invalid_geoms.keys():
                         invalid_geom = self.build_invalid_geom()
                         ax.projection.invalid_geoms[feature_crs] = invalid_geom
@@ -334,40 +327,6 @@ def draw(self, renderer, *args, **kwargs):
                                                          feature_crs)
                         geom_paths = cpatch.geos_to_path(projected_geom)
                         mapping[key] = geom_paths
-
-                elif (ax.projection != feature_crs) & \
-                     (geom.geom_type in skip_geoms):
-                    if isinstance(geom, sgeom.LineString):
-                        geom = sgeom.MultiLineString([geom])
-                    elif isinstance(geom, sgeom.LinearRing):
-                        geom = [geom]
-                    for subgeom in geom:
-                        x, y = np.array(subgeom.xy)
-                        xyz = ax.projection.transform_points(feature_crs, x, y)
-                        x_proj, y_proj = xyz[:, 0], xyz[:, 1]
-
-                        inds_x = ~np.isfinite(x_proj)
-                        inds_y = ~np.isfinite(y_proj)
-                        inds = inds_x | inds_y
-
-                        x_proj = np.ma.masked_array(x_proj, mask=inds)
-                        clump_slices = np.ma.clump_unmasked(x_proj)
-                        geom_paths = []
-                        for clump_slice in clump_slices:
-                            if clump_slice.stop - clump_slice.start > 1:
-                                xt = x[clump_slice]
-                                yt = y[clump_slice]
-                                cleaned_geom = sgeom.LineString(zip(xt, yt))
-                                projected_geom = \
-                                    ax.projection.project_geometry(
-                                            cleaned_geom, feature_crs
-                                    )
-                                geom_paths.extend(
-                                    cpatch.geos_to_path(projected_geom))
-                        if key not in mapping.keys():
-                            mapping[key] = geom_paths
-                        else:
-                            mapping[key].extend(geom_paths)
                 else:
                     geom_paths = cpatch.geos_to_path(geom)
                     mapping[key] = geom_paths

lib/cartopy/tests/mpl/test_features.py

@@ -69,22 +69,41 @@ def test_wfs():
 
 
 @pytest.mark.natural_earth
-@ImageTesting(['UTM_all_zones'])
-def test_utm_all():
+@ImageTesting(['UTM_coastlines'])
+def test_utm_coastlines():
     # this test fails because of padding (?) differences
     # probably solved with an rcParam
+    # if a feature contains points that project to infinity the plot will have
+    # artifacts
+    # this tests coastlines and should replace the example in the gallery
     zones = range(1, 61)
     fig = plt.figure(figsize=(18, 6))
     for zone in zones:
         ax = fig.add_subplot(1, len(zones), zone,
                              projection=ccrs.UTM(zone=zone,
                              southern_hemisphere=True))
         ax.add_feature(cfeature.LAND, facecolor="tab:blue")
+        ax.coastlines()
+
+
+@pytest.mark.natural_earth
+@ImageTesting(['tmerc_oceans'])
+def test_tmerc_oceans():
+    # if a feature contains points that project to infinity the plot will have
+    # artifacts
+    # this tests polygons
+    fig = plt.figure(figsize=(10, 7))
+    proj = ccrs.TransverseMercator(central_longitude=18.14159, approx=False)
+    ax = fig.add_subplot(1, 1, 1, projection=proj)
+    ax.add_feature(cfeature.OCEAN.with_scale("10m"), facecolor="tab:blue")
+    ax.coastlines()
 
 
 @pytest.mark.natural_earth
 @ImageTesting(['nonplatecarree_with_projected_hole'])
 def test_nonpc_hole():
+    # tests the difference algo if a feature in a non-Plate Carree CRS is
+    # passed
     lamaz_box = Polygon(((8.5e5, 1.1e6), (-8.5e5, 1.1e6),
                          (-5.6e5, -1.2e6), (5.6e5, -1.2e6),
                          (8.5e5, 1.1e6)))
@@ -99,3 +118,56 @@ def test_nonpc_hole():
                                            central_latitude=50.)
     ax2 = fig.add_subplot(1, 2, 2, projection=proj2)
     ax2.add_geometries([lamaz_box], crs=proj1, color="tab:blue")
+
+
+@pytest.mark.natural_earth
+@ImageTesting(['azimuthal_interior_rings'])
+def test_azi_interior_rings():
+    # in cartopy <= 0.18 the following test produces an entirely blue plot
+    # because the algo inverts ~6800 interior rings wrt the boundary instead of
+    # treating them as holes
+    fig = plt.figure(figsize=(10, 7))
+    proj = ccrs.AzimuthalEquidistant(central_longitude=15, central_latitude=62)
+    ax = fig.add_subplot(1, 1, 1, projection=proj)
+    ax.add_feature(cfeature.OCEAN.with_scale("10m"), facecolor="tab:blue")
+    ax.coastlines()
+
+
+@pytest.mark.natural_earth
+@ImageTesting(['lambert_equalarea_oceans'])
+def test_lam_eq_oceans():
+    # Note: this test doesn't crash in 0.18 and does in 7e077e589d2a
+    # in 0.18 it produces an entirely blue plot
+    fig = plt.figure(figsize=(10, 7))
+    proj = ccrs.LambertAzimuthalEqualArea(central_longitude=-62,
+                                          central_latitude=-15)
+    ax = fig.add_subplot(1, 1, 1, projection=proj)
+    ax.add_feature(cfeature.OCEAN.with_scale("10m"), facecolor="tab:blue")
+    ax.coastlines()
+
+
+@pytest.mark.natural_earth
+@ImageTesting(['lambert_exterior_rings'])
+def test_azi_exterior_rings():
+    # in cartopy <= 0.18 the following test produces an entirely blue plot
+    # because the algo doesn't subtract all interior rings from multiple
+    # exterior rings that span the entire domain
+    fig = plt.figure(figsize=(10, 7))
+    proj = ccrs.LambertAzimuthalEqualArea(central_longitude=-62,
+                                          central_latitude=15)
+    ax = fig.add_subplot(1, 1, 1, projection=proj)
+    ax.add_feature(cfeature.OCEAN.with_scale("10m"), facecolor="tab:blue")
+    ax.coastlines()
+
+
+@pytest.mark.natural_earth
+@ImageTesting(['lambert_ring_flip'])
+def test_lam_ring_flip():
+    # sometimes interior rings become exterior rings; test the containment
+    # logic
+    fig = plt.figure(figsize=(10, 7))
+    proj = ccrs.LambertAzimuthalEqualArea(central_longitude=62,
+                                          central_latitude=-50)
+    ax = fig.add_subplot(1, 1, 1, projection=proj)
+    ax.add_feature(cfeature.OCEAN.with_scale("10m"), facecolor="tab:blue")
+    ax.coastlines()