Update repair.py

Incorporate changes from pull request #91.  New function is called "make_valid_polygons" rather than "make_valid" to avoid conflict with Shapely's make_valid function, and default value for force_polygons has been set to True.

Default value for force_polygons in autorepair/quick_repair is still set to False and make_valid is used throughout by default; if force_polygons is set to True, then make_valid_polygons is used throughout instead of make_valid.

maup/repair.py

@@ -4,7 +4,7 @@
 import pandas
 
 from geopandas import GeoSeries
-from shapely.geometry import Polygon, MultiPolygon, LineString, MultiLineString
+from shapely.geometry import Polygon, MultiPolygon, LineString, MultiLineString, GeometryCollection
 from shapely.ops import unary_union
 
 from .adjacencies import adjacencies
@@ -38,6 +38,32 @@ class AreaCroppingWarning(UserWarning):
     pass
 
 
+def make_valid_polygons(geometries, force_polygons=True):
+    """Extended make valid function with optional filter for non polygonal
+    data from the output."""
+    geometries = geometries.make_valid()
+    if force_polygons:
+        # Map function causes CRS information to be dropped.
+        crs = geometries.crs
+        geometries = geometries.map(trim_valid)
+        if crs is not None:
+            geometries = geometries.set_crs(crs)
+    return geometries
+
+
+def trim_valid(value):
+    """Ensures that the results from the "make_valid_polygons" method won't be
+    geometry collections with mixed types by filtering out non-polygons.
+    """
+    if isinstance(value, GeometryCollection):
+        # List comprehension excluding non-Polygons
+        value = [item for item in value.geoms
+                 if isinstance(item, (Polygon, MultiPolygon))]
+        # Re-aggregegating multiple Polygons into single MultiPolygon object.
+        value = value[0] if len(value) == 1 else MultiPolygon(value)
+    return value
+
+
 def holes_of_union(geometries):
     """Returns any holes in the union of the given geometries."""
     geometries = get_geometries(geometries)
@@ -68,7 +94,7 @@ def holes(geometry):
         raise TypeError("geometry must be a Polygon or MultiPolygon to have holes")
 
 
-def close_gaps(geometries, relative_threshold=0.1):
+def close_gaps(geometries, relative_threshold=0.1, force_polygons=False):
     """Closes gaps between geometries by assigning the hole to the polygon
     that shares the greatest perimeter with the hole.
 
@@ -78,15 +104,18 @@ def close_gaps(geometries, relative_threshold=0.1):
     gaps while closing the tiny gaps that can occur as artifacts of
     geospatial operations. Set `relative_threshold=None` to attempt close all
     gaps. Due to floating point precision issues, all gaps may not be closed.
+
+    Optional "force_polygons" argument included to apply an automatic filter
+    non-polygonal fragments during operation.
     """
     geometries = get_geometries(geometries)
     gaps = holes_of_union(geometries)
     return absorb_by_shared_perimeter(
-        gaps, geometries, relative_threshold=relative_threshold
+        gaps, geometries, relative_threshold=relative_threshold, force_polygons=force_polygons
     )
 
 
-def resolve_overlaps(geometries, relative_threshold=0.1):
+def resolve_overlaps(geometries, relative_threshold=0.1, force_polygons=False):
     """For any pair of overlapping geometries, assigns the overlapping area to the
     geometry that shares the greatest perimeter with the overlap. Returns the GeoSeries
     of geometries, which will have no overlaps.
@@ -98,10 +127,16 @@ def resolve_overlaps(geometries, relative_threshold=0.1):
     that might indicate deeper issues and should be handled on a case-by-case
     basis. Set `relative_threshold=None` to attempt to resolve all overlaps. Due
     to floating point precision issues, all overlaps may not be resolved.
+
+    Optional "force_polygons" argument included to apply an automatic filter
+    non-polygonal fragments during operation.
     """
     geometries = get_geometries(geometries)
     inters = adjacencies(geometries, warn_for_islands=False, warn_for_overlaps=False)
-    overlaps = inters[inters.area > 0].make_valid()
+    if force_polygons:
+        overlaps = make_valid_polygons(inters[inters.area > 0])
+    else:
+        overlaps = inters[inters.area > 0].make_valid()
 
     if relative_threshold is not None:
         left_areas, right_areas = split_by_level(geometries.area, overlaps.index)
