Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
concave hulls using shapely and scipy
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
@geodancer

This comment has been minimized.

Copy link

@geodancer geodancer commented Apr 19, 2018

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

@geodancer

This comment has been minimized.

Copy link

@geodancer geodancer commented Apr 28, 2018

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

@caffeine-potent

This comment has been minimized.

Copy link

@caffeine-potent caffeine-potent commented Aug 29, 2018

Nice Job. Going to use this as a reference.

@simoncozens

This comment has been minimized.

Copy link

@simoncozens 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))
    return cascaded_union(triangles), edge_points
@vitale232

This comment has been minimized.

Copy link

@vitale232 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:
image

@nys09

This comment has been minimized.

Copy link

@nys09 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

This comment has been minimized.

Copy link

@drdarina 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

This comment has been minimized.

Copy link

@vic1309 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

This comment has been minimized.

Copy link
Owner Author

@dwyerk 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

This comment has been minimized.

Copy link

@vic1309 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

This comment has been minimized.

Copy link
Owner Author

@dwyerk dwyerk commented Jul 30, 2020

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

@vic1309

This comment has been minimized.

Copy link

@vic1309 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

This comment has been minimized.

Copy link
Owner Author

@dwyerk dwyerk commented Jul 31, 2020

Looks like a great presentation!

@dd-debug

This comment has been minimized.

Copy link

@dd-debug dd-debug commented Oct 24, 2020

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

@dwyerk

This comment has been minimized.

Copy link
Owner Author

@dwyerk 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

This comment has been minimized.

Copy link

@DanielLauferPhysics 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:

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!

@potatoheadbig

This comment has been minimized.

Copy link

@potatoheadbig 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

This comment has been minimized.

Copy link
Owner Author

@dwyerk 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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment