# How does ST_Area in PostGIS work?

I am performing a simple calculation on a polygon known to be approximately 6226 km^2 in area. It is stored in a Geography (WGS84 SRID) column.

The query is:

``````select st_astext(col), st_area(col) area from table
``````

and returns:

``````"POLYGON((-180 58.282525588539,-178.916399160189 57.4759784390599,-178.191728834624 58.5761461944577,-180 58.282525588539))" | 5807028547.33813
``````

The area returned (5807028547.33813) appears to be mm^2 and not km^2? The documentation http://postgis.net/docs/ST_Area.html states "by default area is determined on a spheroid with units in square meters"

Is this a documentation error, or is the above correct and I'm fundamentally misunderstanding the functionality?

SELECT
st_astext(col)
,st_area(col, false) AS area
FROM table

ST_Area(geometry) calculates the polygon area as WGS1984, WITHOUT projecting to equal area sphere/ellipsis (if you use the sql-type Geometry instead of Geography). The result is measured in the unit in the geometry's SRID.

ST_Area(geography) calculates the polygon area as WGS1984, WITH projecting to equal area sphere/ellipsis (if you use the sql-type Geography instead of Geometry). The result is measured in square meters. To get from m2 to km2, you need to divide m2 by 10002 (1000 meters in a kilometer - it squares because it's an area, so 1000*1000 aka 10002).

ST_Area(geometry, true/false) calculates the area (in m2) with coordinates projected into CylindricalEqualAreaworld coordinate system (preserving area - makes sense if you want to calculate the area).

The difference between true/false is the accuracy.
ST_Area(geog,false) uses a faster but less accurate sphere.

Say, when I use this polygon:

``````var poly = [
[47.3612503, 8.5351944],
[47.3612252, 8.5342631],
[47.3610145, 8.5342755],
[47.3610212, 8.5345227],
[47.3606405, 8.5345451],
[47.3606350, 8.5343411],
[47.3604067, 8.5343545],
[47.3604120, 8.5345623],
[47.3604308, 8.5352457],
[47.3606508, 8.5352328],
[47.3606413, 8.5348784],
[47.3610383, 8.5348551],
[47.3610477, 8.5352063],
[47.3612503, 8.5351944]
];
``````

I get the following results:

``````ST_Area(g) =             5.21556075001092E-07
ST_Area(g, false)     6379.25032051953
ST_Area(g, true)      6350.65051177517
``````

I think the important part to be taken from the docs is this:

For geometry, a 2D Cartesian area is determined with units specified by the SRID.
For geography, by default area is determined on a spheroid with units in square meters.

So you need to be careful to choose geography, and NOT geometry.
If you use geometry, you NEED to use the true/false overloads of ST_Area.

In C#, I get more-or-less the same as true with KnownCoordinateSystems.Projected.World.CylindricalEqualAreaworld, and false seems to be an earth-mean-radius-world, something close to WorldSpheroid.CylindricalEqualAreasphere or WorldSpheroid.EckertIVsphere, but it's off by 2m2, so it seems to do its own thing.

``````using DotSpatial.Projections;
using DotSpatial.Topology;

namespace TestSpatial
{

static class Program
{

// https://stackoverflow.com/questions/46159499/calculate-area-of-polygon-having-wgs-coordinates-using-dotspatial
// pfff wrong...
public static void TestPolygonArea()
{
// this feature can be see visually here http://www.allhx.ca/on/toronto/westmount-park-road/25/
string feature = "-79.525542519049552,43.691278124243432 -79.525382520578987,43.691281097414787 -79.525228855617627,43.69124858593392 -79.525096151437353,43.691183664769774 -79.52472799258571,43.690927163079735 -79.525379447437814,43.690771996666641 -79.525602330675355,43.691267524226838 -79.525542519049552,43.691278124243432";
feature = "47.3612503,8.5351944 47.3612252,8.5342631 47.3610145,8.5342755 47.3610212,8.5345227 47.3606405,8.5345451 47.3606350,8.5343411 47.3604067,8.5343545 47.3604120,8.5345623 47.3604308,8.5352457 47.3606508,8.5352328 47.3606413,8.5348784 47.3610383,8.5348551 47.3610477,8.5352063 47.3612503,8.5351944";

string[] coordinates = feature.Split(' ');
// System.Array.Reverse(coordinates);

// dotspatial takes the x,y in a single array, and z in a separate array.  I'm sure there's a
// reason for this, but I don't know what it is.'
double[] xy = new double[coordinates.Length * 2];
double[] z = new double[coordinates.Length];
for (int i = 0; i < coordinates.Length; i++)
{
double lon = double.Parse(coordinates[i].Split(',')[0]);
double lat = double.Parse(coordinates[i].Split(',')[1]);
xy[i * 2] = lon;
xy[i * 2 + 1] = lat;
z[i] = 0;
}

double area = CalculateArea(xy);
System.Console.WriteLine(area);
}

public static double CalculateArea(double[] latLonPoints)
{
// source projection is WGS1984
ProjectionInfo projFrom = KnownCoordinateSystems.Geographic.World.WGS1984;
// most complicated problem - you have to find most suitable projection
ProjectionInfo projTo = KnownCoordinateSystems.Projected.UtmWgs1984.WGS1984UTMZone37N;
projTo = KnownCoordinateSystems.Projected.Europe.EuropeAlbersEqualAreaConic; // 6350.9772005155683
// projTo= KnownCoordinateSystems.Geographic.World.WGS1984; // 5.215560750019806E-07
projTo = KnownCoordinateSystems.Projected.WorldSpheroid.EckertIVsphere; // 6377.26664171461
projTo = KnownCoordinateSystems.Projected.World.EckertIVworld; // 6391.5626849671826
projTo = KnownCoordinateSystems.Projected.World.CylindricalEqualAreaworld; // 6350.6506013739854
projTo = KnownCoordinateSystems.Projected.WorldSpheroid.CylindricalEqualAreasphere; // 6377.2695087222382
projTo = KnownCoordinateSystems.Projected.WorldSpheroid.EquidistantCylindricalsphere; // 6448.6818862780929
projTo = KnownCoordinateSystems.Projected.World.Polyconicworld; // 8483.7701716953889
projTo = KnownCoordinateSystems.Projected.World.EquidistantCylindricalworld; // 6463.1380225215107
projTo = KnownCoordinateSystems.Projected.World.EquidistantConicworld; // 8197.4427198320627
projTo = KnownCoordinateSystems.Projected.World.VanderGrintenIworld; // 6537.3942984174937
projTo = KnownCoordinateSystems.Projected.World.WebMercator; // 6535.5119516421109
projTo = KnownCoordinateSystems.Projected.World.Mercatorworld; // 6492.7180733950809
projTo = KnownCoordinateSystems.Projected.SpheroidBased.Lambert2; // 9422.0631835013628
projTo = KnownCoordinateSystems.Projected.SpheroidBased.Lambert2Wide; // 9422.0614012926817
projTo = KnownCoordinateSystems.Projected.TransverseMercator.WGS1984lo33; // 6760.01638841012
projTo = KnownCoordinateSystems.Projected.Europe.EuropeAlbersEqualAreaConic; // 6350.9772005155683
projTo = KnownCoordinateSystems.Projected.UtmOther.EuropeanDatum1950UTMZone37N; // 6480.7883094931021

// ST_Area(g, false)     6379.25032051953
// ST_Area(g, true)      6350.65051177517
// ST_Area(g)            5.21556075001092E-07

// prepare for ReprojectPoints (it's mutate array)
double[] z = new double[latLonPoints.Length / 2];
// double[] pointsArray = latLonPoints.ToArray();

Reproject.ReprojectPoints(latLonPoints, z, projFrom, projTo, 0, latLonPoints.Length / 2);

// assemblying new points array to create polygon
System.Collections.Generic.List<Coordinate> points =
new System.Collections.Generic.List<Coordinate>(latLonPoints.Length / 2);

for (int i = 0; i < latLonPoints.Length / 2; i++)
points.Add(new Coordinate(latLonPoints[i * 2], latLonPoints[i * 2 + 1]));

Polygon poly = new Polygon(points);
return poly.Area;
}

static void Main(string[] args)
{
TestPolygonArea();

System.Console.WriteLine(System.Environment.NewLine);
System.Console.WriteLine(" --- Press any key to continue --- ");
}
}
}
``````

e.g. you get a close-fit to false with the mean-radius:

``````// https://gis.stackexchange.com/a/816/3997
function polygonArea()
{
var poly = [
[47.3612503, 8.5351944],
[47.3612252, 8.5342631],
[47.3610145, 8.5342755],
[47.3610212, 8.5345227],
[47.3606405, 8.5345451],
[47.3606350, 8.5343411],
[47.3604067, 8.5343545],
[47.3604120, 8.5345623],
[47.3604308, 8.5352457],
[47.3606508, 8.5352328],
[47.3606413, 8.5348784],
[47.3610383, 8.5348551],
[47.3610477, 8.5352063],
[47.3612503, 8.5351944]
];

var area = 0.0;
var len = poly.length;

if (len > 2)
{

var p1, p2;

for (var i = 0; i < len - 1; i++)
{

p1 = poly[i];
p2 = poly[i + 1];

area += Math.radians(p2[0] - p1[0]) *
(
2
);
}

// https://en.wikipedia.org/wiki/Earth_ellipsoid
// The radius you are using, 6378137.0 m corresponds to the equatorial radius of the Earth.
var equatorial_radius = 6378137; // m
var polar_radius = 6356752.3142; // m
var mean_radius = 6371008.8; // m
var authalic_radius = 6371007.2; // m (radius of perfect sphere with same surface as reference ellipsoid)
var volumetric_radius = 6371000.8 // m (radius of a sphere of volume equal to the ellipsoid)
// geodetic latitude φ

// https://en.wikipedia.org/wiki/Semi-major_and_semi-minor_axes
// https://en.wikipedia.org/wiki/World_Geodetic_System
var a = 6378137; // m
var b = 6356752.3142; // m
// where a and b are, respectively, the equatorial radius and the polar radius.

var R1 = Math.pow(a * a * Math.cos(siteLatitude), 2) + Math.pow(b * b * Math.sin(siteLatitude), 2)
var R2 = Math.pow(a * Math.cos(siteLatitude), 2) + Math.pow(b * Math.sin(siteLatitude), 2);

var R = Math.sqrt(R1 / R2);
// var merid_radius = ((a * a) * (b * b)) / Math.pow(Math.pow(a * Math.cos(siteLatitude), 2) + Math.pow(b * Math.sin(siteLatitude), 2), 3/2)

// console.log(R);

} // End if len > 0

// R:                 6368.48180842528 m2

// http://postgis.net/docs/doxygen/2.2/dc/d52/geography__measurement_8c_a1a7c48d59bcf4ed56522ab26c142f61d.html
// ST_Area(false)     6379.25032051953
// ST_Area(true)      6350.65051177517

// return area;
return area.toFixed(2);
}
``````

WebMercator is the coordinate system used by Google-Maps.
The official name for this coordinate system is EPSG:3857.

What exactly PostGIS does, is documented here:
https://postgis.net/docs/ST_Area.html

And details in source-code can be found here:
http://postgis.net/docs/doxygen/2.2/dc/d52/geography__measurement_8c_a1a7c48d59bcf4ed56522ab26c142f61d.html

Albers-Projection:

Cylindrical Equal-Area-Projection:

• Stefan, this is weird. Your remark and actual calculation results for the "ST_Area(geography)" case, so without specifying the "use_spheroid" boolean, seem to contradict the official PostGIS documentation. In the documentation, it says "For geography, by default area is determined on a spheroid with units in square meters.", so without specifying the boolean, it should default to giving area in m2 based on spheroid. This is also the result shown in the last example in the grey box on the Help page, where clearly no boolean is used in the ST_Area command, yet the result shows "sqm_spheroid"?.. Feb 19 '19 at 12:34
• Stefan, was your source data actually in geography, or converted to geography, or just geometry happening to be in WGS1984 lat/long? Feb 19 '19 at 12:37
• @Marco_B: Definitely Geometry. Otherwise, the 5.21556075001092E-07 cannot be explained. I see, i should have used geography instead of geometry. see github.com/ststeiger/AnySqlWebAdmin/blob/master/AnySqlWebAdmin/… Feb 20 '19 at 7:37
• Thanks for the response, that explains it then. Maybe it is good to edit your initial post, to change the now confusing and somewhat inaccurate statement of "ST_Area(geography) calculates the polygon area as WGS1984, without projecting to equal area sphere/ellipsis.", and that this statement only applies to the situation of the shapes being in Geometry, not Geography storage, and that in the case of Geography, the shapes will be using a spheroid / sphere. Feb 20 '19 at 10:48
• @Marco_B: I had already started doing that when you wrote that comment. Now done. Feb 20 '19 at 23:01

The calculation gives the right output. As you cited the documentation states "by default area is determined on a spheroid with units in square meters".

The result of your query is 5807028547.33813 m^2. To get the area in km^2 you have to devide the result by 1,000,000.

`5807028547.33813 m^2 / 1,000,000 = 5807.02854733813 km^2`

5807.02854733813 km^2 corresponds approximately to your expected 6226 km^2.

• This is correct. For reference: postgis.net/docs/ST_Area.html Nov 7 '15 at 19:59
• But what is the basis for dividing by 100,000? km^2 = m^2 / 1,000,000 Nov 7 '15 at 21:03
• It's a typo. The value should be 1.0e+006, but the result was right (5807 km^2) mm^2 would have been 1.0e+006 larger still. Nov 7 '15 at 23:59
• @Vince Your are right. I corrected the typo.
– yxcv
Nov 8 '15 at 6:36
• @ J Tilesone The question of the basis for that division you can ask on math.stackexchange.com.
– yxcv
Nov 8 '15 at 6:39