@@ -119,11 +154,11 @@ def resolve_overlaps(geometries, relative_threshold=0.1):
     with_overlaps_removed = geometries.apply(lambda x: x.difference(unary_union(to_remove)))
 
     return absorb_by_shared_perimeter(
-        overlaps, with_overlaps_removed, relative_threshold=None
+        overlaps, with_overlaps_removed, relative_threshold=None, force_polygons=force_polygons
     )
 
 
-def quick_repair(geometries, relative_threshold=0.1):
+def quick_repair(geometries, relative_threshold=0.1, force_polygons=False):
     """
     New name for autorepair function from Maup 1.x.
     Uses simplistic algorithms to repair most gaps and overlaps.
@@ -137,12 +172,14 @@ def quick_repair(geometries, relative_threshold=0.1):
     For a more careful repair that takes adjacencies and higher-order overlaps
     between geometries into account, consider using smart_repair instead.
     """
-    return autorepair(geometries, relative_threshold=relative_threshold)
+    return autorepair(geometries, relative_threshold=relative_threshold, force_polygons=force_polygons)
 
 
-def autorepair(geometries, relative_threshold=0.1):
+def autorepair(geometries, relative_threshold=0.1, force_polygons=False):
     """
-    Uses simplistic algorithms to repair most gaps and overlaps.
+    Uses simplistic algorithms to repair most gaps and overlaps. Additional
+    optional argument provided that can drop any fragmented line segments that
+    may occur while using shapely's make valid function. Should work by default.
 
     The default `relative_threshold` is `0.1`. This default is chosen to include
     tiny overlaps that can be safely auto-fixed while preserving major overlaps
@@ -155,9 +192,16 @@ def autorepair(geometries, relative_threshold=0.1):
     """
     geometries = get_geometries(geometries)
 
-    geometries = remove_repeated_vertices(geometries).make_valid()
-    geometries = resolve_overlaps(geometries, relative_threshold=relative_threshold).make_valid()
-    geometries = close_gaps(geometries, relative_threshold=relative_threshold).make_valid()
+    if force_polygons:
+        geometries = make_valid_polygons(remove_repeated_vertices(geometries))
+        geometries = make_valid_polygons(resolve_overlaps(geometries,
+                                                          relative_threshold=relative_threshold))
+        geometries = make_valid_polygons(close_gaps(geometries,
+                                                    relative_threshold=relative_threshold))
+    else:
+        geometries = remove_repeated_vertices(geometries).make_valid()
+        geometries = resolve_overlaps(geometries, relative_threshold=relative_threshold).make_valid()
+        geometries = close_gaps(geometries, relative_threshold=relative_threshold).make_valid()
 
     return geometries
 
@@ -199,19 +243,23 @@ def crop_to(source, target):
 
 
 @require_same_crs
-def expand_to(source, target):
+def expand_to(source, target, force_polygons=False):
     """
     Expands the source geometries to the target geometries.
     """
-    geometries = get_geometries(source).make_valid()
+
+    if force_polygons:
+        geometries = make_valid_polygons(get_geometries(source))
+    else:
+        geometries = get_geometries(source).make_valid()
 
     source_union = unary_union(geometries)
 
     leftover_geometries = get_geometries(target).apply(lambda x: x - source_union)
     leftover_geometries = leftover_geometries[~leftover_geometries.is_empty].explode(index_parts=False)
 
     geometries = absorb_by_shared_perimeter(
-        leftover_geometries, get_geometries(source), relative_threshold=None
+        leftover_geometries, get_geometries(source), relative_threshold=None, force_polygons=force_polygons
     )
 
     return geometries
@@ -351,14 +399,17 @@ def split_by_level(series, multiindex):
 
 
 @require_same_crs
-def absorb_by_shared_perimeter(sources, targets, relative_threshold=None):
+def absorb_by_shared_perimeter(sources, targets, relative_threshold=None, force_polygons=False):
     if len(sources) == 0:
         return targets
 
     if len(targets) == 0:
         raise IndexError("targets must be nonempty")
 
-    inters = intersections(sources, targets, area_cutoff=None).make_valid()
+    if force_polygons:
+        inters = make_valid_polygons(intersections(sources, targets, area_cutoff=None))
+    else:
+        inters = intersections(sources, targets, area_cutoff=None).make_valid()
 
     assignment = assign_to_max(inters.length)