python - 从 xyz 坐标查找多边形的面积

Question

我正在尝试使用 shapely.geometry.Polygon 模块来查找多边形的面积，但它在 xy 平面上执行所有计算。这对我的一些多边形来说很好，但其他多边形也有 z 维度，所以它不是我想要的。

是否有一个包可以根据 xyz 坐标给出平面多边形的面积，或者有一个包或算法可以将多边形旋转到 xy 平面这样我就可以使用 shapely.geometry.Polygon().area 了？

多边形表示为 [(x1,y1,z1),(x2,y2,z3),...(xn,yn,zn)] 形式的元组列表 .

Answer 1

Here is the derivation of a formula for calculating the area of a 3D planar polygon

这是实现它的 Python 代码:

#determinant of matrix a
def det(a):
    return a[0][0]*a[1][1]*a[2][2] + a[0][1]*a[1][2]*a[2][0] + a[0][2]*a[1][0]*a[2][1] - a[0][2]*a[1][1]*a[2][0] - a[0][1]*a[1][0]*a[2][2] - a[0][0]*a[1][2]*a[2][1]

#unit normal vector of plane defined by points a, b, and c
def unit_normal(a, b, c):
    x = det([[1,a[1],a[2]],
             [1,b[1],b[2]],
             [1,c[1],c[2]]])
    y = det([[a[0],1,a[2]],
             [b[0],1,b[2]],
             [c[0],1,c[2]]])
    z = det([[a[0],a[1],1],
             [b[0],b[1],1],
             [c[0],c[1],1]])
    magnitude = (x**2 + y**2 + z**2)**.5
    return (x/magnitude, y/magnitude, z/magnitude)

#dot product of vectors a and b
def dot(a, b):
    return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]

#cross product of vectors a and b
def cross(a, b):
    x = a[1] * b[2] - a[2] * b[1]
    y = a[2] * b[0] - a[0] * b[2]
    z = a[0] * b[1] - a[1] * b[0]
    return (x, y, z)

#area of polygon poly
def area(poly):
    if len(poly) < 3: # not a plane - no area
        return 0

    total = [0, 0, 0]
    for i in range(len(poly)):
        vi1 = poly[i]
        if i is len(poly)-1:
            vi2 = poly[0]
        else:
            vi2 = poly[i+1]
        prod = cross(vi1, vi2)
        total[0] += prod[0]
        total[1] += prod[1]
        total[2] += prod[2]
    result = dot(total, unit_normal(poly[0], poly[1], poly[2]))
    return abs(result/2)

为了测试它，这里有一个倾斜的 10x5 正方形:

>>> poly = [[0, 0, 0], [10, 0, 0], [10, 3, 4], [0, 3, 4]]
>>> poly_translated = [[0+5, 0+5, 0+5], [10+5, 0+5, 0+5], [10+5, 3+5, 4+5], [0+5, 3+5, 4+5]]
>>> area(poly)
50.0
>>> area(poly_translated)
50.0
>>> area([[0,0,0],[1,1,1]])
0

最初的问题是我过于简单化了。它需要计算垂直于平面的单位向量。面积是点积与所有叉积之和的一半，而不是所有叉积量值之和的一半。

这可以稍微清理一下(矩阵和向量类会使它更好，如果你有它们，或者行列式/叉积/点积的标准实现)，但它在概念上应该是合理的。