If I have polygons in geographic coordinates (WGS84), how do I measure the total area each takes on the surface of the earth, taking into account the curvature of the earth?
5 Answers
Here's a link to some code that'll yield the area of a simple polygon (originally from the World Wind Forum): http://forum.worldwindcentral.com/showthread.php?t=20724. This solves the problem on a sphere, roughly based on the relationship:
S = area of polygon; theta is the sum of interior angles in radians; n is the number of vertices; r is the radius of the sphere.
See also (source of formula image): http://www.geom.uiuc.edu/docs/reference/CRC-formulas/node59.html
I would be delighted to see links and/or code for polygon area on an oblate spheroid.
-
1@glennon Code for area of a geodesic polygon on a spheroid is available in GeographicLib. The code contains a link to the paper where the algorithm is derived.– cffkCommented Aug 7, 2011 at 20:21
-
1This method, based on the Gauss-Bonnet formula, is of great mathematical interest but is almost worthless for small figures on the earth: there is enormous cancellation to find the minuscule "excess" within parentheses, often leading to catastrophic loss of precision--it's easy to lose almost all precision entirely, even when working in double precision arithmetic. For routine GIS calculations, it just doesn't work.– whuberCommented Apr 22, 2013 at 22:15
Here's the source for the simplified calculation that we make in OpenLayers. This method comes from "Some Algorithms for Polygons on a Sphere" (Robert. G. Chamberlain and William H. Duquette, NASA JPL Publication 07-03). The code linked to above is for determining the area of a linear ring (with geographic coordinates). Areas for Polygons and MultiPolygons are summed up from the rings.
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;
}
Ring components are two element arrays of x, y (lon, lat) coords in the above code. The OpenLayers.Util.rad method just converts degrees to radians (deg * PI / 180).
-
1The radius you are using, 6378137.0 m corresponds to the equatorial radius of the Earth. Is there a reason to that? Wouldn't using a medium radius be more precise ?– FredBCommented Aug 29, 2013 at 15:39
-
Math.radians = function (degrees) { return degrees * Math.PI / 180.0; }; x = latitude, y = longitude Commented Jul 18, 2018 at 14:28
-
@StefanSteiger I'm pretty sure x is longitude but even so it varies slightly from google map's computeArea Commented Sep 30, 2020 at 9:29
-
@Anton Duzenko: Yea you're right, I discovered that, too - some time ago. I was always wondering why one would do long/lat, coordinates instead of lat/long, but then I realized my mistake. Understanding lat/long is a silly/tricky thing. The funny thing is, as long as you just determine if a point is inside a polygon,it doesn't seem to matter,as long as you consequently do it the wrong way - perfect symmetry ;) Commented Sep 30, 2020 at 9:52
PostGIS 1.5 introduced a new GEOGRAPHY type. The GEOGRAPHY
type allows for unprojected coordinates on a spheroid to be stored in a PostGIS table, and some analysis functions to be performed upon them.
ST_Area
queries can be performed upon GEOGRAPHY type polygons in order to calculate their area in square meters.
The following query outputs the area of all polygons using the spheroid (currently only the WGS-84
spheroid is supported), assuming they are stored using the GEOGRAPHY
type:
SELECT ST_Area(the_geom) FROM table_of_polygons;
The algorithm used to calculate area on a spheroid can be derived from the source-code.
You'll need to convert your geographic coordinates into a projection that has a coordinate system that allows you to use Cartesian math to calculate area.
I believe UTM is the accepted standard projection, as it is very simple to select a zone based on your latitude and longitude, and also the distortion is minimal, even across zones. So, if you have a polygon the size of Texas, you can use UTM Zone 14 N and it will still be fairly accurate.
If your polygons are over the North or South pole, then you should use UPS instead, as the UTM projections are less accurate over the poles, and you will quickly traverse them as the boundaries get smaller (since they follow lines of longitude)
Once your points are in a Cartesian-friendly coordinate system, you can treat them like polygons on a grid and calculate area.
A PolygonArea class was added to GeographicLib in 2011-07. This computes the true ellipsoidal area of a polygon whose edges are geodesics. Unlike PostGIS, The method does not entail numerical integration. For documenation (and a link to the paper where the formulas are derived), see
http://geographiclib.sf.net/html/classGeographicLib_1_1PolygonAreaT.html
(Link fixed to reflect generalization of PolygonArea to a template class.)