Find area of polygon from xyz coordinates

Question

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)].

Allinclusive

a polygon is a strictly 2 dimensional figure. What exactly are you trying to calculate? – Racoon

I'm trying to find the surface areas of roofs and walls of a building from 'xyz' coordinates of the vertices. – Allinclusive

I haven't found any module to do that, but you could simply cast down each face, to an xy plane, and calculate that with the module you have been using – Racoon

What do you mean by "cast down"? – Allinclusive

Just rotate the shape until it's flat on the z plane. – Bifocals

Answer 1

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.

Margret

Thanks, Tom. I'd found that page and also some sample code for applying Stoke's theorem to a 2D polygon but was having trouble making it work for 3D. Your implementation looks good to me. I'm just adapting it to work with the way my data's structured which is [(x1,y1,z1),(x2,y2,z2),...]. – Allinclusive

The area function should be the same. cross_product_magnitude would change to x = a[1] * b[2] - a[2] * b[1] etc. – Margret

Yes, I've got that - but it's throwing out results that are way too large. Do I need to move the shape so one vertex is at the origin? – Allinclusive

You shouldn't have to. I think I messed up somewhere, I'll look into it. – Margret

Did my edit not work then? It works in my code here. Also, your code there is still missing the abs(). I've added another answer with the edited code, and using numpy for the cross-product. – Allinclusive

Actually my edit didn't get the right result on your tests. It worked for walls in my application but not for sloping surfaces. Yours is working just fine. I just need to add abs() when calling it as I'm summing up all the surface areas and some of the areas come out negative: poly_area = abs(area(poly)). – Allinclusive

Added abs to my answer for completeness. It has to do with Stokes' Theorem and how areas can be signed depending on the direction you integrate around the perimeter. – Margret

unit_normal and det are not needed. It is enough to get norm(total) (much simpler to compute; norm should be written). – Pavla

Why is the unit normal computed via determinant? Can’t you just do a cross product of the polygon first two edges + normalization ? – Plagioclase

Answer 2

This is the final code I've used. It doesn't use shapely, but implements Stoke's theorem to calculate the area directly. It builds on @Tom Smilack's answer which shows how to do it without numpy.

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)

Allinclusive

I am looking to implement this solution but what is not clear is why the unit_normal function implements the first 3 points of the polygon. poly is a list of 3d points i.e., a list of tuples as posted in the original question. or is the response applicabel to a 3-point polygon only? thanks – Solo

From what I remember, the unit normal vector is the same for any three (non-colinear) points on a polygon, we can just take the first three points and calculate it from that – Allinclusive

Answer 3

#pythonn code for polygon area in 3D (optimised version)

def polygon_area(poly):
    #shape (N, 3)
    if isinstance(poly, list):
        poly = np.array(poly)
    #all edges
    edges = poly[1:] - poly[0:1]
    # row wise cross product
    cross_product = np.cross(edges[:-1],edges[1:], axis=1)
    #area of all triangles
    area = np.linalg.norm(cross_product, axis=1)/2
    return sum(area)



if __name__ == "__main__":
    poly = [[0+5, 0+5, 0+5], [10+5, 0+5, 0+5], [10+5, 3+5, 4+5], [0+5, 3+5, 4+5]]
    print(polygon_area(poly))

Henriettahenriette

Doesn't work with this "L" shaped polygon [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 0.7071067811865475, 0.7071067811865476], [2.0, 0.7071067811865475, 0.7071067811865476], [2.0, 1.414213562373095, 1.4142135623730951], [0.0, 1.414213562373095, 1.4142135623730951]] Are you sure this will work for concave polygons? – Desalinate

Answer 4

The area of a 2D polygon can be calculated using Numpy as a one-liner...

poly_Area(vertices) = np.sum( [0.5, -0.5] * vertices * np.roll( np.roll(vertices, 1, axis=0), 1, axis=1) )

Versatile

