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.

The most optimized solution that covers all possible cases, would be to use a geometry package, like shapely, scikit-geometry or pygeos. All of them use C++ geometry packages under the hood. The first one is easy to install via pip:

pip install shapely

and simple to use:

from shapely.geometry import Polygon
pgon = Polygon(zip(x, y)) # Assuming the OP's x,y coordinates

print(pgon.area)

To build it from scratch or understand how the underlying algorithm works, check the shoelace formula:

# e.g. corners = [(2.0, 1.0), (4.0, 5.0), (7.0, 8.0)]
def Area(corners):
    n = len(corners) # of corners
    area = 0.0
    for i in range(n):
        j = (i + 1) % n
        area += corners[i][0] * corners[j][1]
        area -= corners[j][0] * corners[i][1]
    area = abs(area) / 2.0
    return area

Since this works for simple polygons:

If you have a polygon with holes : Calculate the area of the outer ring and subtrack the areas of the inner rings
If you have self-intersecting rings : You have to decompose them into simple sectors

By analysis of Mahdi's answer, I concluded that the majority of time was spent doing np.roll(). By removing the need of the roll, and still using numpy, I got the execution time down to 4-5µs per loop compared to Mahdi's 41µs (for comparison Mahdi's function took an average of 37µs on my machine).

def polygon_area(x,y):
    correction = x[-1] * y[0] - y[-1]* x[0]
    main_area = np.dot(x[:-1], y[1:]) - np.dot(y[:-1], x[1:])
    return 0.5*np.abs(main_area + correction)

By calculating the correctional term, and then slicing the arrays, there is no need to roll or create a new array.

Benchmarks:

10000 iterations
PolyArea(x,y): 37.075µs per loop
polygon_area(x,y): 4.665µs per loop

Timing was done using the time module and time.clock()

Calculate area of polygon given (x,y) coordinates

Related

Recent Posts