41

The title basically says it all. I need to calculate the area inside a polygon on the Earth's surface using Python. Calculating area enclosed by arbitrary polygon on Earth's surface says something about it, but remains vague on the technical details:

If you want to do this with a more "GIS" flavor, then you need to select an unit-of-measure for your area and find an appropriate projection that preserves area (not all do). Since you are talking about calculating an arbitrary polygon, I would use something like a Lambert Azimuthal Equal Area projection. Set the origin/center of the projection to be the center of your polygon, project the polygon to the new coordinate system, then calculate the area using standard planar techniques.

So, how do I do this in Python?

1
  • If you choose an area preserving projection, the sides of the polygon will no longer be straight Jan 13, 2011 at 20:13

10 Answers 10

43

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.

7
  • 6
    Strictly speaking, that GeoJSON should have a fifth closing point, [-102.05, 41.0]. The spec is a bit vague on these things sometimes.
    – Brad Koch
    Jul 7, 2014 at 18:00
  • 7
    what's the unit of that area result? is it in square meters or square inches?
    – jyf1987
    Nov 9, 2016 at 7:42
  • 4
    In case anyone else is wondering, it is in square meters. You can google area of colorado and get 104,185 square miles. Converting that to square meters gives roughly the same result.
    – spanishgum
    Feb 7, 2017 at 21:38
  • 1
    Does this work just as well for polygons that might not be simple looking quadrilaterals?
    – Ali
    Jun 9, 2017 at 9:36
  • So, basically I can use any equal-area projection to convert the coordinates and use Green's theorem to compute the area, am I right? If so, then one can alternatively use basemap to do the projection (e.g. 'aea', 'cea'). In terms of Green's theorem, a one-liner area=np.abs(0.5*np.sum(y[:-1]*np.diff(x) - x[:-1]*np.diff(y))) frees you the need of the shapely module
    – Jason
    Dec 23, 2017 at 9:31
27

The easiest way to do this (in my opinion), is to project things into (a very simple) equal-area projection and use one of the usual planar techniques for calculating area.

First off, I'm going to assume that a spherical earth is close enough for your purposes, if you're asking this question. If not, then you need to reproject your data using an appropriate ellipsoid, in which case you're going to want to use an actual projection library (everything uses proj4 behind the scenes, these days) such as the python bindings to GDAL/OGR or (the much more friendly) pyproj.

However, if you're okay with a spherical earth, it quite simple to do this without any specialized libraries.

The simplest equal-area projection to calculate is a sinusoidal projection. Basically, you just multiply the latitude by the length of one degree of latitude, and the longitude by the length of a degree of latitude and the cosine of the latitude.

def reproject(latitude, longitude):
    """Returns the x & y coordinates in meters using a sinusoidal projection"""
    from math import pi, cos, radians
    earth_radius = 6371009 # in meters
    lat_dist = pi * earth_radius / 180.0

    y = [lat * lat_dist for lat in latitude]
    x = [long * lat_dist * cos(radians(lat)) 
                for lat, long in zip(latitude, longitude)]
    return x, y

Okay... Now all we have to do is to calculate the area of an arbitrary polygon in a plane.

There are a number of ways to do this. I'm going to use what is probably the most common one here.

def area_of_polygon(x, y):
    """Calculates the area of an arbitrary polygon given its verticies"""
    area = 0.0
    for i in range(-1, len(x)-1):
        area += x[i] * (y[i+1] - y[i-1])
    return abs(area) / 2.0

Hopefully that will point you in the right direction, anyway...

3
  • Note that lines will be distorted with a sinusoidal projection. It might be wise to interpolate along great circles between points so that your projected polygon doesn't get distorted.
    – Spacedman
    Jan 13, 2011 at 17:17
  • 2
    @spacedman - Of course! I was just trying to show a simple approximation. It's not intended to be entirely accurate. The distortion of the sides of the polygon will only matter if he's trying to calculate a polygon with very large sides and few verticies, though. Jan 13, 2011 at 17:21
  • Thank you so much for that reprojection function, you saved me after a day of being lost! Side note to others: after the reprojection, you can also use Python's shapely package to compute the area with shapely.geometry.Polygon(list(zip(x, y))).area to get the area in sq. meters.
    – Ali
    Jun 9, 2017 at 11:04
9

