various fixes in attempt to process Berlin

geoplanar/gap.py

@@ -1,13 +1,13 @@
 #!/usr/bin/env python3
 #
+from collections import defaultdict
+
 import geopandas
 import numpy as np
 import pandas as pd
 import shapely
-from packaging.version import Version
-from collections import defaultdict
 from esda.shape import isoperimetric_quotient
-
+from packaging.version import Version
 
 __all__ = ["gaps", "fill_gaps", "snap"]
 
@@ -53,43 +53,49 @@ def gaps(gdf):
         crs=gdf.crs,
     )
     if GPD_GE_014:
-        poly_idx, _ = gdf.sindex.query(polygons, predicate="covers")
+        poly_idx, _ = gdf.representative_point().sindex.query(
+            polygons, predicate="contains"
+        )
     else:
-        poly_idx, _ = gdf.sindex.query_bulk(polygons, predicate="covers")
+        poly_idx, _ = gdf.representative_point().sindex.query_bulk(
+            polygons, predicate="cotains"
+        )
 
     return polygons.drop(poly_idx).reset_index(drop=True)
 
 
-def fill_gaps(gdf, gap_df=None, strategy='largest', inplace=False):
+def fill_gaps(gdf, gap_df=None, strategy="largest", inplace=False):
     """Fill gaps in a GeoDataFrame by merging them with neighboring polygons.
 
+    Note that as part of the gap filling heuristics, trim_overlaps is called internally.
+
     Parameters
     ----------
     gdf : GeoDataFrame
         A GeoDataFrame containing polygon or multipolygon geometries.
 
     gap_df : GeoDataFrame, optional
-        A GeoDataFrame containing the gaps to be filled. If None, gaps will be 
+        A GeoDataFrame containing the gaps to be filled. If None, gaps will be
         automatically detected within `gdf`.
 
     strategy : {'smallest', 'largest', 'compact', None}, default 'largest'
         Strategy to determine how gaps are merged with neighboring polygons:
           - 'smallest': Merge each gap with the smallest neighboring polygon.
           - 'largest' : Merge each gap with the largest neighboring polygon.
           - 'compact' : Merge each gap with the neighboring polygon that results in
-                        the new polygon having the highest compactness 
+                        the new polygon having the highest compactness
                         (isoperimetric quotient).
           - None      : Merge each gap with the first available neighboring polygon
                         (order is indeterminate).
 
     inplace : bool, default False
-        If True, modify the input GeoDataFrame in place. Otherwise, return a new 
+        If True, modify the input GeoDataFrame in place. Otherwise, return a new
         GeoDataFrame with the gaps filled.
 
     Returns
     -------
     GeoDataFrame or None
-        A new GeoDataFrame with gaps filled if `inplace` is False. Otherwise, 
+        A new GeoDataFrame with gaps filled if `inplace` is False. Otherwise,
         modifies `gdf` in place and returns None.
     """
     if gap_df is None:
@@ -99,35 +105,35 @@ def fill_gaps(gdf, gap_df=None, strategy='largest', inplace=False):
         gdf = gdf.copy()
 
     if not GPD_GE_014:
-        gap_idx, gdf_idx = gdf.sindex.query_bulk(
-            gap_df.geometry, predicate="intersects"
+        gap_idx, gdf_idx = gdf.boundary.sindex.query_bulk(
+            gap_df.boundary, predicate="overlaps"
         )
     else:
-        gap_idx, gdf_idx = gdf.sindex.query(gap_df.geometry, predicate="intersects")
+        gap_idx, gdf_idx = gdf.boundary.sindex.query(
+            gap_df.boundary, predicate="overlaps"
+        )
 
     to_merge = defaultdict(set)
 
     for g_ix in range(len(gap_df)):
         neighbors = gdf_idx[gap_idx == g_ix]
 
-        if strategy == 'compact':
+        if strategy == "compact":
             # Find the neighbor that results in the highest IQ
             gap_geom = shapely.make_valid(gap_df.geometry.iloc[g_ix])
             best_iq = -1
             best_neighbor = None
             neighbor_geometries = gdf.geometry.iloc[neighbors].apply(shapely.make_valid)
             for neighbor, neighbor_geom in zip(neighbors, neighbor_geometries):
-                combined_geom = shapely.union_all(
-                    [neighbor_geom, gap_geom]
-                )
+                combined_geom = shapely.union_all([neighbor_geom, gap_geom])
                 iq = isoperimetric_quotient(combined_geom)
                 if iq > best_iq:
                     best_iq = iq
                     best_neighbor = neighbor
             to_merge[best_neighbor].add(g_ix)
         elif strategy is None:  # don't care which polygon we attach cap to
             to_merge[gdf.index[neighbors[0]]].add(g_ix)
-        elif strategy == 'largest':
+        elif strategy == "largest":
             # Attach to the largest neighbor
             to_merge[gdf.iloc[neighbors].area.idxmax()].add(g_ix)
         else:
@@ -141,7 +147,10 @@ def fill_gaps(gdf, gap_df=None, strategy='largest', inplace=False):
                 [gdf.geometry.loc[k]] + [gap_df.geometry.iloc[i] for i in v]
             )
         )
