Calculate Area of Polygon Given (X,Y) Coordinates

Calculate area of polygon given (x,y) coordinates

Implementation of Shoelace formula could be done in Numpy. Assuming these vertices:

import numpy as np
x = np.arange(0,1,0.001)
y = np.sqrt(1-x**2)

We can redefine the function in numpy to find the area:

def PolyArea(x,y):
    return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))

And getting results:

print PolyArea(x,y)
# 0.26353377782163534

Avoiding for loop makes this function ~50X faster than PolygonArea:

%timeit PolyArea(x,y)
# 10000 loops, best of 3: 42 µs per loop
%timeit PolygonArea(zip(x,y))
# 100 loops, best of 3: 2.09 ms per loop.

Timing is done in Jupyter notebook.

How to calculate area of polygon from list of points with python?

use the shapely module available for both Python 2.7 and 3

In [41]: from shapely.geometry import Polygon

In [48]: coords  = ((-1, 0), (-1, 1), (0, 0.5), (1, 1), (1, 0), (-1, 0))

In [49]: polygon = Polygon(coords)

In [50]: polygon.area
Out[50]: 1.5

How do I calculate the area of a 2d polygon?

Here is the standard method, AFAIK. Basically sum the cross products around each vertex. Much simpler than triangulation.

Python code, given a polygon represented as a list of (x,y) vertex coordinates, implicitly wrapping around from the last vertex to the first:

def area(p):
    return 0.5 * abs(sum(x0*y1 - x1*y0
                         for ((x0, y0), (x1, y1)) in segments(p)))

def segments(p):
    return zip(p, p[1:] + [p[0]])

David Lehavi comments: It is worth mentioning why this algorithm works: It is an application of Green's theorem for the functions −y and x; exactly in the way a planimeter works. More specifically:

Formula above =

integral_over_perimeter(-y dx + x dy) =
integral_over_area((-(-dy)/dy+dx/dx) dy dx) =
2 Area

How to calculate the area of a polygon on the earth's surface using python?

Let's say you have a representation of the state of Colorado in GeoJSON format

{"type": "Polygon", 
 "coordinates": [[
   [-102.05, 41.0], 
   [-102.05, 37.0], 
   [-109.05, 37.0], 
   [-109.05, 41.0]
 ]]}

All coordinates are longitude, latitude. You can use pyproj to project the coordinates and Shapely to find the area of any projected polygon:

co = {"type": "Polygon", "coordinates": [
    [(-102.05, 41.0),
     (-102.05, 37.0),
     (-109.05, 37.0),
     (-109.05, 41.0)]]}
lon, lat = zip(*co['coordinates'][0])
from pyproj import Proj
pa = Proj("+proj=aea +lat_1=37.0 +lat_2=41.0 +lat_0=39.0 +lon_0=-106.55")

That's an equal area projection centered on and bracketing the area of interest. Now make new projected GeoJSON representation, turn into a Shapely geometric object, and take the area:

x, y = pa(lon, lat)
cop = {"type": "Polygon", "coordinates": [zip(x, y)]}
from shapely.geometry import shape
shape(cop).area  # 268952044107.43506

It's a very close approximation to the surveyed area. For more complex features, you'll need to sample along the edges, between the vertices, to get accurate values. All caveats above about datelines, etc, apply. If you're only interested in area, you can translate your feature away from the dateline before projecting.

Find area of polygon from xyz coordinates

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.

How calculate the area of a polygon in R, if I only have X,Y coordinates?

This is a convex hull problem. Other questions cover similar territory. We can gather the plotting and area calculation problems here for good measure.

To plot the convex hull of a point cloud, we can use the chull function:

library(tidyverse)

data <- tibble(x = runif(50), y = runif(50))

convex_hull <- data %>% slice(chull(x, y))

ggplot(data, aes(x = x, y = y)) + 
  geom_point() + 
  geom_polygon(data = convex_hull,
               alpha = 0.25)

Library splancs has an areapl function to calculate the area of a non-selfintersecting polygon, such as the convex_hull above. Hence, we may do:

library(splancs)
as.matrix(convex_hull) %>% areapl

---

Related Topics