A bit late perhaps, but here is a different method, using Girard's theorem. It states that the area of a polygon of great circles is R**2 times the sum of the angles between the polygons minus (N-2)*pi where N is number of corners.

I thought this would be worth posting, since it doesn't rely on any other libraries than numpy, and it is a quite different method than the others. Of course, this only works on a sphere, so there will be some inaccuracy when applying it to the Earth.

First, I define a function to compute the bearing angle from point 1 along a great circle to point 2:

import numpy as np
from numpy import cos, sin, arctan2

d2r = np.pi/180

def greatCircleBearing(lon1, lat1, lon2, lat2):
    dLong = lon1 - lon2

    s = cos(d2r*lat2)*sin(d2r*dLong)
    c = cos(d2r*lat1)*sin(d2r*lat2) - sin(lat1*d2r)*cos(d2r*lat2)*cos(d2r*dLong)

    return np.arctan2(s, c)

Now I can use this to find the angles, and then the area (In the following, lons and lats should of course be specified, and they should be in the right order. Also, the radius of the sphere should be specified.)

N = len(lons)

angles = np.empty(N)
for i in range(N):

    phiB1, phiA, phiB2 = np.roll(lats, i)[:3]
    LB1, LA, LB2 = np.roll(lons, i)[:3]

    # calculate angle with north (eastward)
    beta1 = greatCircleBearing(LA, phiA, LB1, phiB1)
    beta2 = greatCircleBearing(LA, phiA, LB2, phiB2)

    # calculate angle between the polygons and add to angle array
    angles[i] = np.arccos(cos(-beta1)*cos(-beta2) + sin(-beta1)*sin(-beta2))

area = (sum(angles) - (N-2)*np.pi)*R**2

With the Colorado coordinates given in another reply, and with Earth radius 6371 km, I get that the area is 268930758560.74808

6
  • I tried your method and got a negative value. Besides, its abs (139013699.103) is quite different from the value (809339.212) I got from using a "grid" method, which is taking all the grid cells inside the polygon from a rectangular Mercator grid, and summing up the grid cell areas. The area of each cell is the product of it zonal and meridional intervals (converted from lat,lon to km).
    – Jason
    Dec 23, 2017 at 7:58
  • Could you clarify your comment? Where are the numbers you refer to coming from? I tried the code I posted again, and it gives the reported result...
    – sulkeh
    Dec 23, 2017 at 13:30
  • I applied it on a contour found using matplotlib's contour function, in fact I have a few tens of them and everyone reported a negative area using your method.
    – Jason
    Dec 23, 2017 at 15:08
  • Ok. I still don't understand where the numbers that you mention come from. I would happy to fix any mistake if you can be more specific. Can you post an example of lons/lats that give a negative area?
    – sulkeh
    Dec 23, 2017 at 15:32
  • I exported one contour coordinates here pastebin.com/2fMkFgvu maybe you could have a try. One thing I noticed that if I sub-sample the 135-point contour by taking every 5 points (x=x[::5]; y=y[::5]), I could get a positive value that is close to the 'true' value.
    – Jason
    Dec 23, 2017 at 15:55
9

Or simply use a library: https://github.com/scisco/area

from area import area
>>> obj = {'type':'Polygon','coordinates':[[[-180,-90],[-180,90],[180,90],[180,-90],[-180,-90]]]}
>>> area(obj)
511207893395811.06

...returns the area in square meters.

2
  • 1
    I like that this method lets me use latitude/longitude points directly, and accounts for the curvature of the earth without needing to project the points on a Cartesian grid. I compared the results from this method and the standard Green's theorem to the areas of different states on Google and found that this library gave the most accurate results. Thanks!
    – blaylockbk
    Mar 8, 2019 at 22:35
  • 2
    AFAI understand, this library is calculating geodesic area assuming the earth is a perfect sphere with radius defined as WGS84_RADIUS = 6378137 Just keep that in mind.
    – rysqui
    Nov 14, 2019 at 9:55
6

You can compute the area directly on the sphere, instead of using an equal-area projection.

Moreover, according to this discussion, it seems that Girard's theorem (sulkeh's answer) does not give accurate results in certain cases, for example "the area enclosed by a 30º lune from pole to pole and bounded by the prime meridian and 30ºE" (see here).

