{{ message }}

Instantly share code, notes, and snippets.

# dwyerk/concave_hulls.ipynb

Created Apr 12, 2014
concave hulls using shapely and scipy
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.

### geodancer commented Apr 28, 2018 • edited

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

### 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 I hope that helps address your question!

### 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.