Treatment of interior rings from transform between feature and projection domains

lib/cartopy/mpl/feature_artist.py

@@ -17,10 +17,16 @@
 import numpy as np
 import matplotlib.artist
 import matplotlib.collections
+import matplotlib.figure as figure
+
 
 import cartopy.mpl.patch as cpatch
+from cartopy.crs import PlateCarree
 from .style import merge as style_merge, finalize as style_finalize
 
+import shapely.ops as ops
+import shapely.geometry as sgeom
+
 
 class _GeomKey:
     """
@@ -130,6 +136,120 @@ def __init__(self, feature, **kwargs):
 
         self._feature = feature
 
+    def xy_samples(self, lon_sample, lat_sample):
+        feature_crs = self._feature.crs
+        if not isinstance(feature_crs, PlateCarree):
+            x_lims = feature_crs.x_limits
+            y_lims = feature_crs.y_limits
+
+            x_samples, x_sep = np.linspace(*x_lims, 400,
+                                           endpoint=True, retstep=True)
+
+            if 0. not in x_samples:
+                x_samples = np.insert(x_samples,
+                                      np.searchsorted(x_samples, 0.), 0.)
+
+            y_samples, y_sep = np.linspace(*y_lims, 400,
+                                           endpoint=True, retstep=True)
+
+            if 0. not in y_samples:
+                y_samples = np.insert(y_samples,
+                                      np.searchsorted(y_samples, 0.), 0.)
+
+            return x_samples, x_sep, y_samples, y_sep
+        else:
+            # We need to sample the central longitude and latitude, the
+            # antipode, and -1 * (central latitude)
+            x_sample, y_sample = lon_sample, lat_sample
+
+            x_sep = y_sep = 1
+
+            x_offset = x_sample - np.floor(x_sample)
+            x_samples = np.arange(-180. + x_offset, 180. + x_offset, x_sep)
+
+            y_offset = y_sample - np.floor(y_sample)
+            y_samples = np.arange(-90. + y_offset, 90. + y_offset, y_sep)
+
+            if x_offset == 0.:
+                x_samples = np.append(x_samples, 180.)
+            if y_offset == 0.:
+                y_samples = np.append(y_samples, 90.)
+
+            if -y_sample not in y_samples:
+                ind = np.searchsorted(y_samples, -y_sample)
+                y_samples = np.insert(y_samples, ind, -y_sample)
+
+            return x_samples, x_sep, y_samples, y_sep
+
+    def build_invalid_geom(self):
+        ax = self.axes
+        # the following may work after improvements to multipolygon creation
+        # from projection results
+#         feature_crs = self._feature.crs
+#         projection_crs = ax.projection
+#         projection_domain = projection_crs.domain
+#         feature_proj_domain = feature_crs.project_geometry(projection_domain,
+#         projection_crs)
+#         feature_domain = feature_crs.domain
+#         invalid_geom = feature_domain.difference(feature_proj_domain)
+        lon_sample = ax.projection._lon_sample
+        lat_sample = ax.projection._lat_sample
+
+        x_samples, x_sep, y_samples, y_sep = self.xy_samples(
+            lon_sample, lat_sample)
+        x_grid, y_grid = np.meshgrid(x_samples, y_samples)
+        xyz = ax.projection.transform_points(self._feature.crs, x_grid, y_grid)
+        x_proj_grid, y_proj_grid = xyz[:, :, 0], xyz[:, :, 1]
+        inds_x = ~np.isfinite(x_proj_grid)
+        inds_y = ~np.isfinite(y_proj_grid)
+
+        inds = inds_x | inds_y
+        inds_bin = inds.astype(int)
+
+        fig_ = figure.Figure()
+        ax_ = fig_.add_subplot()
+        # Method: create contours around the points that project to infinity
+        # Then contract the exterior rings and expand the interior rings to
+        # buffer the contours
+        # Create Polygons from the rings
+        fcontour = ax_.contourf(x_grid, y_grid,
+                                inds_bin, levels=[0.5, 1.5])
+
+        paths = fcontour.collections[0].get_paths()
+        invalid_geoms = []
+        for path in paths:
+            poly = path.to_polygons()
+            if poly:
+                # 0th path is polygon exterior
+                # Subsequent paths are interior rings
+                exterior = poly[0]
+                exterior = sgeom.LinearRing(exterior)
+                offset = max(x_sep, y_sep)
+                exterior = exterior.parallel_offset(offset, "right",
+                                                    join_style=3)
+                # Point of discussion:
+                # parallel_offset can output a MultiLineString from LinearRing
+                # input. Do we want to catch this behavior?
+                if isinstance(exterior, sgeom.MultiLineString):
+                    raise NotImplementedError
+                exterior.coords = list(exterior.coords)[::-1]
+                interiors = poly[1:]
+                interiors_shrunk = []
+                for interior in interiors:
+                    interior = sgeom.LinearRing(interior)
+                    interior_shrunk = interior.parallel_offset(offset, "right",
+                                                               join_style=3)
+                    # Again the possibility of a MultiLineString
+                    if not interior_shrunk.is_empty:
+                        interior_shrunk.coords = list(
+                                                interior_shrunk.coords)[::-1]
+                        interiors_shrunk.append(interior_shrunk)
+                invalid_geom = sgeom.Polygon(exterior, holes=interiors_shrunk)
+                invalid_geoms.append(invalid_geom)
+
+        invalid_geom = ops.unary_union(invalid_geoms)
+        return invalid_geom
+
     @matplotlib.artist.allow_rasterization
     def draw(self, renderer, *args, **kwargs):
         """
