Merge pull request #60 from sjsrey/compact

Add option for maximizing compactness when filling gaps

Author: sjsrey

ci/310-latest.yaml

@@ -4,6 +4,7 @@ channels:
 dependencies:
   - python=3.10
   - libpysal
+  - esda
   - geopandas
   - packaging
   - pytest

ci/310-oldest.yaml

@@ -7,6 +7,7 @@ dependencies:
   - geopandas=0.10.2
   - scipy=1.8
   - pandas=1.5
+  - esda
   - packaging
   - pytest
   - pytest-cov

ci/311-latest.yaml

@@ -4,6 +4,7 @@ channels:
 dependencies:
   - python=3.11
   - libpysal
+  - esda
   - geopandas
   - packaging
   - pytest

ci/312-dev.yaml

@@ -15,4 +15,5 @@ dependencies:
       - --pre --index-url https://pypi.anaconda.org/scientific-python-nightly-wheels/simple --extra-index-url https://pypi.org/simple
       - pandas
       - git+https://github.com/geopandas/geopandas.git@main
-      - git+https://github.com/pysal/libpysal.git@main
+      - git+https://github.com/pysal/libpysal.git@main
+      - git+https://github.com/pysal/esda.git@main

ci/312-latest.yaml

@@ -4,6 +4,7 @@ channels:
 dependencies:
   - python=3.12
   - libpysal
+  - esda
   - packaging
   - pytest
   - pytest-cov

geoplanar/gap.py

@@ -5,6 +5,9 @@
 import pandas as pd
 import shapely
 from packaging.version import Version
+from collections import defaultdict
+from esda.shape import isoperimetric_quotient
+
 
 __all__ = ["gaps", "fill_gaps", "snap"]
 
@@ -57,51 +60,38 @@ def gaps(gdf):
     return polygons.drop(poly_idx).reset_index(drop=True)
 
 
-def fill_gaps(gdf, gap_df=None, largest=True, inplace=False):
-    """Fill any gaps in a geodataframe.
+def fill_gaps(gdf, gap_df=None, strategy='largest', inplace=False):
+    """Fill gaps in a GeoDataFrame by merging them with neighboring polygons.
 
     Parameters
     ----------
