To give you an idea of what would be involved to code it... consider the following
import numpy as np
def pnt_on_poly(pnt, poly):
"""Find closest point location on a polygon/polyline.
Parameters
----------
pnt : 1D ndarray array
XY pair representing the point coordinates.
poly : 2D ndarray array
A sequence of XY pairs in clockwise order is expected. The first and
last points may or may not be duplicates, signifying sequence closure.
Returns
-------
A list of [x, y, distance] for the intersection point on the line
Notes
-----
e_dist is represented by _e_2d and pnt_on_seg by its equivalent below.
This may be as simple as finding the closest point on the edge, but if
needed, an orthogonal projection onto a polygon/line edge will be done.
This situation arises when the distance to two sequential points is the
same
"""
def _e_2d_(a, p):
""" array points to point distance... mini e_dist"""
diff = a - p[np.newaxis, :]
return np.sqrt(np.einsum('ij,ij->i', diff, diff))
def _pnt_on_seg_(seg, pnt):
"""mini pnt_on_seg function normally required by pnt_on_poly"""
x0, y0, x1, y1, dx, dy = *pnt, *seg[0], *(seg[1] - seg[0])
dist_ = dx*dx + dy*dy
u = ((x0 - x1)*dx + (y0 - y1)*dy)/dist_
u = max(min(u, 1), 0)
xy = np.array([dx, dy])*u + [x1, y1]
return xy
def _line_dir_(orig, dest):
"""mini line direction function"""
orig = np.atleast_2d(orig)
dest = np.atleast_2d(dest)
dxy = dest - orig
ang = np.degrees(np.arctan2(dxy[:, 1], dxy[:, 0]))
return ang
pnt = np.asarray(pnt)
poly = np.asarray(poly)
if np.all(poly[0] == poly[-1]):
poly = poly[:-1]
d = _e_2d_(poly, pnt)
key = np.argsort(d)[0]
if key == 0:
seg = np.vstack((poly[-1:], poly[:3]))
elif (key + 1) >= len(poly):
seg = np.vstack((poly[-2:], poly[:1]))
else:
seg = poly[key-1:key+2]
n1 = _pnt_on_seg_(seg[:-1], pnt)
d1 = np.linalg.norm(n1 - pnt)
n2 = _pnt_on_seg_(seg[1:], pnt)
d2 = np.linalg.norm(n2 - pnt)
if d1 <= d2:
dest = [n1[0], n1[1]]
ang = _line_dir_(pnt, dest)
return [n1[0], n1[1], np.asscalar(d1), np.asscalar(ang)]
else:
dest = [n2[0], n2[1]]
ang = _line_dir_(pnt, dest)
return [n2[0], n2[1], np.asscalar(d2), np.asscalar(ang)]
a = np.array([[0, 0], [0, 100], [100, 100], [100, 80], [20, 80],
[20, 20], [100, 20], [100, 0], [0, 0]])
pnt = [150, 150]
pnt_on_poly(pnt, a)
but you would have to do a general sort of the closest point in one polygon to the extent of the other, then limit your candidate points before doing the calculation. One gains a greater appreciation of what goes on behinds the scenes for such a 'simple' question