How to compute geodesic area in GeoDjango?

Question

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.

Answer 1

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[2] # 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...

Answer 2

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

Answer 3

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.

Answer 4

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