-    gdf.loc[list(to_merge.keys()), gdf.geometry.name] = new_geom
+    if isinstance(gdf, geopandas.GeoDataFrame):
+        gdf.loc[list(to_merge.keys()), gdf.geometry.name] = new_geom
+    else:
+        gdf.loc[list(to_merge.keys())] = new_geom
 
     return gdf
 
@@ -182,6 +191,14 @@ def _snap(geometry, reference, threshold, segment_length):
     shapely.Polygon
         snapped geometry
     """
+    # special case - corner touch. If the area of intersection of buffers is below
+    # 110% of what a clean corner would cover, assume it is a corner that shall
+    # not be snapped.
+    if (
+        geometry.buffer(threshold).intersection(reference.buffer(threshold)).area
+        < threshold * 4 * 1.1
+    ):
+        return geometry
 
     # extract the shell and holes from the first geometry
     shell, holes = geometry.exterior, geometry.interiors
@@ -203,13 +220,17 @@ def _snap(geometry, reference, threshold, segment_length):
     # re-create the polygon with new coordinates and original holes and simplify
     # to remove any extra vertices.
     polygon = shapely.Polygon(coords, holes=holes)
-    simplified = shapely.make_valid(shapely.simplify(polygon, segment_length / 100))
+    simplified = shapely.make_valid(polygon)
     # the function may return invalid and make_valid may return non-polygons
     # the largest polygon is the most likely the one we want
     if simplified.geom_type != "Polygon":
         parts = _get_parts(simplified)
         simplified = parts[np.argmax(shapely.area(parts))]
 
+
+    # TODO: we need to add the same points to the exterior of reference, otherwise
+    # floating point precision may end in overlaps.
+
     return simplified
 
 
@@ -238,6 +259,10 @@ def snap(geometry, threshold):
     GeoSeries
         GeoSeries with snapped geometries
     """
+    # importing here to avoid circular imports
+    from .planar import fix_npe_edges
+
+    geometry = geometry.copy()
     if not GPD_GE_100:
         raise ImportError("geopandas 1.0.0 or higher is required.")
 
@@ -278,7 +303,10 @@ def snap(geometry, threshold):
             if previous_geom == geom:
                 new_geoms.append(
                     _snap(
-                        snapped_geom, ref, threshold=threshold, segment_length=threshold
+                        snapped_geom,
+                        ref,
+                        threshold=threshold,
+                        segment_length=threshold,
                     )
                 )
             else:
@@ -288,8 +316,14 @@ def snap(geometry, threshold):
                 new_geoms.append(snapped_geom)
                 previous_geom = geom
 
-        snapped = geometry.geometry.copy()
-        snapped.iloc[pairs_to_snap.get_level_values("source")] = new_geoms
-    else:
-        snapped = geometry.geometry.copy()
-    return snapped
+        if isinstance(geometry, geopandas.GeoDataFrame):
+            gs = geometry.geometry
+            gs.iloc[pairs_to_snap.get_level_values("source")] = new_geoms
+            geometry.geometry = gs
+        else:
+            geometry.iloc[pairs_to_snap.get_level_values("source")] = new_geoms
+
+    fixed_topology = fix_npe_edges(geometry)
+
+    fixed_topology.geometry = fixed_topology.simplify_coverage(threshold / 100)
+    return fixed_topology

geoplanar/overlap.py

@@ -4,8 +4,9 @@
 import geopandas
 import libpysal
 import numpy as np
-from packaging.version import Version
+import shapely
 from esda.shape import isoperimetric_quotient
