GeoPandas Enveloping Polygon – Creating Minimal Enveloping Polygon to Set of Geometries in GeoPandas

concave hullgeopandaspolygonpolygon-creationpython

To be able to do operations on a set of geometries in a GeoPandas GeoDataFrame, I need to be able to determine whether objects are on the outer "rim" of the set. The set of geometries is as follows:

enter image description here

To do this, I would like to create a polygon that perfectly matches the outer bound of the set of geometrical objects. I first thought about using the convex hull of the set:

convex_hull = Sectioned_geostore_obstacles_geometry.unary_union.convex_hull
convex_hull = geopandas.GeoDataFrame({'geometry': convex_hull, 'convex_hull':[1]})

ax = Sectioned_geostore_obstacles_geometry['Gondola'].plot(color='red')
convex_hull.plot(ax=ax, color='green', alpha=0.5)

which results in

enter image description here

but this isn't quite right since what I am looking for isn't convex. The second idea is to use the envelope:

envelope = Sectioned_geostore_obstacles_geometry.unary_union.envelope
envelope = geopandas.GeoDataFrame({'geometry': envelope, 'convex_hull':[1]})

ax = Sectioned_geostore_obstacles_geometry['Gondola'].plot(color='red')
envelope.plot(ax=ax, color='green', alpha=0.5)

which is

enter image description here

Again, this isn't it. Yet another attempt is to use the cascade_union functionality from shapely:

from shapely.ops import cascaded_union

polygons = list(Sectioned_geostore_obstacles_geometry.Gondola)
boundary = gpd.GeoSeries(cascaded_union(polygons))

which is:

enter image description here

But, this isn't it either as it returns a MultiPolygon instead of the minimal developing polygon. Basically, I need the envelope to shrink to follow the contour of the set of objects.

To test this, I add the following example data:

test_df =  geopandas.GeoSeries([Polygon([(0,0), (2,0), (2,2), (0,2)]),
                              Polygon([(2,2), (4,2), (4,4), (2,4)])])
test_df = geopandas.GeoDataFrame({'geometry': test_df, 'df1':[1,2]})

convex_hull = test_df.unary_union.convex_hull
convex_hull = geopandas.GeoDataFrame({'geometry': convex_hull, 'convex_hull':[1]})

ax1 = test_df['geometry'].plot(color='red')
convex_hull.plot(ax=ax1, color='green', alpha=0.5)

envelope = test_df.unary_union.envelope
envelope = geopandas.GeoDataFrame({'geometry': envelope, 'convex_hull':[1]})

ax2 = test_df['geometry'].plot(color='red')
envelope.plot(ax=ax2, color='green', alpha=0.5)

enter image description here

enter image description here

Best Answer

What you need is a concave hull. Create a list of all polygons coordinates and concave hull them. This takes about 30 s for two polygon groups so try it on a subset if you have a very large dataset.

import geopandas as gpd
import alphashape #pip install alphashape
import re

df = gpd.read_file(r'/home/bera/Desktop/GIStest/buildings_two_groups.shp')

def giveVertices(frame):
    """A function to list all vertices in a polygon geometry"""
    vertices = [float(coord) for coord in re.findall('\d+\.\d+', frame.geometry.wkt)]
    return vertices
    
df['vertices'] = df.apply(giveVertices, axis=1) #Each polygons vertices as a list [418957.7407420929, ..., 418968.9631302697]
df2 = df.groupby('group')['vertices'].apply(list).reset_index() #The vertices of all polygons in each group, as a list of lists

def unnest(frame):
    """A function to turn a list of lists, into a list of tuples [(x coordinate, y coordinate), ...]"""
    flatList = [item for sublist in frame['vertices'] for item in sublist]
    listOfTuples = [(x,y) for x,y in zip(flatList[::2], flatList[1::2])]
    return listOfTuples

df2['vertices'] = df2.apply(unnest, axis=1) #One list of all polygons vertices in each group

def createConcavehull(frame):
    """Create a concave hull of each groups vertices"""
    alpha = alphashape.optimizealpha(frame['vertices'])/2 #I divide by two because a higher alpha made the hull
    #   look more like a star and not follow the building outlines so nicely.
    hull = alphashape.alphashape(frame['vertices'], alpha)
    return hull

df2['hull'] = df2.apply(createConcavehull, axis=1)
result = gpd.GeoDataFrame(df2, geometry='hull', crs="EPSG:3006")
result = result[['group','hull']] #Drop the vertices column, which is a list that shapefile outputs cant handle
result.to_file(r'/home/bera/Desktop/GIStest/buildings_con_hull_two_groups_2.shp')

enter image description here