-
-    gdf :  GeoDataFrame with polygon (multipolygon) GeoSeries
-
-
-    gap_df:  GeoDataFrame with gaps to fill
-             If None, gaps will be determined
-
-    largest: boolean (Default: True)
-          Merge each gap with its largest (True), or smallest (False) neighbor.
-          If None, merge with any neighbor non-deterministically but performantly.
-
-    inplace: boolean (default: False)
-          Change the geoseries of current dataframe
-
+    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 
+        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 
+                        (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 
+        GeoDataFrame with the gaps filled.
 
     Returns
     -------
-
-    _gaps : GeoDataFrame with gap polygons
-
-    Examples
-    --------
-    >>> p1 = box(0,0,10,10)
-    >>> p2 = Polygon([(10,10), (12,8), (10,6), (12,4), (10,2), (20,5)])
-    >>> gdf = geopandas.GeoDataFrame(geometry=[p1,p2])
-    >>> gdf.area
-    0    100.0
-    1     32.0
-    dtype: float64
-    >>> h = geoplanar.gaps(gdf)
-    >>> h.area
-    array([4., 4.])
-    >>> gdf1 = geoplanar.fill_gaps(gdf)
-    >>> gdf1.area
-    0    108.0
-    1     32.0
-    dtype: float64
+    GeoDataFrame or None
+        A new GeoDataFrame with gaps filled if `inplace` is False. Otherwise, 
+        modifies `gdf` in place and returns None.
     """
-    from collections import defaultdict
-
     if gap_df is None:
         gap_df = gaps(gdf)
 
@@ -116,13 +106,32 @@ def fill_gaps(gdf, gap_df=None, largest=True, inplace=False):
         gap_idx, gdf_idx = gdf.sindex.query(gap_df.geometry, predicate="intersects")
 
     to_merge = defaultdict(set)
+
     for g_ix in range(len(gap_df)):
         neighbors = gdf_idx[gap_idx == g_ix]
-        if largest is None:  # don't care which polygon is a gap attached to
+
+        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]
+                )
+                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 largest:
+        elif strategy == 'largest':
+            # Attach to the largest neighbor
             to_merge[gdf.iloc[neighbors].area.idxmax()].add(g_ix)
         else:
+            # Attach to the smallest neighbor
             to_merge[gdf.iloc[neighbors].area.idxmin()].add(g_ix)
 
     new_geom = []

geoplanar/overlap.py

@@ -5,6 +5,7 @@
 import libpysal
 import numpy as np
 from packaging.version import Version
+from esda.shape import isoperimetric_quotient
 
 __all__ = [
     "overlaps",
@@ -32,7 +33,7 @@ def overlaps(gdf):
     return gdf.sindex.query_bulk(gdf.geometry, predicate="overlaps")
 
 
-def trim_overlaps(gdf, largest=True, inplace=False):
+def trim_overlaps(gdf, strategy='largest', inplace=False):
     """Trim overlapping polygons
 
     Note
@@ -46,11 +47,14 @@ def trim_overlaps(gdf, largest=True, inplace=False):
 
     gdf:  geodataframe with polygon geometries
 
-    largest: boolean (Default: True)
-            True: trim the larger of the pair of overlapping polygons,
-            False: trim the smaller polygon.
-            If None, trim either polygon non-deterministically but performantly.
-
+    strategy : {'smallest', 'largest', 'compact', None}, default 'largest'
+        Strategy to determine which polygon to trim.
+          - 'smallest': Trim the smallest polygon.
+          - 'largest' : Trim the largest polygon.
+          - 'compact' : Trim the polygon yielding the most compact modified polygon.
+                            (isoperimetric quotient).
+          - None      : Trim either polygon non-deterministically but performantly.
+    
     Returns
     -------
 
@@ -67,35 +71,43 @@ def trim_overlaps(gdf, largest=True, inplace=False):
 
     geom_col_idx = gdf.columns.get_loc(gdf.geometry.name)
 
-    if largest is None:  # don't care which polygon to trim
+    if strategy is None:  # don't care which polygon to trim
         for i, j in intersections:
             if i != j:
                 left = gdf.geometry.iloc[i]
                 right = gdf.geometry.iloc[j]
-                right = gdf.geometry.iloc[j].difference(gdf.geometry.iloc[i])
-                gdf.iloc[j, geom_col_idx] = right
-    elif largest:
+                gdf.iloc[j, geom_col_idx] = right.difference(left)
+    elif strategy=='largest':
         for i, j in intersections:
             if i != j:
                 left = gdf.geometry.iloc[i]
                 right = gdf.geometry.iloc[j]
-                if left.area < right.area:
-                    right = gdf.geometry.iloc[j].difference(gdf.geometry.iloc[i])
-                    gdf.iloc[j, geom_col_idx] = right
+                if left.area > right.area:  # trim left
+                    gdf.iloc[i, geom_col_idx] = left.difference(right)
                 else:
-                    left = gdf.geometry.iloc[i].difference(gdf.geometry.iloc[j])
-                    gdf.iloc[i, geom_col_idx] = left
-    else:
+                    gdf.iloc[j, geom_col_idx] = right.difference(left)
+    elif strategy=='smallest':
         for i, j in intersections:
             if i != j:
                 left = gdf.geometry.iloc[i]
                 right = gdf.geometry.iloc[j]
-                if left.area > right.area:
-                    right = gdf.geometry.iloc[j].difference(gdf.geometry.iloc[i])
-                    gdf.iloc[j, geom_col_idx] = right
+                if left.area < right.area:  # trim left
+                    gdf.iloc[i, geom_col_idx] = left.difference(right)
                 else:
-                    left = gdf.geometry.iloc[i].difference(gdf.geometry.iloc[j])
-                    gdf.iloc[i, geom_col_idx] = left
+                    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
     return gdf
 
 

geoplanar/tests/test_compact.py

@@ -0,0 +1,43 @@
+import pytest
+from shapely.geometry import box
+import geopandas as gpd
+from geoplanar import fill_gaps  # Updated module name
+from packaging.version import Version
+import pytest
+
+@pytest.mark.skipif(Version(gpd.__version__) == Version("0.10.2"), reason="Missing pygeos")    
+def test_fill_gaps_compact_visual_case():
+    """Test fill_gaps behavior with compact=True and compact=False using specific geometries."""
+    # Define four polygons with a visible gap
+    p1 = box(0, 0, 10, 10)     # Left rectangle
+    p2 = box(10, 0, 40, 2)     # Bottom thin rectangle
+    p3 = box(10, 3, 40, 10)    # Top thin rectangle
+    p4 = box(40, 0, 100, 10)   # Right rectangle
+
+    # Create GeoDataFrame
+    gdf = gpd.GeoDataFrame(geometry=[p1, p2, p3, p4])
+
+    # Fill gaps with compact
+    filled_gdf_compact = fill_gaps(gdf, strategy='compact')
+    compact_geom_count = len(filled_gdf_compact)
+
+    # Fill gaps with largest (default)
+    filled_gdf_default = fill_gaps(gdf, strategy='largest')
+    default_geom_count = len(filled_gdf_default)
+
+    # Assert the number of geometries is the same after filling gaps
+    assert compact_geom_count == default_geom_count, (
+        "The number of geometries after filling gaps should remain consistent."
+    )
+
+    # Assert the geometries differ when compact=True vs largest=True
+    assert not filled_gdf_compact.equals(filled_gdf_default), (
+        "The resulting geometries should differ between compact=True and largest=True."
+    )
+
+    # Verify that gaps are filled completely
+    filled_gaps_compact = filled_gdf_compact.unary_union
+    filled_gaps_default = filled_gdf_default.unary_union
+    assert filled_gaps_compact.is_valid, "The geometries with compact=True must be valid."
+    assert filled_gaps_default.is_valid, "The geometries with largest=True must be valid."
+

geoplanar/tests/test_gap.py

@@ -36,17 +36,11 @@ def test_fill_gaps(self):
         assert_equal(gdf1.area.values, numpy.array([108.0, 32.0]))
 
     def test_fill_gaps_smallest(self):
-        gdf1 = fill_gaps(self.gdf, largest=False)
-        assert_equal(gdf1.area.values, numpy.array([100.0, 40.0]))
-
-        gdf1 = fill_gaps(self.gdf_str, largest=False)
+        gdf1 = fill_gaps(self.gdf, strategy='smallest')
         assert_equal(gdf1.area.values, numpy.array([100.0, 40.0]))
 
     def test_fill_gaps_none(self):
-        gdf1 = fill_gaps(self.gdf, largest=None)
-        assert_equal(gdf1.area.values, numpy.array([108.0, 32.0]))
-
-        gdf1 = fill_gaps(self.gdf_str, largest=None)
+        gdf1 = fill_gaps(self.gdf, strategy=None)
         assert_equal(gdf1.area.values, numpy.array([108.0, 32.0]))
 
     def test_fill_gaps_gaps_df(self):

geoplanar/tests/test_overlap.py

@@ -4,6 +4,8 @@
 import numpy
 from numpy.testing import assert_equal
 from shapely.geometry import box
+from packaging.version import Version
+import pytest
 
 from geoplanar.overlap import (
     is_overlapping,
@@ -34,29 +36,30 @@ def test_trim_overlaps(self):
         assert_equal(gdf1.area.values, numpy.array([96.0, 8.0]))
 
     def test_trim_overlaps_smallest(self):
-        gdf1 = trim_overlaps(self.gdf, largest=False)
-        assert_equal(gdf1.area.values, numpy.array([100.0, 4.0]))
-
-        gdf1 = trim_overlaps(self.gdf_str, largest=False)
+        gdf1 = trim_overlaps(self.gdf, strategy='smallest')
         assert_equal(gdf1.area.values, numpy.array([100.0, 4.0]))
 
     def test_trim_overlaps_random(self):
-        gdf1 = trim_overlaps(self.gdf, largest=None)
-        assert_equal(gdf1.area.values, numpy.array([100.0, 4.0]))
-
-        gdf1 = trim_overlaps(self.gdf_str, largest=None)
+        gdf1 = trim_overlaps(self.gdf, strategy=None)
         assert_equal(gdf1.area.values, numpy.array([100.0, 4.0]))
 
+    @pytest.mark.skipif(Version(geopandas.__version__) == Version("0.10.2"), reason="Missing pygeos")    
     def test_trim_overlaps_multiple(self):
-        gdf1 = trim_overlaps(self.gdf2, largest=False)
-        assert_equal(gdf1.area.values, numpy.array([100.0, 100.0, 0.0]))
+        gdf1 = trim_overlaps(self.gdf2, strategy='largest')
+        assert_equal(gdf1.area.values, numpy.array([96, 96.0, 8.0]))
 
-        gdf1 = trim_overlaps(self.gdf2, largest=None)
+        gdf1 = trim_overlaps(self.gdf2, strategy=None)
         assert_equal(gdf1.area.values, numpy.array([100.0, 100.0, 0.0]))
 
         gdf1 = trim_overlaps(self.gdf2)
         assert_equal(gdf1.area.values, numpy.array([96.0, 96.0, 8.0]))
 
+        gdf1 = trim_overlaps(self.gdf2, strategy='smallest')
+        assert_equal(gdf1.area.values, numpy.array([100.0, 100.0, 0.0]))
+
+        gdf = trim_overlaps(self.gdf2, strategy='compact')
+        assert_equal(gdf1.area.values, numpy.array([100.0, 100.0, 0.0]))
+
     def test_merge_overlaps(self):
         gdf1 = merge_overlaps(self.gdf, 10, 0)
         assert_equal(gdf1.area.values, numpy.array([104]))