The coordinates of the points in poly1 and poly2 are not an exact match (despite the fact that they look like the same point on the plot).
I guess you want the algorithm to work with "close enough" values (i.e. if the distance between the points is less than a given tolerance, e.g. 0.1, they are considered the same point).
You can join the two polygons by finding a common edge and removing it.
To find the common edge, we first determine which points of the polygons are common to both polygons.
I had a look at this post and chose the cKDTree method to do that.
It works by
- finding for each point of one polygon its closest neighbouring point in the other polygon
- comparing the distance between those two points. If the distance is less than our set tolerance, we consider them the same point and common to both polygons.
Once you've determined, which points are common, you can verify whether they are adjacent. If they are, you've found the edge to remove.
The resulting polygon will be composed of
- all points from
poly1
- those points from
poly2 which do not form the common edge
Here is the code:
import matplotlib.pyplot as plt, numpy as np
from mpl_toolkits.mplot3d import Axes3D
from scipy.spatial import cKDTree
poly1 = np.array([[120787.075999871, 491779.675000143, -2.0699999332428], [120784.319999829, 491781.831000042, 5.96999979019165], [120784.319999829, 491781.831000042, -2.0699999332428], [120787.075999871, 491779.675000143, -2.0699999332428]])
poly2 = np.array([[120787.075999871, 491779.675000143, -2.03999996185303], [120787.075999871, 491779.675000143, 5.90999984741211], [120784.319999829, 491781.831000042, 5.96999979019165], [120787.075999871, 491779.675000143, -2.03999996185303]])
def is_close(a, b, tolerance):
# Get closest distances for each pt in a
dist = cKDTree(b).query(a, k=1)[0] # k=1 selects closest one neighbor
# Check the distances against the given tolerance value
return dist <= tolerance
def find_consecutive_true_values(arr):
i = 0
while i < len(arr) - 1:
if arr[i] and arr[i+1]:
return i
i+=1
raise Exception('No common edge found')
# Find points in poly1, which are closer than given tolerance to any point in poly2
# and vice versa
tolerance = 0.1
points_in_poly1_close_to_poly2 = is_close(poly1, poly2, tolerance)
points_in_poly2_close_to_poly1 = is_close(poly2, poly1, tolerance)
# Scan each array for two adjacent true values (points at those two indices
# form an edge which is common to both polygons and which we want to remove).
# Idx1 (resp. idx2) will contain the index of the first point of that common edge in poly1 (resp. poly2)
idx1 = find_consecutive_true_values(points_in_poly1_close_to_poly2)
idx2 = find_consecutive_true_values(points_in_poly2_close_to_poly1)
#Split poly1 into two parts:
# first part contains points from the start up to the first point of the common edge (inclusive)
# second part contains points from the second point of the common edge to the end
poly1_part1 = poly1[:idx1+1]
poly1_part2 = poly1[idx1+1:]
#Remove common edge from poly2, depending on where it is located, we end up with one or two parts
if idx2 == len(poly2) - 2:
poly2_part1 = poly2[1:len(poly2) - 2]
poly2_part2 = None
elif idx2 == 0:
poly2_part1 = poly2[2:len(poly2) - 1]
poly2_part2 = None
else:
poly2_part1 = poly2[idx2+2:]
poly2_part2 = poly2[1:idx2]
#Create the resulting polygon by concatenating the individual parts (poly2_part2 may be empty)
if(poly2_part2 is None):
poly = np.concatenate((poly1_part1, poly2_part1, poly1_part2))
else:
poly = np.concatenate((poly1_part1, poly2_part1, poly2_part2, poly1_part2))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(poly[:,0], poly[:,1], poly[:,2])
ax.view_init(45, 45)
plt.show()
(the code is far from idiomatic, if you know your Python, feel free to edit it :) )