@@ -164,6 +284,7 @@ def draw(self, renderer, *args, **kwargs):
         # Project (if necessary) and convert geometries to matplotlib paths.
         stylised_paths = OrderedDict()
         key = ax.projection
+
         for geom in geoms:
             # As Shapely geometries cannot be relied upon to be
             # hashable, we have to use a WeakValueDictionary to manage
@@ -182,23 +303,35 @@ def draw(self, renderer, *args, **kwargs):
             mapping = FeatureArtist._geom_key_to_path_cache.setdefault(
                 geom_key, {})
             geom_paths = mapping.get(key)
-            if geom_paths is None:
-                if ax.projection != feature_crs:
-                    projected_geom = ax.projection.project_geometry(
-                        geom, feature_crs)
-                else:
-                    projected_geom = geom
-                geom_paths = cpatch.geos_to_path(projected_geom)
-                mapping[key] = geom_paths
-
             if not self._styler:
                 style = prepared_kwargs
             else:
-                # Unfreeze, then add the computed style, and then re-freeze.
-                style = style_merge(dict(prepared_kwargs), self._styler(geom))
+                # Unfreeze, then add the computed style,
+                # and then re-freeze.
+                style = style_merge(dict(prepared_kwargs),
+                                    self._styler(geom))
                 style = _freeze(style)
-
-            stylised_paths.setdefault(style, []).extend(geom_paths)
+            if geom_paths is None:
+                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
+                    else:
+                        invalid_geom = ax.projection.invalid_geoms[feature_crs]
+                    if not geom.is_valid:
+                        geom = geom.buffer(0)
+                    cleaned_geom = geom.difference(invalid_geom)
+                    if not cleaned_geom.is_empty:
+                        projected_geom = ax.projection.project_geometry(
+                                                         cleaned_geom,
+                                                         feature_crs)
+                        geom_paths = cpatch.geos_to_path(projected_geom)
+                        mapping[key] = geom_paths
+                else:
+                    geom_paths = cpatch.geos_to_path(geom)
+                    mapping[key] = geom_paths
+            if geom_paths:
+                stylised_paths.setdefault(style, []).extend(geom_paths)
 
         transform = ax.projection._as_mpl_transform(ax)
 
@@ -214,6 +347,5 @@ def draw(self, renderer, *args, **kwargs):
             c.set_clip_path(ax.patch)
             c.set_figure(ax.figure)
             c.draw(renderer)
-
         # n.b. matplotlib.collection.Collection.draw returns None
         return None

lib/cartopy/mpl/geoaxes.py

