How to merge almost touching polygons

Question

How can I merge polygons that are almost touching, using Python? For example, given these polygons and some distance threshold T:

The algorithm would produce something like this:

So basically I need to remove any gaps smaller than T between the two polygons.

The original vertices of the polygons should not be unnecessarily modified. So for example simply buffering both polygons and then taking a union of them is not good.

A library that already implements this functionality would also be fine.

EDIT: As another example, ArcGIS polygon aggregation seems to do exactly this. However there isn't much hint about how it is implemented, and using ArcGIS is not an option.

You could define a 'merging' distance, and merge each two vertices which are less than the merging distance away. — Willem Hendriks, Mar 13 '23 at 12:04
Unfortunately that won't work. Take polygons P1 and P2. Vertex V of P1 can be infinitely close to P2 if it sits right next to the center point of any of P2's **edges**, while still not being within the merging distance of your algorithm. — anttikoo, Mar 14 '23 at 08:12
Would also take edge distance into consideration solve this? — Willem Hendriks, Mar 14 '23 at 08:14
Maybe use a morphological dilation https://stackoverflow.com/a/74997349/2836621 — Mark Setchell, Mar 14 '23 at 08:58
@WillemHendriks I suppose it takes care of that particular problem. However many questions remain. Once you have merged all such vertices (I guess merging vertices means adding an edge between them) and edges (split the edge, insert vertice and then merge?), you have a mishmash of edges between the two polygons. Some of which should be removed, some should represent the exterior of the end result, and some might even represent an interior (a hole that was introduced) in the end result. Not sure where to go from there. — anttikoo, Mar 14 '23 at 08:58
@MarkSetchell seems like the raster equivalent of buffering the vector. Which is not ok, as it impacts all of the vertices of the input polygons, not only those in the "gap" between. — anttikoo, Mar 14 '23 at 09:01
Why is the bridging between the two shapes diagonal at (6,4) instead of vertical - if the existing vertices should remain unchanged? — Mark Setchell, Mar 14 '23 at 09:24
Maybe you can use **OpenCV** *"Contour Approximation"* https://docs.opencv.org/4.x/dd/d49/tutorial_py_contour_features.html — Mark Setchell, Mar 14 '23 at 10:17
@MarkSetchell the latter image shows just an example of what the algorithm might produce. i stated that vertices should not be _unnecessarily_ modified. For example buffering the whole polygon yields unnecessary modifications, because that changes also those vertices that are nowhere near the gap to be bridged. As far as I know, OpenCV only deals with raster data (images), so it has no use here. — anttikoo, Mar 14 '23 at 11:29

score 0 · Answer 1 · answered Mar 13 '23 at 17:06

Here is what I have tried:

I have used your threshold T to draw a circle, where the center is base on the nearest point between the points in polygon_1 ad all points in polygon_2. Then find the intersections, and find the convex hull for the intersection zone. Then find the union of the three polygons, codes :

import numpy as np
import matplotlib.pyplot as plt
from shapely.ops import unary_union
from shapely.geometry import Polygon
from scipy.spatial import ConvexHull


polygon_1 = np.array([[1, 6, 7.8, 7.8, 1, 1], [4, 4, 6, 11, 11, 4]]).T
polygon_2 = np.array([[6, 14, 14, 11, 8.2, 9, 6, 6], [0.5, 0.5, 12, 12, 10, 4.2, 3.5, 0.5]]).T

poly_1 = Polygon(polygon_1)
poly_2 = Polygon(polygon_2)

# define threshold T
T = 2.5

#Initialize a variable to create a circle
theta = np.linspace(0, 2*np.pi, 100)



intersection_points  = []
coord_circles = np.zeros((0, 2))
for poly in polygon_1:
    # Calculate the distance between the points in polygon_1 ad all points in polygon_2
    norm = np.sqrt((poly[0] - polygon_2[:, 0])**2 + (poly[1] - polygon_2[:, 1])**2)
    
    # If the distance is smaaller than the threshold
    if np.min(norm) < T:
        # Find the nearest point in polygon_2
        point_near =polygon_2[ np.argmin(norm)]
        x = T*np.cos(theta) + point_near[0]
        y = T*np.sin(theta) + point_near[1]
        
        # Crear a polygon for the circle and find intersection points with polygon_1 and polygon_2
        poly_circle = Polygon(np.array([x, y]).T)
        intersec_1 = poly_1.boundary.intersection(poly_circle.boundary)
        intersection_points.append(intersec_1.geoms[0].coords._coords[0])
        intersection_points.append(intersec_1.geoms[1].coords._coords[0])
        intersec_2 = poly_2.boundary.intersection(poly_circle.boundary)
        intersection_points.append(intersec_2.geoms[0].coords._coords[0])
        intersection_points.append(intersec_2.geoms[1].coords._coords[0])
        
intersection_points = np.array(intersection_points)   

# Find the Convex hulls 
hull = ConvexHull(intersection_points)

# Creat a polygon for the convex hull of the intersections 
poly_3 = Polygon(np.c_[intersection_points[hull.vertices,0], intersection_points[hull.vertices,1]])


# Find the union of the these three polygon
a = unary_union([poly_1, poly_2, poly_3])

result_x, result_y = a.exterior.coords.xy
result_x = result_x.tolist()
result_y = result_y.tolist()

plt.subplot(121) 
plt.plot(polygon_1[:, 0], polygon_1[:, 1], 'b')
plt.plot(polygon_2[:, 0], polygon_2[:, 1], 'r')
plt.axis('equal')
plt.title("Before merging")
plt.subplot(122) 
plt.plot(result_x, result_y, 'k-o')
plt.axis('equal')
plt.title("After merging")

The you get figure :

Thank you for taking the time to come up with this. But the convex part breaks down if it is given a suitable concave polygon. For example `poly_1 = LineString([(1, 1), (3, 3), (1, 5)]).buffer(0.2)` and `poly_2 = LineString([(2, 1), (4, 3), (2, 5)]).buffer(0.2)`. — anttikoo, Mar 14 '23 at 08:30

How to merge almost touching polygons

1 Answers1