A more precise solution would be to perform line integral directly on the sphere. The comparison below shows this method is more precise.

Like all other answers, I should mention the caveat that we assume a spherical earth, but I assume that for non-critical purposes this is enough.

Python implementation

Here is a Python 3 implementation which uses line integral and Green's theorem:

def polygon_area(lats, lons, radius = 6378137):
    """
    Computes area of spherical polygon, assuming spherical Earth. 
    Returns result in ratio of the sphere's area if the radius is specified.
    Otherwise, in the units of provided radius.
    lats and lons are in degrees.
    """
    from numpy import arctan2, cos, sin, sqrt, pi, power, append, diff, deg2rad
    lats = np.deg2rad(lats)
    lons = np.deg2rad(lons)

    # Line integral based on Green's Theorem, assumes spherical Earth

    #close polygon
    if lats[0]!=lats[-1]:
        lats = append(lats, lats[0])
        lons = append(lons, lons[0])

    #colatitudes relative to (0,0)
    a = sin(lats/2)**2 + cos(lats)* sin(lons/2)**2
    colat = 2*arctan2( sqrt(a), sqrt(1-a) )

    #azimuths relative to (0,0)
    az = arctan2(cos(lats) * sin(lons), sin(lats)) % (2*pi)

    # Calculate diffs
    # daz = diff(az) % (2*pi)
    daz = diff(az)
    daz = (daz + pi) % (2 * pi) - pi

    deltas=diff(colat)/2
    colat=colat[0:-1]+deltas

    # Perform integral
    integrands = (1-cos(colat)) * daz

    # Integrate 
    area = abs(sum(integrands))/(4*pi)

    area = min(area,1-area)
    if radius is not None: #return in units of radius
        return area * 4*pi*radius**2
    else: #return in ratio of sphere total area
        return area

I wrote a somewhat more explicit version (and with many more references and TODOs...) in the sphericalgeometry package there.

Numerical Comparison

Colorado will be the reference, since all previous answers were evaluated on its area. Its precise total area is 104,093.67 square miles (from the US Census Bureau, p. 89, see also here), or 269601367661 square meters. I found no source for the actual methodology of the USCB, but I assume it is based on summing actual measurements on ground, or precise computations using WGS84/EGM2008.

Method                 | Author     | Result       | Variation from ground truth
--------------------------------------------------------------------------------
Albers Equal Area      | sgillies   | 268952044107 | -0.24%
Sinusoidal             | J. Kington | 268885360163 | -0.26%
Girard's theorem       | sulkeh     | 268930758560 | -0.25%
Equal Area Cylindrical | Jason      | 268993609651 | -0.22%
Line integral          | Yellows    | 269397764066 | **-0.07%**

Conclusion: using direct integral is more precise.

Performance

I have not benchmarked the different methods, and comparing pure Python code with compiled PROJ projections would not be meaningful. Intuitively less computations are needed. On the other hand, trigonometric functions may be computationally intensive.

2
  • I now prefer this solution over my own, for its not relying on any extra packages. Basemap has been deprecated.
    – Jason
    Jun 8, 2021 at 3:25
  • The border of Colorado is not a perfect square, there is a small "bulge" on the western side from imperfect surveys when the border was establish. The census area probably accounts for this and is why all areas here underestimate it slightly.
    – Shawn
    Oct 6, 2021 at 0:36
6

Here is a solution that uses basemap, instead of pyproj and shapely, for the coordinate conversion. The idea is the same as suggested by @sgillies though. NOTE that I've added the 5th point so that the path is a closed loop.

import numpy
from mpl_toolkits.basemap import Basemap

coordinates=numpy.array([
[-102.05, 41.0], 
[-102.05, 37.0], 
[-109.05, 37.0], 
[-109.05, 41.0],
[-102.05, 41.0]])

lats=coordinates[:,1]
lons=coordinates[:,0]

lat1=numpy.min(lats)
lat2=numpy.max(lats)
lon1=numpy.min(lons)
lon2=numpy.max(lons)

bmap=Basemap(projection='cea',llcrnrlat=lat1,llcrnrlon=lon1,urcrnrlat=lat2,urcrnrlon=lon2)
xs,ys=bmap(lons,lats)

area=numpy.abs(0.5*numpy.sum(ys[:-1]*numpy.diff(xs)-xs[:-1]*numpy.diff(ys)))
area=area/1e6

