diff --git a/ci/37.yaml b/ci/37.yaml index a5a418604..8a9f43db5 100644 --- a/ci/37.yaml +++ b/ci/37.yaml @@ -15,6 +15,7 @@ dependencies: - pytest - pytest-cov - pytest-xdist + - networkx=2.6 # optional - geopandas>=0.7.0 - joblib diff --git a/libpysal/cg/alpha_shapes.py b/libpysal/cg/alpha_shapes.py index 2ae419009..27baa5960 100644 --- a/libpysal/cg/alpha_shapes.py +++ b/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 diff --git a/libpysal/cg/tests/test_ashapes.py b/libpysal/cg/tests/test_ashapes.py index a55f099c9..f9e53a216 100644 --- a/libpysal/cg/tests/test_ashapes.py +++ b/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)) \ No newline at end of file