问题描述
我正在尝试使用 shapely.geometry.Polygon 模块来查找多边形的面积，但它在 xy 平面上执行所有计算.这对于我的一些多边形来说很好，但其他多边形也有 z 维度，所以它并没有完全按照我的意愿行事.
是否有一个包可以从 xyz 坐标中给我一个平面多边形的面积，或者一个包或算法来将多边形旋转到 xy 平面这样我就可以使用 shapely.geometry.Polygon().area?
多边形表示为形式为 [(x1,y1,z1),(x2,y2,z3),...(xn,yn,zn)] 的元组列表.
这是实现它的 Python 代码:
#矩阵a的行列式定义 det(a):返回 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]#由点a、b和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]])幅度 = (x**2 + y**2 + z**2)**.5返回(x/幅度，y/幅度，z/幅度)#向量a和b的点积定义点(a，b):返回 a[0]*b[0] + a[1]*b[1] + a[2]*b[2]#向量a和b的叉积定义交叉(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]返回 (x, y, z)#多边形poly的面积定义区域(多边形):如果 len(poly) 为了测试它，这里有一个倾斜的 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]]>>>面积(多边形)50.0>>>区域(poly_translated)50.0>>>区域([[0,0,0],[1,1,1]])0问题最初是我过于简单化了.它需要计算垂直于平面的单位向量.面积是其点积与所有叉积之和的一半，而不是所有叉积大小之和的一半.
这可以稍微清理一下(如果你有矩阵和向量类，或者行列式/交叉积/点积的标准实现，矩阵和向量类会使它更好)，但它在概念上应该是合理的.
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.
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?
The polygons are represented as a list of tuples in the form [(x1,y1,z1),(x2,y2,z3),...(xn,yn,zn)].
Here is the derivation of a formula for calculating the area of a 3D planar polygon
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)
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 坐标中查找多边形区域的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！