Removed debugging junk, added comments about algorithm, reduced duplicate code

Author: kdpenner

lib/cartopy/crs.py

@@ -312,84 +312,17 @@ def _project_multiline(self, geometry, src_crs):
         else:
             return []
 
-    def _project_multipolygon(self, geometry, src_crs):
-        # Project the polygon exterior/interior rings.
-        # Each source ring will result in either a ring, or one or more
-        # lines.
-        rings = []
-        multi_lines = []
-        for polygon in geometry:
-            if src_crs.is_geodetic():
-                is_ccw = True
-            else:
-                is_ccw = polygon.exterior.is_ccw
-            for src_ring in [polygon.exterior] + list(polygon.interiors):
-                p_rings, p_mline = self._project_linear_ring(src_ring, src_crs)
-                if p_rings:
-                    rings.extend(p_rings)
-                if len(p_mline) > 0:
-                    multi_lines.append(p_mline)
-
-            # Convert any lines to rings by attaching them to the boundary.
-            if multi_lines:
-                rings.extend(self._attach_lines_to_boundary(multi_lines,
-                                                            is_ccw))
-
-        # Resolve all the inside vs. outside rings, and convert to the
-        # final MultiPolygon.
-        exteriors, interiors = self._rings_to_multi_polygon(rings, is_ccw)
-        if exteriors or interiors:
-            if not interiors:
-                multi_poly = sgeom.MultiPolygon(exteriors)
-            elif not exteriors:
-                multi_poly_holes = ops.unary_union(interiors)
-                boundary_poly = self.domain
-                multi_poly = boundary_poly.difference(multi_poly_holes)
-                if isinstance(multi_poly, sgeom.Polygon):
-                    multi_poly = sgeom.MultiPolygon([multi_poly])
-            else:
-                multi_poly_bits = ops.unary_union(exteriors)
-                multi_poly_holes = ops.unary_union(interiors)
-                holes_envelope = multi_poly_holes.envelope
-                if holes_envelope.contains(multi_poly_bits.envelope):
-                    print("everything flipped")
-                    multi_poly = holes_envelope.difference(
-                            multi_poly_holes.difference(multi_poly_bits)
-                            )
-                else:
-                    multi_poly = multi_poly_bits.difference(multi_poly_holes)
-                if isinstance(multi_poly, sgeom.Polygon):
-                    multi_poly = sgeom.MultiPolygon([multi_poly])
-        else:
-            multi_poly = sgeom.MultiPolygon()
-#         geoms = []
-#         for geom in geometry.geoms:
-#             r = self._project_polygon(geom, src_crs)
-#             if r:
-#                 geoms.extend(r.geoms)
-#         if geoms:
-#             result = sgeom.MultiPolygon(geoms)
-#         else:
-#             result = sgeom.MultiPolygon()
-        return multi_poly
-
-    def _project_polygon(self, polygon, src_crs):
-        """
-        Return the projected polygon(s) derived from the given polygon.
-
-        """
+    def _make_rings_lines(self, polygon, rings, multi_lines, src_crs):
         # Determine orientation of polygon.
         # TODO: Consider checking the internal rings have the opposite
         # orientation to the external rings?
         if src_crs.is_geodetic():
             is_ccw = True
         else:
             is_ccw = polygon.exterior.is_ccw
-        # Project the polygon exterior/interior rings.
+        # Project the polygon exterior/interior rings
         # Each source ring will result in either a ring, or one or more
         # lines.
-        rings = []
-        multi_lines = []
         for src_ring in [polygon.exterior] + list(polygon.interiors):
             p_rings, p_mline = self._project_linear_ring(src_ring, src_crs)
             if p_rings:
@@ -400,7 +333,9 @@ def _project_polygon(self, polygon, src_crs):
         # Convert any lines to rings by attaching them to the boundary.
         if multi_lines:
             rings.extend(self._attach_lines_to_boundary(multi_lines, is_ccw))
