Update to smart_repair with some minor algorithmic changes to improve robustness, and updated Pandas syntax to remove FutureWarnings

Update to smart_repair with some minor algorithmic changes to improve robustness, and updated Pandas syntax to remove FutureWarnings

Author: jnclelland

maup/intersections.py

@@ -9,7 +9,7 @@
 @require_same_crs
 def intersections(sources, targets, output_type="geoseries", area_cutoff=None):
     """Computes all of the nonempty intersections between two sets of geometries.
-    By default, the returned :meth:`~geopandas.GeoSeries` will have a MultiIndex, where the
+    By default, the returned `~geopandas.GeoSeries` will have a MultiIndex, where the
     geometry at index *(i, j)* is the intersection of ``sources[i]`` and ``targets[j]``
     (if it is not empty).
     If output_type == "geodataframe", the return type is a range-indexed GeoDataFrame
@@ -20,8 +20,8 @@ def intersections(sources, targets, output_type="geoseries", area_cutoff=None):
     :param targets: geometries
     :type targets: :class:`~geopandas.GeoSeries` or :class:`~geopandas.GeoDataFrame`
     :rtype: :class:`~geopandas.GeoSeries`
-    :param area_cutoff: (optional) if provided, only return intersections with area 
-    greater than ``area_cutoff``
+    :param area_cutoff: (optional) if provided, only return intersections with
+        area greater than ``area_cutoff``
     :type area_cutoff: Number or None
     """
 

maup/repair.py

@@ -4,14 +4,8 @@
 import pandas
 
 from geopandas import GeoSeries
-from shapely.geometry import (
-    Polygon,
-    MultiPolygon,
-    LineString,
-    MultiLineString,
-    GeometryCollection,
-)
-from shapely import union_all
+from shapely.geometry import Polygon, MultiPolygon, LineString, MultiLineString, GeometryCollection
+from shapely.ops import unary_union
 
 from .adjacencies import adjacencies
 from .assign import assign_to_max
@@ -63,9 +57,8 @@ def trim_valid(value):
     """
     if isinstance(value, GeometryCollection):
         # List comprehension excluding non-Polygons
-        value = [
-            item for item in value.geoms if isinstance(item, (Polygon, MultiPolygon))
-        ]
+        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
@@ -77,11 +70,9 @@ def holes_of_union(geometries):
     if not all(
         isinstance(geometry, (Polygon, MultiPolygon)) for geometry in geometries
     ):
-        raise TypeError(
-            f"Must be a Polygon or MultiPolygon (got types {set([x.geom_type for x in geometries])})!"
-        )
+        raise TypeError(f"Must be a Polygon or MultiPolygon (got types {set([x.geom_type for x in geometries])})!")
 
-    union = union_all(geometries)
+    union = unary_union(geometries)
     series = holes(union)
     series.crs = geometries.crs
     return series
@@ -120,10 +111,7 @@ def close_gaps(geometries, relative_threshold=0.1, force_polygons=False):
     geometries = get_geometries(geometries)
     gaps = holes_of_union(geometries)
     return absorb_by_shared_perimeter(
-        gaps,
-        geometries,
-        relative_threshold=relative_threshold,
-        force_polygons=force_polygons,
+        gaps, geometries, relative_threshold=relative_threshold, force_polygons=force_polygons
     )
 
 
@@ -163,15 +151,10 @@ def resolve_overlaps(geometries, relative_threshold=0.1, force_polygons=False):
     to_remove = GeoSeries(
         pandas.concat([overlaps.droplevel(1), overlaps.droplevel(0)]), crs=overlaps.crs
     )
-    with_overlaps_removed = geometries.apply(
-        lambda x: x.difference(union_all(to_remove))
-    )
+    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,
-        force_polygons=force_polygons,
+        overlaps, with_overlaps_removed, relative_threshold=None, force_polygons=force_polygons
     )
 
 
@@ -189,9 +172,7 @@ def quick_repair(geometries, relative_threshold=0.1, force_polygons=False):
     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, force_polygons=force_polygons
-    )
+    return autorepair(geometries, relative_threshold=relative_threshold, force_polygons=force_polygons)
 
 
 def autorepair(geometries, relative_threshold=0.1, force_polygons=False):
@@ -213,28 +194,16 @@ def autorepair(geometries, relative_threshold=0.1, force_polygons=False):
 
     if force_polygons:
         geometries = make_valid_polygons(remove_repeated_vertices(geometries))
-        geometries = make_valid_polygons(
-            resolve_overlaps(
-                geometries,
-                relative_threshold=relative_threshold,
-                force_polygons=force_polygons,
-            )
-        )
-        geometries = make_valid_polygons(
-            close_gaps(
-                geometries,
-                relative_threshold=relative_threshold,
-                force_polygons=force_polygons,
-            )
-        )
+        geometries = make_valid_polygons(resolve_overlaps(geometries,
+                                                          relative_threshold=relative_threshold,
+                                                          force_polygons=force_polygons))
+        geometries = make_valid_polygons(close_gaps(geometries,
+                                                    relative_threshold=relative_threshold,
+                                                    force_polygons=force_polygons))
     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()
+        geometries = resolve_overlaps(geometries, relative_threshold=relative_threshold).make_valid()
+        geometries = close_gaps(geometries, relative_threshold=relative_threshold).make_valid()
 
     return geometries
 
@@ -244,9 +213,7 @@ def remove_repeated_vertices(geometries):
     Removes repeated vertices. Vertices are considered to be repeated if they
     appear consecutively, excluding the start and end points.
     """
