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 Projpa = 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 shapeshape(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.

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

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 npfrom numpy import cos, sin, arctan2d2r = np.pi/180def 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