Shapely：在三维空间中线段和多边形的交点

Question

上次我使用shapely时，真的感觉很好，就像轻松导入和使用一样。然而最近，在尝试在三维空间中找到线段和三角形之间的交点时，我遇到了一个相当不直观的模块行为。让我们定义一个线段和一个三角形如下：

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平面上且垂直于线段时，函数的行为符合预期。然而，当三角形“倾斜”时，它返回一条线。这暂时让我相信它返回了线与三角形边界框之间的交点。上面的红线证明了不同。

因此，解决这个问题的方法是阅读这个非常有启发性的网页，并调整他们的C++代码以适应Shapely对象。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

不再很重要和实时的......

所以最终我的问题是：

• 当应用于3D对象时，intersection返回什么，为什么是一条线？

• 在shapely中是否有一个可以满足我需求的函数？或者任何可选参数、调整或黑魔法技巧？

• 是否有其他库可以完成这个工作，同时满足我简单和懒惰的梦想？

- snake_charmer

Answer 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)

- ewcz

---

现在是2022年，这个答案仍然有效 :-( - Stanley F.