+from packaging.version import Version
 
 __all__ = [
     "overlaps",
@@ -29,11 +30,21 @@ def overlaps(gdf):
     array-like: Pairs of indices with overlapping geometries.
     """
     if GPD_GE_014:
-        return gdf.sindex.query(gdf.geometry, predicate="overlaps")
-    return gdf.sindex.query_bulk(gdf.geometry, predicate="overlaps")
+        overlaps = gdf.sindex.query(gdf.geometry, predicate="overlaps")
+    else:
+        overlaps = gdf.sindex.query_bulk(gdf.geometry, predicate="overlaps")
+
+    # predicate occasionally marks pairs as overlapping if their intersection is a
+    # 1-d multipart geometry. Filter those out.
+    dimension = shapely.get_dimensions(
+        gdf.geometry.iloc[overlaps[0]].intersection(
+            gdf.geometry.iloc[overlaps[1]], align=False
+        )
+    )
+    return overlaps.T[dimension == 2].T
 
 
-def trim_overlaps(gdf, strategy='largest', inplace=False):
+def trim_overlaps(gdf, strategy="largest", inplace=False):
     """Trim overlapping polygons
 
     Note
@@ -54,17 +65,21 @@ def trim_overlaps(gdf, strategy='largest', inplace=False):
           - 'compact' : Trim the polygon yielding the most compact modified polygon.
                             (isoperimetric quotient).
           - None      : Trim either polygon non-deterministically but performantly.
-    
+
     Returns
     -------
 
     gdf: geodataframe with corrected geometries
 
     """
+
+    # TODO: difference may hit floating point precision and still result in overlappign
+    # geometries. In that case, fix_npe_edges fails.
+
     if GPD_GE_014:
-        intersections = gdf.sindex.query(gdf.geometry, predicate="intersects").T
+        intersections = gdf.sindex.query(gdf.geometry, predicate="overlaps").T
     else:
-        intersections = gdf.sindex.query_bulk(gdf.geometry, predicate="intersects").T
+        intersections = gdf.sindex.query_bulk(gdf.geometry, predicate="overlaps").T
 
     if not inplace:
         gdf = gdf.copy()
@@ -77,7 +92,7 @@ def trim_overlaps(gdf, strategy='largest', inplace=False):
                 left = gdf.geometry.iloc[i]
                 right = gdf.geometry.iloc[j]
                 gdf.iloc[j, geom_col_idx] = right.difference(left)
-    elif strategy=='largest':
+    elif strategy == "largest":
         for i, j in intersections:
             if i != j:
                 left = gdf.geometry.iloc[i]
@@ -86,7 +101,7 @@ def trim_overlaps(gdf, strategy='largest', inplace=False):
                     gdf.iloc[i, geom_col_idx] = left.difference(right)
                 else:
                     gdf.iloc[j, geom_col_idx] = right.difference(left)
-    elif strategy=='smallest':
+    elif strategy == "smallest":
         for i, j in intersections:
             if i != j:
                 left = gdf.geometry.iloc[i]
@@ -95,19 +110,19 @@ def trim_overlaps(gdf, strategy='largest', inplace=False):
                     gdf.iloc[i, geom_col_idx] = left.difference(right)
                 else:
                     gdf.iloc[j, geom_col_idx] = right.difference(left)
-    elif strategy=='compact':
-         for i, j in intersections:
-             if i != j:
-                 left = gdf.geometry.iloc[i]
-                 right = gdf.geometry.iloc[j]
-                 left_c = left.difference(right)
-                 right_c = right.difference(left)
-                 iq_left = isoperimetric_quotient(left_c)
-                 iq_right = isoperimetric_quotient(right_c)
-                 if iq_left > iq_right:  # trimming left is more compact than right
-                     gdf.iloc[i, geom_col_idx] = left_c
-                 else:
-                     gdf.iloc[j, geom_col_idx] = right_c
+    elif strategy == "compact":
+        for i, j in intersections:
+            if i != j:
+                left = gdf.geometry.iloc[i]
+                right = gdf.geometry.iloc[j]
+                left_c = left.difference(right)
+                right_c = right.difference(left)
+                iq_left = isoperimetric_quotient(left_c)
+                iq_right = isoperimetric_quotient(right_c)
+                if iq_left > iq_right:  # trimming left is more compact than right
+                    gdf.iloc[i, geom_col_idx] = left_c
+                else:
+                    gdf.iloc[j, geom_col_idx] = right_c
     return gdf