# How to compute geodesic area in GeoDjango?

I have some polygons which can be anywhere in the world. I would like to compute their areas.

The GeoDjango docs for GEOSGeometry.area don't specify whether it is geometric or geodesic area, but given that GEOSGeometry.distance() says that "GEOS does not perform a spherical calculation even if the SRID specifies a geographic coordinate system", I suspect that it is a simple geometric area.

I'm willing to write a little code, but the shorter the better, of course.

We do have a working geodesic distance function already.

## 4 Answers

if you use PostGIS Geography Type for your table, you can calculate your area as you calculate on the plane surface.

The geography type provides native support for spatial features represented on "geographic" coordinates (sometimes called "geodetic" coordinates, or "lat/lon", or "lon/lat"). Geographic coordinates are spherical coordinates expressed in angular units (degrees).

creating geography(not geometry) table:

CREATE TABLE mypoly (
id SERIAL PRIMARY KEY,
name VARCHAR(64),
the_geom GEOGRAPHY(POINT,4326)
);


hereafter is that using psycopg2 for reaching postgis table and querying some sql with ST_Area.

ST_Area — Returns the area of the surface if it is a polygon or multi-polygon. For "geometry" type area is in SRID units. For "geography" area is in square meters.

in your view.py, write this code:

import psycopg2

def getGeodesicArea:

res = []
con = psycopg2.connect("dbname='mydb' user='reid' host='127.0.0.1' password='reid'")
cur = con.cursor()
cur.execute("SELECT ST_Area(the_geom) FROM mypoly")
rows = cur.fetchall()

for row in rows:
respre = row # find your result.
res.append(respre)

con.commit()
cur.close()
con.close()


if you are searching a pythonic way for fetching your result as json, check out this SO answer.

UPDATE

There is a javascript code in openlayers api for calculating geodetic areas. i think you can convert it to python for django application.

getGeodesicArea: function(projection) {
var ring = this; // so we can work with a clone if needed
if(projection) {
var gg = new OpenLayers.Projection("EPSG:4326");
if(!gg.equals(projection)) {
ring = this.clone().transform(projection, gg);
}
}
var area = 0.0;
var len = ring.components && ring.components.length;
if(len > 2) {
var p1, p2;
for(var i=0; i<len-1; i++) {
p1 = ring.components[i];
p2 = ring.components[i+1];
area += OpenLayers.Util.rad(p2.x - p1.x) *
(2 + Math.sin(OpenLayers.Util.rad(p1.y)) +
Math.sin(OpenLayers.Util.rad(p2.y)));
}
area = area * 6378137.0 * 6378137.0 / 2.0;
}
return area;
}


i hope it helps you...

• Sadly, there's no PostGIS in the system. The polygons are computed in Python and we need the areas without a database roundtrip. – Reid Nov 14 '12 at 16:11
• check out my answer, i have updated it with a javascript code. – Aragon Nov 14 '12 at 16:16
• Can you summarize the algorithm that JavaScript code is using? For example is it doing the standard area-of-a-polygon computation with (non-ellipsoidal) geodetic distances between the points? Concise explanations of the principles are a lot more helpful to me than lengthy code to wade through.... – Reid Nov 14 '12 at 16:25
• it calculates the approximate area of the polygon were it projected onto the earth. it is referenced here and more information about Some Algorithms for Polygons on a Sphere here... – Aragon Nov 14 '12 at 16:28

The problem of computing the area of a polygon on an ellipsoid of revolution where the polygon edges are geodesics is solved in Section 6 of

C. F. F. Karney, Algorithms for geodesics, J. Geodesy (2012); DOI: http://dx.doi.org/10.1007/s00190-012-0578-z; Addenda: http://geographiclib.sf.net/geod-addenda.html

Code (C++, Python, Javascript, and Matlab) that implements this solution is available as part of GeographicLib, see

http://geographiclib.sf.net

A somewhat higher-level view of Aragon's answer is as follows. In this paper:

Robert Chamberlain & William Duquette, “Some algorithms for polygons on a sphere”. Proc. Association of American Geographers, 2007. (PDF)

there is a formula for computing area of a polygon on a sphere. Here's two different versions copied from the paper; one or the other might be easier to implement in a particular case.  $\lambda$ is latitude and $\phi$ is longitude. $R$ is Earth radius, and $A$ is the area you're trying to find.

I assume a sufficiently talented mathematician could extend this to a standard geographic ellipsoid, but spherical is good enough for me.

Another option is to simply transform to an equal-area projection and compute the geometric area.

• Accepting because this is what we ended up doing and it was accurate enough for our purposes; I believe it was much less than 1% off for the area of Colorado. – Reid Dec 8 '12 at 1:47