Instantly share code, notes, and snippets.

# dwyerk/concave_hulls.ipynb

Created April 12, 2014 23:20
Show Gist options
• Save dwyerk/10561690 to your computer and use it in GitHub Desktop.
concave hulls using shapely and scipy
Display the source blob
Display the rendered blob
Raw
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.

### geodancer commented Apr 19, 2018

Please, may I have the import shapefile "concave_demo_points.shp" ? I want to play, mimic, learn.

Never mind. I found the data and recreated the examples. Thank you! On to the "learn" part. :-)

### caffeine-potent commented Aug 29, 2018

Nice Job. Going to use this as a reference.

### simoncozens commented Jan 25, 2019

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))

### vitale232 commented Jul 23, 2019

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:

### nys09 commented Feb 1, 2020

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.

### drdarina commented Apr 7, 2020

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.

### vic1309 commented Jul 29, 2020

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

### dwyerk commented Jul 30, 2020

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

### vic1309 commented Jul 30, 2020

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

### dwyerk commented Jul 30, 2020

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

### vic1309 commented Jul 31, 2020

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 :)

### dwyerk commented Jul 31, 2020

Looks like a great presentation!

### dd-debug commented Oct 24, 2020

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

### dwyerk commented Oct 27, 2020

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!

### DanielLauferPhysics commented Dec 29, 2020

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:

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

### potatoheadbig commented Aug 13, 2021

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.

### dwyerk commented Sep 1, 2021

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.

### hhopson-primrose commented Mar 24, 2023

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

### corsair20141 commented Nov 22, 2023

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

Is the `.vertices` method deprecated for Delaunay?

### dwyerk commented Nov 27, 2023

Yeah looks like it scipy/scipy@99cc995

Switch to tri.simplices

### corsair20141 commented Nov 28, 2023

Yeah looks like it scipy/scipy@99cc995

Switch to tri.simplices

That fixed it. Thanks!