5

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

Taras
  • 32,823
  • 4
  • 66
  • 137

2 Answers2

6

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

BERA
  • 72,339
  • 13
  • 72
  • 161
1

Nowadays you can get the convex hull with GeoPandas very easily:

convex_hull = buildings[buildings.is_valid].unary_union.convex_hull 
convex_hull = gpd.GeoDataFrame(geometry=[convex_hull], crs=buildings.crs)

enter image description here

Phil
  • 161
  • 6
  • it doesn't seem to be working with the latest geopandas version. I've installed the latest version of geopandas but this function is not included – foxhq Aug 31 '23 at 09:30