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:

Before Polygon Aggregation

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

After Polygon Aggregation

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.


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.

enter image description here

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   

enter image description here

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

enter image description here

  • That's terrific, and I can get it to work (once I've ogr2ogr cleaned the shapefile from polygonize such that Fiona can read it without error). I think I understand what you're doing here, but how would I go about including the remaining polygons? – Joska May 24 '18 at 19:29

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.

  • That's a clever work-around. The slight distortion is likely insignificant. I will give it a try and report back. – Joska May 17 '18 at 19:32

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)
  • This may not actually be more efficient as the morphological operation, while fast, is applied to the full image. I've found that working in vector space can be much faster for similar applications. Also note that eroding and dilating can remove "fuzz" on the blob boundaries (which may be unimportant, but good to be aware of). – Jon Aug 22 '18 at 15:32
  • @Jon is correct; while in this specific case the original data was in raster form, the vector operations with geopandas are actually quite computationally efficient. – Joska Aug 22 '18 at 22:29

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.