Como o ST_Area no PostGIS funciona?

Estou realizando um cálculo simples em um polígono com aproximadamente 6226 km ^ 2 de área. Ele é armazenado em uma coluna Geografia (WGS84 SRID).

A consulta é:

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

e retorna:

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

A área retornada (5807028547.33813) parece ser mm ^ 2 e não km ^ 2? A documentação http://postgis.net/docs/ST_Area.html declara "por padrão a área é determinada em um esferóide com unidades em metros quadrados"

Isso é um erro de documentação ou o que foi descrito acima está correto e estou basicamente entendendo mal a funcionalidade?

O cálculo fornece a saída correta. Como você citou, a documentação indica "a área padrão é determinada em um esferóide com unidades em metros quadrados".

O resultado da sua consulta é 5807028547.33813 m ^ 2 . Para obter a área em km ^ 2, você deve dividir o resultado em 1.000.000.

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

5807.02854733813 km ^ 2 corresponde aproximadamente aos 6226 km ^ 2 esperados.

SELECT área
     st_astext (col)
    , st_area (col, false ) AS area
FROM tabela

ST_Area (geometry) calcula a área do polígono como WGS1984, SEM projetar para área / elipse de área igual (se você usar a geometria do tipo sql em vez de geography). O resultado é medido na unidade no SRID da geometria.

ST_Area (geography) calcula a área do polígono como WGS1984, COM projetando-se para área / elipse de área igual (se você usar a geografia do tipo sql em vez de geometria). O resultado é medido em metros quadrados. Para ir de m ² a km ² , é necessário dividir m ² por 1000 ² (1000 metros em um quilômetro - o quadrado é porque é uma área, então 1000 * 1000 aka 1000 ² ).

ST_Area (geometria, verdadeiro / falso) calcula a área (em m ² ) com coordenadas projetadas no sistema de coordenadas CylindricalEqualAreaworld (preservar área - faz sentido se você deseja calcular a área).

A diferença entre verdadeiro / falso é a precisão.
ST_Area (geog, false) usa uma esfera mais rápida, mas menos precisa.

Digamos, quando eu uso esse polígono:

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]
];

Eu recebo os seguintes resultados:

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

Eu acho que a parte importante a ser tirada dos documentos é a seguinte:

Para geometria , uma área cartesiana 2D é determinada com unidades especificadas pelo SRID.
Para geo grafia , por área padrão é determinado com um esferóide com unidades em metros quadrados .

Portanto, você precisa ter cuidado ao escolher a geografia , e NÃO a geometria.
Se você usa geometria, PRECISA usar sobrecargas verdadeiras / falsas de ST_Area.

Em C #, fico mais ou menos o mesmo que verdadeiro com KnownCoordinateSystems.Projected.World.CylindricalEqualAreaworld, e false parece ser um mundo com raio médio da terra, algo próximo ao WorldSpheroid.CylindricalEqualAreasphere ou WorldSpheroid.EckertIVsphere, mas é 2m ² , então parece fazer o que quiser.

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();
        }
    }
}

por exemplo, você obtém um ajuste perfeito para false com o raio médio:

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

O WebMercator é o sistema de coordenadas usado pelo Google-Maps.
O nome oficial para este sistema de coordenadas é EPSG: 3857.

O que exatamente o PostGIS faz, está documentado aqui:
https://postgis.net/docs/ST_Area.html

E detalhes no código-fonte podem ser encontrados aqui:
http://postgis.net/docs/doxygen/2.2/dc/d52/geography__measurement_8c_a1a7c48d59bcf4ed56522ab26c142f61d.html

e aqui:
http://postgis.net/docs/doxygen/2.2/d1/dc0/lwspheroid_8c_a29d141c632f6b46587dec3a1dbe3d176.html#a29d141c632f6b46587dec3a1dbe3d176

Albers-Projeção:

Projeção cilíndrica de área igual: