18

モジュールを使用してshapely.geometry.Polygonポリゴンの領域を見つけようとしていますが、xy平面上ですべての計算を実行します。これは、一部のポリゴンでは問題ありませんが、他のポリゴンにもzディメンションがあるため、希望どおりに機能していません。

座標から平面ポリゴンの面積を取得するパッケージ、または使用できるようにポリゴンを平面xyzに回転させるパッケージまたはアルゴリズムのいずれかを提供するパッケージはありますか?xyshapely.geometry.Polygon().area

ポリゴンは、形式のタプルのリストとして表されます[(x1,y1,z1),(x2,y2,z3),...(xn,yn,zn)]。

4

7 に答える 7

23

これは、3D 平面ポリゴンの面積を計算する式の導出です。

これを実装する 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

問題はもともと、私が単純化しすぎたことでした。平面に垂直な単位ベクトルを計算する必要があります。面積は、それとすべての外積の合計の内積の半分であり、外積のすべての大きさの合計の半分ではありません。

これは少しクリーンアップできます (行列クラスとベクトル クラスがある場合、または行列式/クロス積/ドット積の標準実装がある場合は、より適切になります) が、概念的には健全なはずです。

于 2012-09-28T15:53:07.680 に答える
8

これは私が使用した最終的なコードです。shapely は使用しませんが、ストークの定理を実装して面積を直接計算します。これは、numpy なしでそれを行う方法を示す @Tom Smilack の回答に基づいています。

import numpy as np

#unit normal vector of plane defined by points a, b, and c
def unit_normal(a, b, c):
    x = np.linalg.det([[1,a[1],a[2]],
         [1,b[1],b[2]],
         [1,c[1],c[2]]])
    y = np.linalg.det([[a[0],1,a[2]],
         [b[0],1,b[2]],
         [c[0],1,c[2]]])
    z = np.linalg.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)

#area of polygon poly
def poly_area(poly):
    if len(poly) < 3: # not a plane - no area
        return 0
    total = [0, 0, 0]
    N = len(poly)
    for i in range(N):
        vi1 = poly[i]
        vi2 = poly[(i+1) % N]
        prod = np.cross(vi1, vi2)
        total[0] += prod[0]
        total[1] += prod[1]
        total[2] += prod[2]
    result = np.dot(total, unit_normal(poly[0], poly[1], poly[2]))
    return abs(result/2)
于 2012-09-29T15:02:32.627 に答える
0

2D ポリゴンの面積は、Numpy をワンライナーとして使用して計算できます...

poly_Area(vertices) = np.sum( [0.5, -0.5] * vertices * np.roll( np.roll(vertices, 1, axis=0), 1, axis=1) )
于 2013-07-28T09:55:57.857 に答える