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


1 Answer 1


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

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.