查找点是否位于多边形内部或边界上（Python - Shapely）

Question

好的，我需要检查一个点是否在多边形内部或位于其边界上，但是对于这种情况我无法使其正常工作。 首先是一些函数和变量：

from matplotlib import pyplot as plt    
from shapely.geometry import Polygon,Point,LinearRing,LineString
from math import sqrt
from shapely import affinity
from descartes.patch import PolygonPatch

GM = (sqrt(5)-1.0)/2.0
W = 8.0
H = W*GM
SIZE = (W, H)

BLUE = '#6699cc'
GRAY = '#999999'
DARKGRAY = '#333333'
YELLOW = '#ffcc33'
GREEN = '#339933'
RED = '#ff3333'
BLACK = '#000000'

COLOR_ISVALID = {
    True: BLUE,
    False: RED,
}

def plot_line(ax, ob, color=GRAY, zorder=1, linewidth=3, alpha=1):
    x, y = ob.xy
    ax.plot(x, y, color=color, linewidth=linewidth, solid_capstyle='round', zorder=zorder, alpha=alpha)

def plot_coords(ax, ob, color=GRAY, zorder=1, alpha=1):
    x, y = ob.xy
    ax.plot(x, y, 'o', color=color, zorder=zorder, alpha=alpha)

def color_isvalid(ob, valid=BLUE, invalid=RED):
    if ob.is_valid:
        return valid
    else:
        return invalid

def color_issimple(ob, simple=BLUE, complex=YELLOW):
    if ob.is_simple:
        return simple
    else:
        return complex

def plot_line_isvalid(ax, ob, **kwargs):
    kwargs["color"] = color_isvalid(ob)
    plot_line(ax, ob, **kwargs)

def plot_line_issimple(ax, ob, **kwargs):
    kwargs["color"] = color_issimple(ob)
    plot_line(ax, ob, **kwargs)

def plot_bounds(ax, ob, zorder=1, alpha=1):
    x, y = zip(*list((p.x, p.y) for p in ob.boundary))
    ax.plot(x, y, 'o', color=BLACK, zorder=zorder, alpha=alpha)

def add_origin(ax, geom, origin):
    x, y = xy = affinity.interpret_origin(geom, origin, 2)
    ax.plot(x, y, 'o', color=GRAY, zorder=1)
    ax.annotate(str(xy), xy=xy, ha='center',
                textcoords='offset points', xytext=(0, 8))

def set_limits(ax, x0, xN, y0, yN):
    ax.set_xlim(x0, xN)
    ax.set_xticks(range(x0, xN+1))
    ax.set_ylim(y0, yN)
    ax.set_yticks(range(y0, yN+1))
    ax.set_aspect("equal")

这是我的代码：

pts_tri = [(23.5, 40.703193977868615), (6.5, 11.258330249197702), (19.650283067421512,12.213053056663625)]
test = (8.258506615, 14.304153051823665)
poly = Polygon(pts_tri)
linearring = LinearRing(list(poly.exterior.coords))
fig = plt.figure(figsize=(10,5))
ax2 = plt.subplot(1,2,1)
for point in pts_tri:
    plt.plot(point[0],point[1],'o',color='black')
plot_coords(ax2, poly.exterior)
patch = PolygonPatch(poly, facecolor=color_isvalid(poly), edgecolor=color_isvalid(poly, valid=BLUE), alpha=0.5, zorder=2)
ax2.add_patch(patch)
plt.plot(test[0],test[1],'o', color='red')
plt.show()
print((poly.contains(Point(test))),(poly.touches(Point(test))),linearring.intersects(Point(test)),poly.overlaps(Point(test)))

plt.show()的结果

而且打印的结果是：False False False False

问题是点位于多边形边界上，但我只能得到False。我该怎么做呢？ 提前感谢。

- Altanir Junior

---

请修改此内容，包含最少可行代码以重现错误。 - Hayden Eastwood

它缺少plt导入。pts_tri是多边形顶点的列表，test是测试点。 - Altanir Junior

我已经看了一下，无法解决 - 明天早上会回来处理。你尝试过在多边形区域更明显地标出点吗？ - Hayden Eastwood

谢谢。我正在使用这段代码来验证岩石分类图上的点，使用了数千个点。当点不在边界上时，shapely测试算法按预期工作。 - Altanir Junior

Answer 1

由于浮点精度问题，Shapely的相交可能会出现意外情况。如果我舍入您的test点...

test = (round(8.258506615, 14), round(14.304153051823665, 14))

然后我可以得到 poly.contains(Point(test)) 为 True 和 poly.intersection(Point(test)) 为 POINT (8.258506615 14.30415305182366)。一个更好的方法可能是只检查距离。

epsilon = 1e-15
print(Point(test).distance(poly) < epsilon)
>>> True

更详细的技术解释可以在这里找到。

- Matt M