print area

The result is 268993.609651 in km^2.

UPDATE: Basemap has been deprecated, so you may want to consider alternative solutions first.

2
  • For clarification, what is the source for that equation used to compute the area?
    – blaylockbk
    Feb 27, 2019 at 22:39
  • 1
    @blaylockbk It's Green's theorem.
    – Jason
    Feb 28, 2019 at 12:08
4

Because the earth is a closed surface a closed polygon drawn on its surface creates TWO polygonal areas. You also need to define which one is inside and which is outside!

Most times people will be dealing with small polygons, and so it's 'obvious' but once you have things the size of oceans or continents, you better make sure you get this the right way round.

Also, remember that lines can go from (-179,0) to (+179,0) in two different ways. One is very much longer than the other. Again, mostly you'll make the assumption that this is a line that goes from (-179,0) to (-180,0) which is (+180,0) and then to (+179,0), but one day... it won't.

Treating lat-long like a simple (x,y) coordinate system, or even neglecting the fact that any coordinate projection is going to have distortions and breaks, can make you fail big-time on spheres.

2
  • Very true! I was just giving a simple example of the most basic case in my answer... It doesn't even remotely try to deal with crossing the prime meridian (in 0-360) or the antimeridian (in -180 to 180). Jan 13, 2011 at 17:17
  • 20
    Yeah. Life was so much easier when the earth was flat.
    – Spacedman
    Jan 13, 2011 at 17:43
4

I know that answering 10 years later has some advantages, but to somebody that looks today at this question it seems fair to provide an updated answer.

pyproj directly calculates areas, without need of calling shapely:

# Modules:
from pyproj import Geod
import numpy as np

# Define WGS84 as CRS:
geod = Geod('+a=6378137 +f=0.0033528106647475126')

# Data for Colorado (no need to close the polygon):
coordinates = np.array([
[-102.05, 41.0], 
[-102.05, 37.0], 
[-109.05, 37.0], 
[-109.05, 41.0]])
lats = coordinates[:,1]
lons = coordinates[:,0]

# Compute:
area, perim = geod.polygon_area_perimeter(lons, lats)

print(abs(area))  # Positive is counterclockwise, the data is clockwise.

The result is: 269154.54988400977 km2, or -0.17% of the reported correct value (269601.367661 km2).

0

According to Yellows' assertion, direct integral is more precise.

But Yellows use an earth radius = 6378 137m, which is the WGS-84 ellipsoid, semi-major axis, while Sulkeh use 6371 000 m.

Using a radius = 6378 137 m in the Sulkeh' method, gives 269533625893 square meters.

Assuming that the true value of Colorado area (from the US Census Bureau) is 269601367661 square meters then the variation from the ground truth of Sulkeh' method is : -0,025%, better than -0.07 with the Line integral method.

So Sulkeh' proposal seems to be the more precise so far.

In order to be able to make a numerical comparison of the solutions, with the assumption of a spherical Earth, all calculations must use the same terrestrial radius.

0

Here is a Python 3 implementation where the function would take a list of tuple-pairs of lats and longs and would return the area enclosed in the projected polygon.It uses pyproj to project the coordinates and then Shapely to find the area of any projected polygon

def calc_area(lis_lats_lons):

import numpy as np
from pyproj import Proj
from shapely.geometry import shape


lons, lats = zip(*lis_lats_lons)
ll = list(set(lats))[::-1]
var = []
for i in range(len(ll)):
    var.append('lat_' + str(i+1))
st = ""
for v, l in zip(var,ll):
    st = st + str(v) + "=" + str(l) +" "+ "+"
st = st +"lat_0="+ str(np.mean(ll)) + " "+ "+" + "lon_0" +"=" + str(np.mean(lons))
tx = "+proj=aea +" + st
pa = Proj(tx)

x, y = pa(lons, lats)
cop = {"type": "Polygon", "coordinates": [zip(x, y)]}

return shape(cop).area 

For a sample set of lats/longs, it gives an area value close to the surveyed approximation value

calc_area(lis_lats_lons = [(-102.05, 41.0),
 (-102.05, 37.0),
 (-109.05, 37.0),
 (-109.05, 41.0)])

Which outputs an area of 268952044107.4342 Sq. Mts.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.