-    return geometries.geometry.apply(
-        lambda x: apply_func_to_polygon_parts(x, dedup_vertices)
-    )
+    return geometries.geometry.apply(lambda x: apply_func_to_polygon_parts(x, dedup_vertices))
 
 
 def snap_to_grid(geometries, n=-7):
@@ -263,19 +230,16 @@ def crop_to(source, target):
     """
     Crops the source geometries to the target geometries.
     """
-    target_union = union_all(get_geometries(target))
-    cropped_geometries = get_geometries(source).apply(
-        lambda x: x.intersection(target_union)
-    )
+    target_union = unary_union(get_geometries(target))
+    cropped_geometries = get_geometries(source).apply(lambda x: x.intersection(target_union))
 
     if (cropped_geometries.area == 0).any():
-        warnings.warn(
-            "Some cropped geometries have zero area, likely due to\n"
-            + "large differences in the union of the geometries in your\n"
-            + "source and target shapefiles. This may become an issue\n"
-            + "when maupping.\n",
-            AreaCroppingWarning,
-        )
+        warnings.warn("Some cropped geometries have zero area, likely due to\n" +
+                      "large differences in the union of the geometries in your\n" +
+                      "source and target shapefiles. This may become an issue\n" +
+                      "when maupping.\n",
+                      AreaCroppingWarning
+                      )
 
     return cropped_geometries
 
@@ -291,18 +255,13 @@ def expand_to(source, target, force_polygons=False):
     else:
         geometries = get_geometries(source).make_valid()
 
-    source_union = union_all(geometries)
+    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
-    )
+    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,
-        force_polygons=force_polygons,
+        leftover_geometries, get_geometries(source), relative_threshold=None, force_polygons=force_polygons
     )
 
     return geometries
