Search code examples
javascriptgisgeospatialpostgiswgs84

Why exactly do I get a slightly different area compared to PostGis ST_Area?

I'm trying to automatically compute the area(s) of a building(s) as displayed on OpenStreetMaps.
In order to do that, I get the coordinates of the building polygon via the overpass API.

e.g. https://overpass-turbo.eu/ 

[out:json];
way(29858799);
out ids geom;

Which outputs me

{
  "version": 0.6,
  "generator": "Overpass API",
  "osm3s": {
    "timestamp_osm_base": "2017-10-09T09:54:02Z",
    "copyright": "The data included in this document is from www.openstreetmap.org. The data is made available under ODbL."
  },
  "elements": [

{
  "type": "way",
  "id": 29858799,
  "bounds": {
    "minlat": 47.3604067,
    "minlon": 8.5342631,
    "maxlat": 47.3612503,
    "maxlon": 8.5352457
  },
  "geometry": [
    { "lat": 47.3612503, "lon": 8.5351944 },
    { "lat": 47.3612252, "lon": 8.5342631 },
    { "lat": 47.3610145, "lon": 8.5342755 },
    { "lat": 47.3610212, "lon": 8.5345227 },
    { "lat": 47.3606405, "lon": 8.5345451 },
    { "lat": 47.3606350, "lon": 8.5343411 },
    { "lat": 47.3604067, "lon": 8.5343545 },
    { "lat": 47.3604120, "lon": 8.5345623 },
    { "lat": 47.3604308, "lon": 8.5352457 },
    { "lat": 47.3606508, "lon": 8.5352328 },
    { "lat": 47.3606413, "lon": 8.5348784 },
    { "lat": 47.3610383, "lon": 8.5348551 },
    { "lat": 47.3610477, "lon": 8.5352063 },
    { "lat": 47.3612503, "lon": 8.5351944 }
  ]
}

  ]
}

Now, before I compute the area myselfs in JavaScript, I put them into PostGIS, and see what area PostGIS gives me:

SELECT     
    ST_Area
       (
       ST_GeomFromText('
POLYGON
(
       (
             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
       )
)'
             ,4326 -- WGS84
             )
             ,false -- 
    )
    -- false: 6379.25032051953
    -- true:  6350.65051177517

Which gives 6379.25032051953 m2 for elliptic, and 6350.65051177517 for spheroid.

Now I try to compute the area in JavaScript
So I take these coordinates and put them into a JS array:

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

and try to compute the area in JavaScript, using this gis post as reference.

Math.radians = function(degrees)
{
    return degrees * Math.PI / 180.0;
};


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

But no matter which radius I choose, I am at least 2 sqm apart from the PostGis output.

More importantly, when I compute the area using the exact radius at latitude X, I come closer to the PostGIS spheric result, while when I choose the spheroid radius, I don't get a lower result as PostGis does - I get a higher result.

I am actually wondering why I get such different results.

Using google I found here http://postgis.net/docs/doxygen/2.2/dc/d52/geography__measurement_8c_a1a7c48d59bcf4ed56522ab26c142f61d.html that PostGIS ST_Area calls geography_area, and now I'm wondering which of the two result-sets is more wrong...

Is there anything wrong with this calculation ? Or does the blame go to PostGIS ?

Ironically, when I compute the area in SQL-Server, I get the PostGIS spherical area (actually 6350.65051472187)...

DECLARE @v_polygon_string varchar(1000);
DECLARE @g Geography;
SET @v_polygon_string = 'POLYGON((
47.3612503 8.5351944, 
47.3610477 8.5352063, 
47.3610383 8.5348551, 
47.3606413 8.5348784, 
47.3606508 8.5352328, 
47.3604308 8.5352457, 
47.3604120 8.5345623, 
47.3604067 8.5343545, 
47.3606350 8.5343411, 
47.3606405 8.5345451, 
47.3610212 8.5345227, 
47.3610145 8.5342755, 
47.3612252 8.5342631, 
47.3612503 8.5351944
)) ';
SET @g = Geography::STGeomFromText(@v_polygon_string,4326);
SELECT @g.STArea()

(but only if the polygon is defined using the left-hand rule [which is why the array is reversed here], otherwise I get System.ArgumentException: 24200: The specified input does not represent a valid geography instance.)

Reference:
ST_Area
Source Code

Solution

  • Ah, found the answer myselfs.
    This is because Google-Maps uses Web-Mercator projection.

    The coordinates are WGS84, and in order to calculate area from that, the coordinates need to be transformed into an area-preserving map-projection.

    Otherwise, you'll get an area whose actual area isn't necessarely preserved (depends on where on the globe), and the computed area will diverge from the actual area.

    The real underlying question is: What does ST_Area actually do ?
    And I found the answer to this question by looking into a geospacial library.

    Actually the result of 6377.28 fits with a precision of 0.02 m2 with a cylindrical equal-area sphere (6377.2695087222382) or an EckertIV-Sphere (6377.26664171461).

    While the result of 6350 fits with a (non-spheric) cylindrical-equal-area-world, or Albers-Projection.

    Details on the stackexchange-GIS-site.

    Albers-Projection: Albers-Orojection

    Cylindrical Equal-Area-Projection: Cylindrical Cylindrical 2