+        return rings, multi_lines, is_ccw
 
+    def _construct_multipolygon(self, rings, is_ccw):
         # Resolve all the inside vs. outside rings, and convert to the
         # final MultiPolygon.
         exteriors, interiors = self._rings_to_multi_polygon(rings, is_ccw)
@@ -418,7 +353,6 @@ def _project_polygon(self, polygon, src_crs):
                 multi_poly_holes = ops.unary_union(interiors)
                 holes_envelope = multi_poly_holes.envelope
                 if holes_envelope.contains(multi_poly_bits.envelope):
-                    print("everything flipped")
                     multi_poly = holes_envelope.difference(
                             multi_poly_holes.difference(multi_poly_bits)
                             )
@@ -431,6 +365,31 @@ def _project_polygon(self, polygon, src_crs):
 
         return multi_poly
 
+    def _project_multipolygon(self, geometry, src_crs):
+        rings = []
+        multi_lines = []
+        polygons = []
+        for polygon in geometry:
+            rings, multi_lines, is_ccw = self._make_rings_lines(
+                                           polygon, rings, multi_lines, src_crs
+                                           )
+            polygons.extend(list(self._construct_multipolygon(rings, is_ccw)))
+            rings = []
+            multi_lines = []
+        return sgeom.MultiPolygon(polygons)
+
+    def _project_polygon(self, polygon, src_crs):
+        """
+        Return the projected polygon(s) derived from the given polygon.
+
+        """
+        rings = []
+        multi_lines = []
+        rings, multi_lines, is_ccw = self._make_rings_lines(
+            polygon, rings, multi_lines, src_crs
+            )
+        return self._construct_multipolygon(rings, is_ccw)
+
     def _attach_lines_to_boundary(self, multi_line_strings, is_ccw):
         """
         Return a list of LinearRings by attaching the ends of the given lines
@@ -627,6 +586,13 @@ def makes_valid_ring(line_string):
         return linear_rings
 
     def _rings_to_multi_polygon(self, rings, is_ccw):
+        # _project_linear_rings sometimes creates multiple exterior rings
+        # encircling the entire domain or almost the entire domain
+        # In this case, the method is:
+        # for each exterior ring find all interior rings contained by it
+        # create polygons
+        # if there are leftover interior rings, store them as Polygons to
+        # be added as holes by upstream multipolygon constructors
         exterior_rings = []
         interior_rings = []
         for ring in rings:
@@ -646,38 +612,35 @@ def _rings_to_multi_polygon(self, rings, is_ccw):
             for i, interior_ring in enumerate(interior_rings[:]):
                 if prep_polygon.contains(interior_ring):
                     holes.append(interior_ring)
-#                     interior_rings.remove(interior_ring)
                     interior_ring_flags.append(i)
                 elif polygon.crosses(interior_ring):
                     # Likely that we have an invalid geometry such as
                     # that from #509 or #537.
                     holes.append(interior_ring)
-#                     interior_rings.remove(interior_ring)
                     interior_ring_flags.append(i)
-#             if polygon.bounds == self.domain.bounds and not holes:
-#                 pass
             else:
                 polygon = sgeom.Polygon(exterior_ring, holes=holes)
                 polygon_bits.append(polygon)
         interior_rings = [ring for i, ring in enumerate(interior_rings)
                           if i not in interior_ring_flags]
 
         extra_holes = []
-        # Any left over "interior" rings need "inverting" with respect
-        # to the boundary.
         if interior_rings:
             boundary_poly = self.domain
-            x3, y3, x4, y4 = boundary_poly.bounds
-            bx = (x4 - x3) * 0.1
-            by = (y4 - y3) * 0.1
-            x3 -= bx
-            y3 -= by
-            x4 += bx
-            y4 += by
             polygon_bits = [poly.buffer(0) for poly in polygon_bits]
             if isinstance(self, PlateCarree):
-                # This should only run when called by
-                # _ViewClippedPathPatch.draw
+                # to set the initial plot extent _gen_axes_patch calls
+                # _ViewClippedPathPatch.draw, which calls
+                # _rings_to_multi_polygon with a geometry from project_geometry
+                #
+                # if set_extent() is not invoked this geometry represents the
+                # region in [-180, 180, -90, 90] that projects to infinity in
+                # the target projection
+                #
+                # the geometry sometimes is an interior ring with no exterior
+                # ring
+                #
+                # if this happens we have to invert the interior ring
                 polygon = sgeom.Polygon(interior_rings[0])
                 # Invert the polygon
                 polygon = boundary_poly.difference(polygon)
@@ -689,22 +652,6 @@ def _rings_to_multi_polygon(self, rings, is_ccw):
                         polygon = polygon.buffer(0)
                     if not polygon.is_empty:
                         extra_holes.append(polygon)
-#         if polygon_bits or extra_holes:
-#             if not extra_holes:
-#                 multi_poly = sgeom.MultiPolygon(polygon_bits)
-#             elif not polygon_bits:
-#                 multi_poly_holes = ops.unary_union(extra_holes)
-#                 multi_poly = boundary_poly.difference(multi_poly_holes)
-#                 if isinstance(multi_poly, sgeom.Polygon):
-#                     multi_poly = sgeom.MultiPolygon([multi_poly])
-#             else:
-#                 multi_poly_bits = ops.unary_union(polygon_bits)
-#                 multi_poly_holes = ops.unary_union(extra_holes)
-#                 multi_poly = multi_poly_bits.difference(multi_poly_holes)
-#                 if isinstance(multi_poly, sgeom.Polygon):
-#                     multi_poly = sgeom.MultiPolygon([multi_poly])
-#         else:
-#             multi_poly = sgeom.MultiPolygon()
         return polygon_bits, extra_holes
 
     def quick_vertices_transform(self, vertices, src_crs):

lib/cartopy/mpl/feature_artist.py

@@ -145,9 +145,17 @@ def xy_samples(self, lon_sample, lat_sample):
             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
@@ -200,6 +208,10 @@ def build_invalid_geom(self):
 
         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])
 
@@ -209,19 +221,25 @@ def build_invalid_geom(self):
             poly = path.to_polygons()
             if poly:
                 # 0th path is polygon exterior
-                # Following paths are interior rings
+                # 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]
@@ -305,10 +323,8 @@ def draw(self, renderer, *args, **kwargs):
                     if feature_crs not in ax.projection.invalid_geoms.keys():
                         invalid_geom = self.build_invalid_geom()
                         ax.projection.invalid_geoms[feature_crs] = invalid_geom
-                        print("ax != feature, missing")
                     else:
                         invalid_geom = ax.projection.invalid_geoms[feature_crs]
-                        print("ax != feature, present")
                     if not geom.is_valid:
                         geom = geom.buffer(0)
                     cleaned_geom = geom.difference(invalid_geom)
@@ -348,7 +364,7 @@ def draw(self, renderer, *args, **kwargs):
                                     )
                                 geom_paths.extend(
                                     cpatch.geos_to_path(projected_geom))
-                        if not key in mapping.keys():
+                        if key not in mapping.keys():
                             mapping[key] = geom_paths
                         else:
                             mapping[key].extend(geom_paths)

lib/cartopy/tests/mpl/baseline_images/mpl/test_features/UTM_all_zones.png

@@ -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,36 @@ def test_wfs():
     feature = cfeature.WFSFeature(url, typename,
                                   edgecolor='red')
     ax.add_feature(feature)
+
+
+@pytest.mark.natural_earth
+@ImageTesting(['UTM_all_zones'])
+def test_utm_all():
+    # this test fails because of padding (?) differences
+    # probably solved with an rcParam
+    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")
+
+
+@pytest.mark.natural_earth
+@ImageTesting(['nonplatecarree_with_projected_hole'])
+def test_nonpc_hole():
+    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")