Saltar al contenido

¿Cómo funciona ST_Area en PostGIS?

Nuestro equipo especializado pasados algunos días de investigación y de recopilar de información, dieron con la respuesta, queremos que resulte de gran utilidad en tu plan.

Solución:

SELECCIONE
st_astext (col)
, st_area (col, falso) AS área
DESDE la mesa

ST_Area (geometría) calcula el área del polígono como WGS1984, SIN proyectando a un área igual a la esfera / puntos suspensivos (si usa la geometría de tipo sql en lugar de la geografía). El resultado se mide en la unidad en el SRID de la geometría.

ST_Area (geografía) calcula el área del polígono como WGS1984, CON proyectando a un área igual a la esfera / puntos suspensivos (si usa Geografía de tipo sql en lugar de Geometría). El resultado se mide en metros cuadrados. Para obtener de m2 a km2, necesitas dividir m2 por 10002 (1000 metros en un kilómetro: cuadra porque es un área, por lo que 1000 * 1000 también conocido como 10002).

ST_Area (geometría, verdadero / falso) calcula el área (en m2) con coordenadas proyectadas en el sistema de coordenadas CylindricalEqualAreaworld (preservar el área, tiene sentido si desea calcular el área).

La diferencia entre verdadero / falso es la precisión.
ST_Area (geog, false) utiliza una esfera más rápida pero menos precisa.

Diga, cuando uso este 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]
];

Obtengo los siguientes resultados:

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

Creo que la parte importante que se debe tomar de los documentos es esta:

Para geometría, un área cartesiana 2D se determina con unidades especificadas por el SRID.
Para geografía, por defecto el área se determina en un esferoide con unidades en metros cuadrados.

Así que debes tener cuidado al elegir geografía, y NO geometría.
Si usa geometría, NECESITA usar las sobrecargas verdadero / falso de ST_Area.

En C #, obtengo más o menos lo mismo que verdadero con KnownCoordinateSystems.Projected.World.CylindricalEqualAreaworld, y false parece ser un mundo de radio medio terrestre, algo cercano a WorldSpheroid.CylindricalEqualAreasphere o WorldSpheroid.EckertIVsphere, pero es de 2 m2, por lo que parece hacer lo suyo.

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 points = 
                new System.Collections.Generic.List(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 ejemplo, obtiene un ajuste cercano a falso con el radio medio:

// 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 es el sistema de coordenadas utilizado por Google-Maps.
El nombre oficial de este sistema de coordenadas es EPSG: 3857.

Lo que hace exactamente PostGIS se documenta aquí:
https://postgis.net/docs/ST_Area.html

Y los detalles en el código fuente se pueden encontrar aquí:
http://postgis.net/docs/doxygen/2.2/dc/d52/geography__measurement_8c_a1a7c48d59bcf4ed56522ab26c142f61d.html

y aquí:

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

Proyección de Albers:
Proyección de Albers

Proyección de Albers 2

Proyección cilíndrica de igual área:
Cilíndrico Cilíndrico 2

El cálculo da la salida correcta. Como citó, la documentación indica que "el área predeterminada se determina en un esferoide con unidades en metros cuadrados".

El resultado de su consulta es 5807028547.33813 m ^ 2. Para obtener el área en km ^ 2, debe dividir el resultado por 1,000,000.

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

5807.02854733813 km ^ 2 corresponde aproximadamente a los 6226 km ^ 2 esperados.

¡Haz clic para puntuar esta entrada!
(Votos: 0 Promedio: 0)



Utiliza Nuestro Buscador

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *