Aggregating Close Polygons with GDAL?

Question

I would like to aggregate polygons that are within 2 meters of each other, like the Aggregate Polygons (Cartography) tool in ArcMap, but using GDAL with Python.

To add some detail in case the problem can be approached earlier in the pipeline, these polygons were created using Polygonize in GDAL from PNG binary masks with the following code:

dst_layername = blockname + "_Piles"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_layername + ".shp" )
dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )

gdal.Polygonize( srcband, None, dst_layer, -1, [], callback=None )

To demonstrate what I mean using the ArcMap tool, here are my polygons before:

And here are my polygons after running the Aggregate Polygons tool:

I wish to do this using non-ESRI tools that I can call from Python. I have been using osgeo and GDAL for my project thus far.

Answer 1

As the solution suggested by Jon, simply use the solutions proposed by Removing small polygons gaps in a Shapely polygon or Python Shapely - filling small gaps between multi polygons with GeoPandas (based on ogr) to read/write the shapefiles and shapely.

1) Aggregate close polygons

import geopandas as gpd
from shapely.geometry import MultiPolygon, JOIN_STYLE
import itertools
eps=5 # width for dilating and eroding (buffer)
dist = 2  # threshold distance
# read the original shapefile
df = gpd.read_file('polygons.shp')
# create new result shapefile
col = ['geometry']
res = gpd.GeoDataFrame(columns=col)
# iterate over pairs of polygons in the GeoDataFrame 
for i, j in list(itertools.combinations(df.index, 2)):
 distance = df.geometry.ix[i].distance( df.geometry.ix[j]) # distance between polygons i and j in the shapefile
 if distance < dist: 
     e = MultiPolygon([ df.geometry.ix[i],df.geometry.ix[j]])
     fx = e.buffer(eps, 1, join_style=JOIN_STYLE.mitre).buffer(-eps, 1, join_style=JOIN_STYLE.mitre)
     res = res.append({'geometry':fx},ignore_index=True)

# save the resulting shapefile   
res.to_file("aggregates.shp")

2) it's easy to complete the script to include the remaining polygons.

Answer 2

I don't know if there are packages that do what you want in a single command, but I would probably just buffer the polygon layer by whatever your threshold distance is, union any intersecting polygons, then perform a negative buffer of the same size. This might distort the polygons ever-so-slightly. The geopandas library has all the required tools as methods: read_file(), buffer(), join(), and to_file() methods should be all you need. You can also use gdal commands to do the same thing, but I am not as familiar with those. See here for buffering a shapefile with gdal.

Answer 3

The answer by Gene is effective and works well with a working GeoPandas installation. However, because the original polygons were created using Polygonize from binary image masks, the problem can be alternatively approached prior to the polygonize step.

Specifically, the use of OpenCV's Dilation and Erosion operations. This is similar in principle to Gene and Jon's answers but acts on the raster rather than resultant vectors. The cv.dilate operation expands the shapes and in doing so merges them, then the cv.erode reduces the shapes to their original size.

Sample code:

import cv2 as cv
import numpy as np

input = 'original_image.png'
output = 'resultant_image.png'

splitpiles = cv.imread(input,0)

#Set the dilation threshold at 20,20 and specify data type
kernel = np.ones((20,20),np.uint8)

#Run the dilation then the erode operations
dilation = cv.dilate(splitpiles,kernel,iterations = 1)
mergedpiles = cv.erode(dilation,kernel,iterations = 1)

# Convert image from Grey to RGB so Polygonize will like it, and save
mergedpiles = cv.cvtColor(mergedpiles,cv.COLOR_GRAY2RGB)
cv.imwrite(output,mergedpiles)