This doesn't work for a 2D polygon in 3D space, e.g. all coplanar but referred to in xyz coordinates. – Allinclusive

Answer 5

Fyi, here is the same algorithm in Mathematica, with a baby unit test

ClearAll[vertexPairs, testPoly, area3D, planeUnitNormal, pairwise];
pairwise[list_, fn_] := MapThread[fn, {Drop[list, -1], Drop[list, 1]}];
vertexPairs[Polygon[{points___}]] := Append[{points}, First[{points}]];
testPoly = Polygon[{{20, -30, 0}, {40, -30, 0}, {40, -30, 20}, {20, -30, 20}}];
planeUnitNormal[Polygon[{points___}]] :=
  With[{ps = Take[{points}, 3]},
   With[{p0 = First[ps]},
    With[{qs = (# - p0) & /@ Rest[ps]},
     Normalize[Cross @@ qs]]]];
area3D[p : Polygon[{polys___}]] :=
  With[{n = planeUnitNormal[p], vs = vertexPairs[p]},
   With[{areas = (Dot[n, #]) & /@ pairwise[vs, Cross]},
    Plus @@ areas/2]];
area3D[testPoly]

Submultiple

The planeUnitNormal computation is not robust in case the first three points are colinear. A smarter algorithm would pick three points that are not colinear (testing by pairwise[...,Cross]=!=0 and throw if it can't find three. – Submultiple

@reb-cabin why throw? If every triple of points is collinear, then the answer is zero. – Mccahill

Answer 6

Same as @Tom Smilack's answer, but in javascript

//determinant of matrix a
function 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
function unit_normal(a, b, c) {
  let x = math.det([
    [1, a[1], a[2]],
    [1, b[1], b[2]],
    [1, c[1], c[2]]
  ]);
  let y = math.det([
    [a[0], 1, a[2]],
    [b[0], 1, b[2]],
    [c[0], 1, c[2]]
  ]);
  let z = math.det([
    [a[0], a[1], 1],
    [b[0], b[1], 1],
    [c[0], c[1], 1]
  ]);
  let magnitude = Math.pow(Math.pow(x, 2) + Math.pow(y, 2) + Math.pow(z, 2), 0.5);
  return [x / magnitude, y / magnitude, z / magnitude];
}
// dot product of vectors a and b
function dot(a, b) {
  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
// cross product of vectors a and b
function cross(a, b) {
  let x = (a[1] * b[2]) - (a[2] * b[1]);
  let y = (a[2] * b[0]) - (a[0] * b[2]);
  let z = (a[0] * b[1]) - (a[1] * b[0]);
  return [x, y, z];
}

// area of polygon poly
function area(poly) {
  if (poly.length < 3) {
    console.log("not a plane - no area");
    return 0;
  } else {
    let total = [0, 0, 0]
    for (let i = 0; i < poly.length; i++) {
      var vi1 = poly[i];
      if (i === poly.length - 1) {
        var vi2 = poly[0];
      } else {
        var vi2 = poly[i + 1];
      }
      let prod = cross(vi1, vi2);
      total[0] = total[0] + prod[0];
      total[1] = total[1] + prod[1];
      total[2] = total[2] + prod[2];
    }
    let result = dot(total, unit_normal(poly[0], poly[1], poly[2]));

    return Math.abs(result/2);
  }

}

Theotokos

"math.det" should just be "det" – Interloper

Answer 7

Thanks for detailed answers, But I am little surprised there is no simple answer to get the area.

So, I am just posting a simplified approach for calculating area using 3d Coordinates of polygon or surface using pyny3d.

#Install pyny3d as:
pip install pyny3d

#Calculate area
import numpy as np
import pyny3d.geoms as pyny

coords_3d = np.array([[0,  0, 0],
                           [7,  0, 0],
                           [7, 10, 2],
                           [0, 10, 2]])
polygon = pyny.Polygon(coords_3d)
print(f'Area is : {polygon.get_area()}')

Upbeat