Concave Hull with geopandas, scipy and numpy

concave_hull.md

Concave Hull in Python

A modified version for calculating concave hull (see source) based on GeoDataFrame of points

from scipy.spatial import Delaunay
import numpy as np
from shapely.ops import polygonize, cascaded_union
from shapely.geometry import MultiLineString
import pygeos
import geopandas as gpd

def concave_hull(points_gdf, alpha=35):
    """
    Compute the concave hull (alpha shape) of a GeoDataFrame of points.
    
    
    Parameters
    ==========
    points_gdf : gpd.GeoDataFrame
      GeoDataFrame of points.
    
    alpha: int
    
      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_gdf) < 4:
        # When you have a triangle, there is no sense
        # in computing an alpha shape.
        return points_gdf.unary_union.convex_hull

    coords = pygeos.coordinates.get_coordinates(points_gdf.geometry.values.data)
    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 = MultiLineString(edge_points)
    triangles = list(polygonize(m))
    return gpd.GeoDataFrame({"geometry": [cascaded_union(triangles)]}, index=[0], crs=points_gdf.crs)
```

Hi! I'm new with Geopandas, could you please give me an example of how "points_gdf" should be? I tried to use your function but it didn't work out. Basically, I have a Pandas dataframe with 3 columns: 'AGRUPAMENTO', 'Longitude' and 'Latitude'. I tried to do like that:

from shapely.geometry import Point

points = []
for xy in zip(Df['Longitude'], Df['Latitude']):
  points(Point(xy))

Gdf = gpd.GeoDataFrame(geometry=points)
Gdf_ConcaveHull = concave_hull(Gdf, alpha=100)

I have also tried:

Gdf = gpd.GeoDataFrame(geometry=gpd.points_from_xy(Df.Longitude, Df.Latitude))
Gdf_ConcaveHull = concave_hull(Gdf, alpha=100)

In both cases, I get:
TypeError: One of the arguments is of incorrect type. Please provide only Geometry objects.

Sorry if my English is bad.

@SaraSampaio You're right. In case you have GeoDataFrame, you need to convert that to Pygeos geometry first. I have slightly modified the script as follows:

def concave_hull(points_gdf, alpha=35):
    # Convert geopandas to pygeos Geometry object
    multipolygon = gdf.geometry.astype(str)[0]
    multipolygon = pygeos.Geometry(multipolygon)    
    coords = pygeos.coordinates.get_coordinates(multipolygon)
    
    # Create Delaunay voxels
    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 = MultiLineString(edge_points)
    triangles = list(polygonize(m))
    return gpd.GeoDataFrame({"geometry": [cascaded_union(triangles)]}, index=[0], crs='epsg:4326')

Reminder: The epsg must match with your EPSG.