Shapely：3D 中线和多边形之间的交点答案

【问题标题】：Shapely: intersection point between line and polygon in 3DShapely：3D 中线和多边形之间的交点
【发布时间】：2018-05-01 18:45:31
【问题描述】：

上次我用shapely，我真的有this nice import-and-fly-feeling。然而最近，我在这个模块中遇到了一个相当不直观的行为，因为我试图在 3D 空间中找到线段和三角形之间的交点。让我们定义一个线段和一个三角形如下：

l = LineString([[1,0.5,0.5],[3,0.5,0.5]])
p = Polygon([[1.2,0.0,0.],[2.2,1.0,0.],[2.8,0.5,1.]])

为了得到他们的交点，我使用了l.intersection(p)，并期望得到一个点，即POINT Z (POINT Z (2 0.5 0.25))。如下图蓝点所示：

相反，我得到的是LINESTRING Z (1.7 0.5 0.25, 2.8 0.5 1) - 下面的红线 - 坦率地说，我对它应该代表什么感到非常困惑。

奇怪的是，当多边形/三角形位于 xz 平面并与线段正交时，函数的行为与预期的一样。但是，当三角形“倾斜”时，它会返回一条线。这让我暂时相信它返回了线和三角形边界框之间的交点。上面的红线证明不是这样。

因此，解决此问题的方法是阅读this very enlightening web-page，并调整他们的C++ 代码以处理形状良好的对象。 intersection 方法可以很好地检查线是否穿过多边形，下面的函数会找到感兴趣的点。

def intersect3D_SegmentPlane(Segment, Plane):

    # Points in Segment: Pn  Points in Plane: Qn
    P0, P1     = np.array(Segment.coords)
    Q0, Q1, Q2 = np.array(Plane.exterior)[:-1]   

    # vectors in Plane
    q1 = Q1 - Q0
    q2 = Q2 - Q0

    # vector normal to Plane
    n  = np.cross(q1, q2)/np.linalg.norm(np.cross(q1, q2))
    u = P1 - P0 # Segment's direction vector 
    w = P0 - Q0 # vector from plane ref point to segment ref point

    ## Tests parallelism
    if np.dot(n, u) == 0:
        print "Segment and plane are parallel"
        print "Either Segment is entirely in Plane or they never intersect."
        return None
    ## if intersection is a point
    else:
        ## Si is the scalar where P(Si) = P0 + Si*u lies in Plane
        Si = np.dot(-n, w) / np.dot(n, u)
        PSi = P0 + Si * u
        return PSi

不再是非常进口和飞行......

所以最后我的问题：

intersection 应用于 3D 对象时返回什么？为什么它是一条线？
shapely 中是否有一个功能可以满足我的要求？或任何可选参数、调整或黑魔法？
是否还有其他库可以完成这项工作，同时实现我对简单和懒惰的梦想？

【问题讨论】：

标签： python vector geometry shapely

【解决方案1】：

不幸的是，正如文档所述：

坐标序列是不可变的。构造实例时可以使用第三个 z 坐标值，但对几何分析没有影响。所有操作都在 x-y 平面上进行。

可以通过以下方式验证这一点：

from shapely.geometry import LineString, Polygon

l = LineString([[1,0.5,0.5],[3,0.5,0.5]])
p = Polygon([[1.2,0.0,0.],[2.2,1.0,0.],[2.8,0.5,1.]])
print(l.intersection(p))
#LINESTRING Z (1.7 0.5 0.25, 2.8 0.5 1)

l = LineString([[1,0.5],[3,0.5]])
p = Polygon([[1.2,0.0],[2.2,1.0],[2.8,0.5]])
print(l.intersection(p))
#LINESTRING (1.7 0.5, 2.8 0.5)

甚至：

from shapely.geometry import LineString, Polygon

l = LineString([[1,0.5,0],[3,0.5,0]])
p = Polygon([[1.2,0.0,1],[2.2,1.0,1],[2.8,0.5,1]])
print(l.intersects(p))
#True (even though the objects are in different z-planes)

【讨论】：