[WIP] Filter holes from alpha shape returns

ci/37.yaml

@@ -15,6 +15,7 @@ dependencies:
   - pytest
   - pytest-cov
   - pytest-xdist
+  - networkx=2.6
   # optional
   - geopandas>=0.7.0
   - joblib

libpysal/cg/alpha_shapes.py

@@ -6,10 +6,12 @@
 
 Author(s):
     Dani Arribas-Bel daniel.arribas.bel@gmail.com
+    Levi John Wolf levi.john.wolf@gmail.com
 """
 
 import numpy as np
 import scipy.spatial as spat
+from scipy import sparse
 
 from ..common import requires, jit, HAS_JIT
 
@@ -354,7 +356,7 @@ def get_single_faces(triangles_is):
 
 
 @requires("geopandas", "shapely")
-def alpha_geoms(alpha, triangles, radii, xys):
+def _alpha_geoms(alpha, triangles, radii, xys):
     """Generate alpha-shape polygon(s) from `alpha` value, vertices of
     `triangles`, the `radii` for all points, and the points themselves
 
@@ -378,9 +380,9 @@ def alpha_geoms(alpha, triangles, radii, xys):
     -------
 
     geoms : GeoSeries
-        Polygon(s) resulting from the alpha shape algorithm. The GeoSeries
-        object remains so even if only a single polygon is returned. There is
-        no CRS included in the object.
+        Polygon(s) resulting from the alpha shape algorithm, in a GeoSeries.
+        The output is a GeoSeries even if only a single polygon is returned. There is
+        no CRS included in the returned GeoSeries.
 
     Examples
     --------
@@ -443,7 +445,8 @@ def alpha_shape(xys, alpha):
     shapes : GeoSeries
          Polygon(s) resulting from the alpha shape algorithm. The GeoSeries
          object remains so even if only a single polygon is returned. There is
-         no CRS included in the object.
+         no CRS included in the object. Note that the returned shape(s) may
+         have holes, as per the definition of the shape in Edselbrunner et al. (1983)
 
     Examples
     --------
@@ -480,7 +483,8 @@ def alpha_shape(xys, alpha):
     c_pts = triangles[:, 2, :]
     radii = r_circumcircle_triangle(a_pts, b_pts, c_pts)
     del triangles, a_pts, b_pts, c_pts
-    geoms = alpha_geoms(alpha, triangulation.simplices, radii, xys)
+    geoms = _alpha_geoms(alpha, triangulation.simplices, radii, xys)
+    geoms = _filter_holes(geoms, xys)
     return geoms
 
 
@@ -526,12 +530,12 @@ def alpha_shape_auto(
 
     This method uses the algorithm proposed by  Edelsbrunner, Kirkpatrick &
     Seidel (1983) to return the tightest polygon that contains all points in
-    `xys`. The algorithm ranks every point based on its radious and iterates
+    `xys`. The algorithm ranks every point based on its radius and iterates
     over each point, checking whether the maximum alpha that would keep the
     point and all the other ones in the set with smaller radii results in a
     single polygon. If that is the case, it moves to the next point;
     otherwise, it retains the previous alpha value and returns the polygon
-    as `shapely` geometry.
+    as `shapely` geometry. Note that this geometry may have holes.
 
     Parameters
     ----------
@@ -616,7 +620,7 @@ def alpha_shape_auto(
     radii_sorted_i = radii.argsort()
     triangles = triangulation.simplices[radii_sorted_i][::-1]
     radii = radii[radii_sorted_i][::-1]
-    geoms_prev = alpha_geoms((1 / radii.max()) - EPS, triangles, radii, xys)
+    geoms_prev = _alpha_geoms((1 / radii.max()) - EPS, triangles, radii, xys)
     if HAS_PYGEOS:
         points = pygeos.points(xys)
     else:
@@ -628,7 +632,7 @@ def alpha_shape_auto(
         alpha = (1 / radi) - EPS
         if verbose:
             print("%.2f%% | Trying a = %f" % ((i + 1) / radii.shape[0], alpha))
-        geoms = alpha_geoms(alpha, triangles, radii, xys)
+        geoms = _alpha_geoms(alpha, triangles, radii, xys)
         if _valid_hull(geoms, points):
             geoms_prev = geoms
             radi_prev = radi
@@ -704,6 +708,36 @@ def _construct_centers(a, b, radius):
         return down_x, down_y
 
 
+def _filter_holes(geoms, points):
+    """
+    Filter hole polygons using a computational geometry solution
+    """
+    if (geoms.interiors.apply(len) > 0).any():
+        from shapely.geometry import Polygon
+
+        # Extract the "shell", or outer ring of the polygon.
+        shells = geoms.exterior.apply(Polygon)
+        # Compute which original geometries are within each shell, self-inclusive
+        inside, outside = shells.sindex.query_bulk(geoms, predicate="within")
+        # Now, create the sparse matrix relating the inner geom (rows)
+        # to the outer shell (cols) and take the sum.
+        # A z-order of 1 means the polygon is only inside if its own exterior. This means it's not a hole.
+        # A z-order of 2 means the polygon is inside of exactly one other exterior. Because
+        #   the hull generation method is restricted to be planar, this means the polygon is a hole.
+        # In general, an even z-order means that the polygon is always exactly matched to one exterior,
+        #   plus some number of intermediate exterior-hole pairs. Therefore, the polygon is a hole.
+        # In general, an odd z-order means that there is an uneven number of exteriors.
+        #   This means the polygon is not a hole.
+        zorder = sparse.csc_matrix((np.ones_like(inside), (inside, outside))).sum(
+            axis=1
+        )
+        zorder = np.asarray(zorder).flatten()
+        # Keep only the odd z-orders
+        to_include = (zorder % 2).astype(bool)
+        geoms = geoms[to_include]
+    return geoms
+
+
 if __name__ == "__main__":
 
     import matplotlib.pyplot as plt

libpysal/cg/tests/test_ashapes.py

@@ -98,3 +98,20 @@ def test_circles(self):
         ashape, radius, centers = alpha_shape_auto(self.vertices, return_circles=True)
         np.testing.assert_allclose(radius, self.circle_radii)
         np.testing.assert_allclose(centers, self.circle_verts)
+
+    def test_holes(self):
+        np.random.seed(seed=100)
+        points = np.random.rand(1000, 2)*100
+        inv_alpha = 3.5
+        geoms = alpha_shape(points, 1/inv_alpha)
+        assert len(geoms) == 1
+        holes = geopandas.GeoSeries(geoms.interiors.explode()).reset_index(drop=True)
+        assert len(holes) == 30
+        # No holes are within the shape (shape has holes already)
+        result = geoms.sindex.query_bulk(holes.centroid, predicate='within')
+        assert result.shape == (2,0)
+        # All holes are within the exterior
+        shell = geopandas.GeoSeries(geoms.exterior.apply(geometry.Polygon))
+        within, outside = shell.sindex.query_bulk(holes.centroid, predicate='within')
+        assert (outside == 0).all()
+        np.testing.assert_array_equal(within, np.arange(30))