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

Question

上次我使用shapely时，我真的有这种美妙的导入和飞翔的感觉。 然而最近，我在这个模块中遇到了一个相当不直观的行为，因为我试图在 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 平面并与线段正交时，该函数的行为与预期一致。然而，当三角形“倾斜”时，它会返回一条线。这暂时让我相信它返回了直线和三角形边界框之间的交点。上面的红线证明事实并非如此。

因此，解决此问题的方法是阅读这个非常有启发性的网页，并调整他们的

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

不再那么进口和飞行......

最后我的问题是：

• 当应用于 3D 对象时，

  intersection

  返回什么？为什么它是一条线？

• shapely中有没有一个函数可以满足我的需求？或任何可选的参数、调整或黑暗魔法？

• 还有其他图书馆可以完成这项工作，同时实现我简单和懒惰的梦想吗？

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)

Answer 2

我找到了这个问题的答案。由于 Shapely 和 Geopandas 默认为 2D，因此它们只是忽略 3 维坐标。相反，我们可以使用 trimesh 来执行相同的操作。有限的工作示例如下。

import trimesh
import numpy as np
#Create a plane for the face 
vertices=[(0,0,0),(4,0,0),(4,4,0),(0,4,0)]
mesh=trimesh.Trimesh(vertices=vertices)
mesh=mesh.convex_hull 
ray_origins = np.array([[2,2,-3],[2,2,100]]) #Create a vertical line
ray_directions=np.array([0,0,1],[0,0,1]])
locations,index_ray,index_tri=mesh.ray.intersects_location(ray_origins,ray_directions)
print(locations)

应给出输出 [2,2,0]。