I'm trying to use the shapely.geometry.Polygon module to find the area of polygons but it performs all calculations on the xy plane. This is fine for some of my polygons but others have a z dimension too so it's not quite doing what I'd like.

有一个包，它要么给我一个平面多边形的面积从 XYZ 坐标，或者一包或算法来旋转多边形的 XY 面，这样我就可以使用 shapely.geometry.Polygon（）。面积？

Is there a package which will either give me the area of a planar polygon from xyz coordinates, or alternatively a package or algorithm to rotate the polygon to the xy plane so that i can use shapely.geometry.Polygon().area?

的多边形重新psented作为元组形式的列表$ P $ [（X1，Y1，Z1），（X2，Y2，Z3），...（XN，YN，锌）] 。

The polygons are represented as a list of tuples in the form [(x1,y1,z1),(x2,y2,z3),...(xn,yn,zn)].

推荐答案

这里有一个公式的推导用于计算3D平面多边形的领域

下面是Python的code实现它：

Here is Python code that implements it:

#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正方形凑到：

And to test it, here's a 10x5 square that leans over:

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

这个问题本来是我过于简单化。它需要计算单位矢量与平面垂直。面积为一半的点积的和所有的交叉乘积的总和，不一半的横产品的所有的大小的总和。

The problem originally was that I had oversimplified. It needs to calculate the unit vector normal to the plane. The area is half of the dot product of that and the total of all the cross products, not half of the sum of all the magnitudes of the cross products.

这可以清理一个位（矩阵，向量类将使它更好，如果有的话，或决定性/跨产品/点产品的标准实现），但它应该是概念上的声音。

This can be cleaned up a bit (matrix and vector classes would make it nicer, if you have them, or standard implementations of determinant/cross product/dot product), but it should be conceptually sound.