How does ST_Area in PostGIS work?

Question

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?

Stefan Steiger · Accepted Answer · 2019-02-20 22:54:48Z

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 m² to km², you need to divide m² by 1000² (1000 meters in a kilometer - it squares because it's an area, so 1000*1000 aka 1000²).

ST_Area(geometry, true/false) calculates the area (in m²) 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 2m², 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;
        }


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

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

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
                    + Math.sin(Math.radians(p1[1]))
                    + Math.sin(Math.radians(p2[1]))
                );
        }

        // https://en.wikipedia.org/wiki/Earth_radius#Equatorial_radius
        // 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 φ
        var siteLatitude = Math.radians(poly[0][0]);


        // 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);

        // https://en.wikipedia.org/wiki/Earth_radius#Radius_at_a_given_geodetic_latitude
        // Geocentric radius
        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);
        // var hrad = polar_radius + (90 - Math.abs(siteLatitude)) / 90 * (equatorial_radius - polar_radius);
        var radius = mean_radius;

        area = area * radius * radius / 2.0;
    } // End if len > 0

    // equatorial_radius: 6391.565558418869 m2
    // mean_radius:       6377.287126172337m2
    // authalic_radius:   6377.283923019292 m2
    // volumetric_radius: 6377.271110415153 m2
    // merid_radius:      6375.314923754325 m2
    // polar_radius:      6348.777989748668 m2
    // R:                 6368.48180842528 m2
    // hrad:              6391.171919886588 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

and here:
http://postgis.net/docs/doxygen/2.2/d1/dc0/lwspheroid_8c_a29d141c632f6b46587dec3a1dbe3d176.html#a29d141c632f6b46587dec3a1dbe3d176

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"?.. — Marco_B, Feb 19, 2019 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? — Marco_B, Feb 19, 2019 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/… — Stefan Steiger, Feb 20, 2019 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. — Marco_B, Feb 20, 2019 at 10:48
@Marco_B: I had already started doing that when you wrote that comment. Now done. — Stefan Steiger, Feb 20, 2019 at 23:01

yxcv · Accepted Answer · 2015-11-08 06:36:07Z

6

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.

edited Nov 8, 2015 at 6:36

answered Nov 7, 2015 at 19:42

yxcv

4091 gold badge3 silver badges6 bronze badges

This is correct. For reference: postgis.net/docs/ST_Area.html
– alphabetasoup
Nov 7, 2015 at 19:59
But what is the basis for dividing by 100,000? km^2 = m^2 / 1,000,000
– J Tileson
Nov 7, 2015 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.
– Vince
Nov 7, 2015 at 23:59
@Vince Your are right. I corrected the typo.
– yxcv
Nov 8, 2015 at 6:36
@ J Tilesone The question of the basis for that division you can ask on math.stackexchange.com.
– yxcv
Nov 8, 2015 at 6:39

Add a comment |

Stack Exchange Network

How does ST_Area in PostGIS work?

2 Answers 2

Your Answer

Not the answer you're looking for? Browse other questions tagged
postgis
area
or ask your own question.

Linked

Hot Network Questions

How does ST_Area in PostGIS work?

2 Answers 2

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged postgisarea or ask your own question.

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
postgis
area
or ask your own question.