1

I am working on a GIS API that work on top of OpenLayers.

I have a feature that calculates the area of a drawn polygon. I have tested it on the browser and the areas I got, appear to be correct. For example, I have measured some country areas and they matched the real areas approximately.

Surely, I have to run some unit tests to prove that the area calculations are correct. I use PostGIS to generate the area value, but I get completely different results than the ones from OpenLayers.

In PostGIS, I do this:

SELECT
st_astext(st_transform(st_geomfromtext('POINT(50 0)', 4326), 900913)),
st_astext(st_transform(st_geomfromtext('POINT(-60 -5)', 4326), 900913)),
st_astext(st_transform(st_geomfromtext('POINT(-30 -30)', 4326), 900913)),
ST_Area(
    ST_GeomFromText(
        'POLYGON(
            (
                50 0,
                -60 -5, 
                -30 -30,
                50 0
            )
        )'
    ), false
) as area

I get these results:

"POINT(5565974.53966368 -7.08115455161362e-10)";
"POINT(-6679169.44759641 -557305.257274577)";
"POINT(-3339584.72379821 -3503549.84350437)";
19582793771831 // area

I introduce these three points in my test, like this:

var points = [
    new OpenLayers.Geometry.Point(5565974.53966368, -7.08115455161362e-10),
    new OpenLayers.Geometry.Point(-6679169.44759641, -557305.257274577),
    new OpenLayers.Geometry.Point(-3339584.72379821, -3503549.84350437)
];
var ring = new OpenLayers.Geometry.LinearRing(points);
var polygon = new OpenLayers.Geometry.Polygon([ring]);

Then, I used the getArea() method from OpenLayers.Measure, but I got an area of: 17050040174394.607 sq. meters

According to PostGIS, the area of this polygon should be: 19582793771831 sq. meters.

It is a big difference. what could I be doing wrong here?

Thanks

PS: I can't publish the whole code here, because it is work related.

EDIT 1 - CLOSE THE LINEAR RING

var points = [
    // Point 0 - Start point
    new OpenLayers.Geometry.Point(5565974.53966368, -7.08115455161362e-10),
    // Point 1
    new OpenLayers.Geometry.Point(-6679169.44759641, -557305.257274577),
    // Point 2
    new OpenLayers.Geometry.Point(-3339584.72379821, -3503549.84350437),
    // Point 3 - End point (the same as start point)
    new OpenLayers.Geometry.Point(5565974.53966368, -7.08115455161362e-10),
];

I tried to close the Linear Ring, but I still get: 17050040174394.607 sq. meters

3
  • 1
    Edit 1 has a different start/end point, with a Cartesian area 17100864718054.7 m², which is closer but still different. Double-check your inputs. – Mike T Sep 23 '14 at 3:09
  • @MikeT, I think the Linear Ring does not need to be closed. It closes automatically. But you said I have different start/end points. I double checked and they seem the same to me. I commented EDIT 1. Would you please explain what you mean by "different start/end points"? Thank you. – joaorodr84 Sep 23 '14 at 10:11
  • @MikeT, You are absolutely right. My mistake. I corrected the EDIT 1, but the value is still 17050040174394.607 sq. meters. – joaorodr84 Sep 23 '14 at 23:37
3

If you look at the source for getGeodesicArea, which is actually in the LinearRing class, you will will see from the comments (and the code) that this calculates an approximate area based on a sphere. If you look at GeographicLib, which includes a Javascript port of the C++ library, you will see that they use ellipsoidal calculations, which will be more accurate, and as Mike T has already pointed out, gives the answer of 19518154994956.3 m², as does Postgis. So, basically, the Postgis calculation is considerably more accurate than the OpenLayers one. To be clear, OpenLayers is a fabulous web mapping library, but not intended for high precision GIS calculations (as noted by the developers in the comments).

EDIT. You can use the geographiclib Javascript directly in your own code, as it is under MIT license, see License. So, if you wish to do this client side, and get an answer that agrees with Postgis, I would recommend using this library instead of OpenLayers for the geodesic area part. The documentation gives an example of using the Javascript functions -- thanks to Mike T for the link.

It is also worth pointing out that the OpenLayers calculation is based on 4326, not 900913. Google developed used Spherical Mercator when they introduced slippy maps, for ease of calculation, not for accurate geodesic calculations. There is a chronic distortion of scale, the further you get from the equator, so I would not trust any area calculations based on 900913, if you were thinking of using it, which is not entirely clear from your code.