@@ -322,14 +281,14 @@ def doctor(source, target=None, silent=False, accept_holes=False):
     False. (Default is accept_holes = False.)
     """
     shapefiles = [source]
-    source_union = union_all(get_geometries(source))
+    source_union = unary_union(get_geometries(source))
 
     health_check = True
 
     if target is not None:
         shapefiles.append(target)
 
-        target_union = union_all(get_geometries(target))
+        target_union = unary_union(get_geometries(target))
         sym_area = target_union.symmetric_difference(source_union).area
 
         if sym_area != 0:
@@ -338,9 +297,7 @@ def doctor(source, target=None, silent=False, accept_holes=False):
             health_check = False
 
     for shp in shapefiles:
-        if not shp.geometry.apply(
-            lambda x: isinstance(x, (Polygon, MultiPolygon))
-        ).all():
+        if not shp.geometry.apply(lambda x: isinstance(x, (Polygon, MultiPolygon))).all():
             if silent is False:
                 print("Some rows do not have geometries.")
             health_check = False
@@ -389,9 +346,7 @@ def apply_func_to_polygon_parts(shape, func):
     elif isinstance(shape, MultiPolygon):
         return MultiPolygon([func(poly) for poly in shape.geoms])
     else:
-        raise TypeError(
-            f"Can only apply {func} to a Polygon or MultiPolygon (got {shape} with type {type(shape)})!"
-        )
+        raise TypeError(f"Can only apply {func} to a Polygon or MultiPolygon (got {shape} with type {type(shape)})!")
 
 
 def dedup_vertices(polygon):
@@ -426,31 +381,16 @@ def dedup_vertices(polygon):
 
 def snap_polygon_to_grid(polygon, n=-7):
     if len(polygon.interiors) == 0:
-        return Polygon(
-            [(round(x, -n), round(y, -n)) for x, y in polygon.exterior.coords]
-        )
+        return Polygon([(round(x, -n), round(y, -n)) for x, y in polygon.exterior.coords])
     else:
-        return Polygon(
-            [(round(x, -n), round(y, -n)) for x, y in polygon.exterior.coords],
-            holes=[
-                [(round(x, -n), round(y, -n)) for x, y in interior_ring.coords]
-                for interior_ring in polygon.interiors
-            ],
-        )
+        return Polygon([(round(x, -n), round(y, -n)) for x, y in polygon.exterior.coords], holes=[[(round(x, -n), round(y, -n)) for x, y in interior_ring.coords] for interior_ring in polygon.interiors])
 
 
 def snap_multilinestring_to_grid(multilinestring, n=-7):
     if multilinestring.geom_type == "LineString":
-        return LineString(
-            [(round(x, -n), round(y, -n)) for x, y in multilinestring.coords]
-        )
+        return LineString([(round(x, -n), round(y, -n)) for x, y in multilinestring.coords])
     elif multilinestring.geom_type == "MultiLineString":
-        return MultiLineString(
-            [
-                LineString([(round(x, -n), round(y, -n)) for x, y in linestring.coords])
-                for linestring in multilinestring.geoms
-            ]
-        )
+        return MultiLineString([LineString([(round(x, -n), round(y, -n)) for x, y in linestring.coords]) for linestring in multilinestring.geoms])
 
 
 def split_by_level(series, multiindex):
@@ -461,9 +401,7 @@ def split_by_level(series, multiindex):
 
 
 @require_same_crs
-def absorb_by_shared_perimeter(
-    sources, targets, relative_threshold=None, force_polygons=False
-):
+def absorb_by_shared_perimeter(sources, targets, relative_threshold=None, force_polygons=False):
     if len(sources) == 0:
         return targets
 
@@ -484,8 +422,7 @@ def absorb_by_shared_perimeter(
         assignment = assignment[under_threshold]
 
     sources_to_absorb = GeoSeries(
-        sources.groupby(assignment).apply(union_all),
-        crs=sources.crs,
+        sources.groupby(assignment).apply(unary_union), crs=sources.crs,
     )
 
     # Note that the following line produces a warning message when sources_to_absorb
@@ -497,7 +434,7 @@ def absorb_by_shared_perimeter(
 
     # This difference in indices is expected since not all target geometries may have sources
     # to absorb, so it would be nice to remove this warning.
-    result = targets.union(sources_to_absorb, align=True)
+    result = targets.union(sources_to_absorb)
 
     # The .union call only returns the targets who had a corresponding
     # source to absorb. Now we fill in all of the unchanged targets.