@@ -772,6 +772,7 @@ def get_extent(self, crs=None):
 
     def _get_extent_geom(self, crs=None):
         # Perform the calculations for get_extent(), which just repackages it.
+
         with self.hold_limits():
             if self.get_autoscale_on():
                 self.autoscale_view()
@@ -834,9 +835,9 @@ def set_extent(self, extents, crs=None):
         # plt.ylim - allowing users to set None for a minimum and/or
         # maximum value
         x1, x2, y1, y2 = extents
-        domain_in_crs = sgeom.polygon.LineString([[x1, y1], [x2, y1],
-                                                  [x2, y2], [x1, y2],
-                                                  [x1, y1]])
+        domain_in_crs = sgeom.Polygon([[x1, y1], [x2, y1],
+                                       [x2, y2], [x1, y2],
+                                       [x1, y1]])
 
         projected = None
 
@@ -850,9 +851,9 @@ def set_extent(self, extents, crs=None):
                            isinstance(self.projection, ccrs.PlateCarree)) or
                           crs == self.projection)
         if try_workaround:
-            boundary = self.projection.boundary
-            if boundary.equals(domain_in_crs):
-                projected = boundary
+            proj_domain = self.projection.domain
+            if proj_domain.equals(domain_in_crs):
+                projected = self.projection.boundary
 
         if projected is None:
             projected = self.projection.project_geometry(domain_in_crs, crs)

lib/cartopy/tests/mpl/test_features.py

@@ -12,6 +12,7 @@
 from cartopy.io.ogc_clients import _OWSLIB_AVAILABLE
 
 from cartopy.tests.mpl import ImageTesting
+from shapely.geometry import Polygon
 
 
 @pytest.mark.filterwarnings("ignore:Downloading")
@@ -65,3 +66,108 @@ def test_wfs():
     feature = cfeature.WFSFeature(url, typename,
                                   edgecolor='red')
     ax.add_feature(feature)
+
+
+@pytest.mark.natural_earth
+@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)))
+
+    fig = plt.figure(figsize=(20, 7))
+    proj1 = ccrs.LambertAzimuthalEqualArea(central_longitude=62.,
+                                           central_latitude=-50.)
+    ax1 = fig.add_subplot(1, 2, 1, projection=proj1)
+    ax1.add_geometries([lamaz_box], crs=proj1, color="tab:blue")
+
+    proj2 = ccrs.LambertAzimuthalEqualArea(central_longitude=-118.,
+                                           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()

lib/cartopy/tests/test_crs.py

@@ -305,3 +305,25 @@ def test_transform_points_outside_domain():
     # the same as the transform_points call with a length-1 array
     result = crs.transform_point(-120, 80, ccrs.PlateCarree())
     assert np.all(np.isnan(result))
+
+
+def test_lonlatgeom_attrs():
+    projs = []
+
+    def recurse_projs(cls):
+        for subcls in cls.__subclasses__():
+            if not subcls.__subclasses__():
+                projs.append(subcls)
+            else:
+                pass
+            recurse_projs(subcls)
+    recurse_projs(ccrs.Projection)
+
+    for proj in projs:
+        if proj.__name__ == "UTM":
+            proj_instance = proj(1)
+        else:
+            proj_instance = proj()
+        assert hasattr(proj_instance, "_lon_sample")
+        assert hasattr(proj_instance, "_lat_sample")
+        assert hasattr(proj_instance, "invalid_geoms")

lib/cartopy/tests/test_line_string.py

@@ -16,6 +16,9 @@
 
 class TestLineString:
     def test_out_of_bounds(self):
+
+        # This test is nonsense
+
         # Check that a line that is completely out of the map boundary produces
         # a valid LineString
         projection = ccrs.TransverseMercator(central_longitude=0, approx=True)

lib/cartopy/tests/test_linear_ring.py

@@ -44,6 +44,9 @@ def assert_intersection_with_boundary(segment_coords):
             assert_intersection_with_boundary(coords[-2:])
 
     def test_out_of_bounds(self):
+
+        # Also nonsense
+
         # Check that a ring that is completely out of the map boundary
         # produces an empty result.
         # XXX Check efficiency?