5
  • @JonhBarça, First of all, thank you for your answer. I use a plygon handler to draw the sketch, and while I am drawing, the current area is calculated and shown. I guess I could use directly 4236 coordinates, but OpenLayer does the conversion automatically, I think. Do you suggest calculating the area "without" OpenLayers, then? You mentioned GeographicLib, but the documentation is quite poor and I don't know how to use it. :s – joaorodr84 Sep 23 '14 at 15:58
  • 1
    If you do view source on the geographiclib page you will see the link to the js source file, geographiclib.sourceforge.net/scripts/geographiclib.js. It says it is copyrighted, but issued under an MIT license. I am no expert on licensing, so maybe email them and ask if it is OK to use it. If you want more accurate area measures on the client than Vincenty, then use this library (assuming you can). – John Powell Sep 23 '14 at 16:02
  • The documentation for GeographicLib is far from "poor". Look again. Furthermore, here is an example. – Mike T Sep 23 '14 at 23:42
  • OK, so you can use it freely. Answer updated. – John Powell Sep 24 '14 at 6:17
  • @MikeT. Agreed, though the Javascript is a lot harder to read than that of OpenLayers, it is more accurate, which is the main thing here. – John Powell Sep 24 '14 at 6:18
3

The geodesic area for your example is 19518154994956.3 m² (using GeographicLib). E.g.:

var points = [
    {lat: 0, lon: 50},
    {lat: -5, lon: -60},
    {lat: -30, lon: -30}
];
var p = GeographicLib.PolygonArea;
var result = p.Area(GeographicLib.Geodesic.WGS84, points, false);
var area = result.area; // 19518154994956.285

Here are multiple interpretations of area in PostGIS:

  • ST_Area(geog::geography, false) = 19582793771831 m² is a calculation on a sphere, and is only approximate.
  • ST_Area(geog::geography, true) for spheroid/geodesic calculation raises an ERROR for this example due to a limitation that will hopefully be fixed with PostGIS 2.2.
  • ST_Area(ST_Transform(geom, 900913)) = 18969178651878 m² is a flat-word Cartesian area on a spherical Mercator projection.

I don't know where the values of getArea() come from with your example, but I would have expected it to be 18969178651878 m². Does it make any difference if you close the linearring?


Looking at the source code for OpenLayers, it looks like getGeodesicArea is approximating the area on a sphere, and not a geodesic of revolution (or spheroid). Furthermore, this function expects coordinate units of degrees. It's not clear what you are doing with the spherical Mercator projected points in the examples in the question. As there are multiple definitions of area, which is regarded by you to be "correct"? The options: are geodesic, spherical or projected.

5
  • It didn't make any difference when I closed the linear ring, as you can see in edit 1. I've got to dig deeper in this. I'll post the advances I make here. I'll stay in touch. Thanks – joaorodr84 Sep 23 '14 at 1:34
  • 1
    The area function is OpenLayers uses the standard approach for calculating an area on planar coordinates, see, getArea in trac.osgeo.org/openlayers/browser/trunk/openlayers/lib/…. and does not consider the projection at all. As a general rule you should always close polygons -- it is expected for valid GeoJSON, WKT, WKB, etc. I would not use 900913 for taking areas -- it is a projection that is designed for ease of use in slippy maps and suffers from chronic scale distortion, so that area calculations will be very unreliable. – John Powell Sep 23 '14 at 10:21
  • @JohnBarça, I use this getArea from Measure.js, which uses the getGeodesicArea from Geometry/Collection.js. It seems to return true geodesic areas. Ijust wanted to know if OpenLayers is wrong or if I am using PostGIS wrong in my tests. – joaorodr84 Sep 23 '14 at 11:21
  • Apologies, I misread your question. OK, I looked at the source for getGeodesicArea and have answered your question. – John Powell Sep 23 '14 at 14:18
  • @MikeT, answering your edit. I use the points to create a Linear Ring and then a Polygon, as you can see. Using the polygon I have an OpenLayers.Measure object, which used the polygon in the call of measureControl.getArea(polygon). The method getArea transforms the projection into 4326, and then ccalls GetGeodesicArea. So, I think the projection might not be the problem. – joaorodr84 Sep 23 '14 at 23:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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