concave hulls using shapely and scipy

concave_hulls.ipynb

Here's a vectorized version of alpha_shape which uses Numpy indexing to avoid the loop:

def alpha_shape(points, alpha):
    """
    Compute the alpha shape (concave hull) of a set
    of points.
    @param points: Iterable container of points.
    @param alpha: alpha value to influence the
        gooeyness of the border. Smaller numbers
        don't fall inward as much as larger numbers.
        Too large, and you lose everything!
    """
    if len(points) < 4:
        # When you have a triangle, there is no sense
        # in computing an alpha shape.
        return geometry.MultiPoint(list(points)).convex_hull

    coords = np.array([point.coords[0] for point in points])
    tri = Delaunay(coords)
    triangles = coords[tri.vertices]
    a = ((triangles[:,0,0] - triangles[:,1,0]) ** 2 + (triangles[:,0,1] - triangles[:,1,1]) ** 2) ** 0.5
    b = ((triangles[:,1,0] - triangles[:,2,0]) ** 2 + (triangles[:,1,1] - triangles[:,2,1]) ** 2) ** 0.5
    c = ((triangles[:,2,0] - triangles[:,0,0]) ** 2 + (triangles[:,2,1] - triangles[:,0,1]) ** 2) ** 0.5
    s = ( a + b + c ) / 2.0
    areas = (s*(s-a)*(s-b)*(s-c)) ** 0.5
    circums = a * b * c / (4.0 * areas)
    filtered = triangles[circums < (1.0 / alpha)]
    edge1 = filtered[:,(0,1)]
    edge2 = filtered[:,(1,2)]
    edge3 = filtered[:,(2,0)]
    edge_points = np.unique(np.concatenate((edge1,edge2,edge3)), axis = 0).tolist()
    m = geometry.MultiLineString(edge_points)
    triangles = list(polygonize(m))
    return cascaded_union(triangles), edge_points

Thanks to both of you! Both of these functions were very useful to me... The original is much easier to read, and the vectorized version is much faster.

From the example, I figured that alpha had to be <=1.0. That's definitely not the case. Here's an example where I set alpha=20 to get the desired result:

Can we use this if the points are in 3D (x,y,z)? I am trying to calculate the volume after fitting a concave hull to the points.

Thanks, this was really helpful! I would suggest adding a parameter for a default distance between coordinates, e.g. filtered = triangles[circums < default_distance / alpha. In my case the points were on average 2000 units apart and it was easier to adjust this parameter than to rescale the data.

Hello! Thank you for the great tutorial.

I would like to know if any of you could estimate the total area of the hull. Any suggestions?

Cheers

Since the polygon itself has no units, you will have to first project the polygon using shapely, then take the area.

This looks like a good example: https://gis.stackexchange.com/a/128072

Dear dwyerk

Following your suggestion, I did the following:

Obtained the (lat, lon) hull values using from shapely.geometry import LineString and then, with the boundary values in hand, I projected them to the Earths surface using Pyproj and finally estimated the area using from shapely.geometry import shape. I can provide a code snippet if any of you want it.

Cheers

I'm glad it worked for you! Definitely keep up the good karma and share an example.

Hello dwyerk

I prepared this example (https://github.com/vic1309/concave_hull_area) with a test data and a Jupyter notebook. Please, feel free to make suggestions and play around with them. A feedback would be just awesome :)

Looks like a great presentation!

Hello dwyerk, I think your work is perfect. Can it be extended to higher dimension (>3)?

Hello dwyerk, I think your work is perfect. Can it be extended to higher dimension (>3)?

Hi dd-debug. There is probably an algorithm for concave hulls in 3D but I think it's more complicated than adjusting this one. Since the triangulation is 2D by nature, you'll need to find a way to work with surfaces. Good luck!

Hey dd-debug, there's already an algorithm for computing N-dimensional concave hulls. You can find discussion of the algorithm and intuition for it in the paper A New Concave Hull Algorithm and Concaveness Measure for n-dimensional Datasets by Jin-Seo Park and Se-Jong Oh, the link for which I've attached here:

https://journal.iis.sinica.edu.tw/paper/1/100295-3.pdf?cd=2217EEBB7C44EDA26.

I've tinkered with an implementation of it that's available on github which I've trial-ran in the past couple of days and I'm getting pleasing results, though I have by no means looked through it in detail for optimization tweaks. You can find that here:

https://gist.github.com/AndreLester/589ea1eddd3a28d00f3d7e47bd9f28fb

I hope that helps address your question!

What if I don't want to scale up and add more points? That is, given the H shape at the first beginning, how to get the concave hull? I've tried multiple alpha shape libraries but none of them worked for a set of points which only include the boundary without any interior points.

You can still try, but you may not be able to find an acceptable alpha value. It's been a few years since I was working on this but if I recall correctly, the delaunay triangulation doesn't generate enough triangles to be later pruned, or at least the space between them is not great enough.

# Area of triangle by Heron's formula
area = math.sqrt(s*(s-a)*(s-b)*(s-c))
circum_r = a*b*c/(4.0*area)
if circum_r < 1.0/alpha:
    # Keep these edges

Those areas need to be low enough to be less than the inverse of the alpha parameter. I'd step through this in a debugger with your shapefile and see what kind of values you're getting for each triangle to see if you can find a better set of parameters.

FYI for anyone running this notebook today (unless Descartes receives an update), you need to patch Descartes as described in this stackoverflow answer or the plot_polygon won't work https://stackoverflow.com/questions/75287534/indexerror-descartes-polygonpatch-wtih-shapely

I get the error "AttributeError: 'Delaunay' object has no attribute 'vertices'"

Is the .vertices method deprecated for Delaunay?

Yeah looks like it scipy/scipy@99cc995

Switch to tri.simplices

Yeah looks like it scipy/scipy@99cc995

Switch to tri.simplices

That fixed it. Thanks!

I am trying to get this to work for two offset sine waves. However, the best I can get is the second image. After that, increasing the alpha parameter creates a segmented geometry. Perhaps my ordering of points isn't in the proper format for this function?