从xyz坐标中查找多边形的区域 [英] Find area of polygon from xyz coordinates

问题描述

我正在尝试使用 shapely.geometry.Polygon 模块查找多边形的面积，但是它在 xy 平面上执行所有计算.这对于我的某些多边形很好，但其他多边形也具有 z 尺寸，因此它并不能完全满足我的要求.

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().area ?

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?

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

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

推荐答案

这里是公式的推导计算3D平面多边形的面积

以下是实现它的Python代码:

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.

这篇关于从xyz坐标中查找多边形的区域的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！