for the record to the ML, see below
---------------------------------------------------------------------------------------------------------
Gesendet: Dienstag, 02. Oktober 2018 um 22:45 Uhr
Von: "Charles Karney"
An: "Markus Metz"
Cc: "Kristian Evers", "GRASS developers list" <grass-dev@lists.osgeo.org>, "Helmut Kudrnovsky"
Betreff: Re: SV: [GRASS-dev] area calculations in several GIS
In comparisons like this, it's probably a good idea to document what
areas the various packages are attempting to compute.
For Planimeter/geod_polygonarea, it is the area of a polygon with edges
given by geodesics.
For GRASS/QGIS, it is the area of a polygon with edges which are
straight lines in the plate-carree projection.
The documentation for GRASS (in lib/gis/area_poly1.c) suggests that this
calculation is a poor man's attempt to compute the area of a polygon
with edges given by rhumb lines.
In fact, the mathematics for computing the area of rhumb polygons is
straightforward, see
https://geographiclib.sourceforge.io/html/rhumb.html#rhumbarea
and GeographicLib's Planimeter utility can perform such calculations
(with the -R flag). For Kristian's sample polygon we have
geodesic area = 14737935340.10 m^2
rhumb area = 14722522188.60 m^2
This explains the apparent discrepancy between Planimeter and
geod_polygonarea reported earlier. The figure given for Planimeter was
the rhumb area.
I recommend that GRASS/QGIS be updated either to compute the area of
either geodesic or rhumb polygons. Surely plate-carree polygons are not
useful?
--Charles
On 10/2/18 1:55 PM, Markus Metz wrote:
On Tue, Oct 2, 2018 at 7:30 PM Charles Karney wrote:
>
> Sorry, I'm coming late to this conversation. The area of the polygon
> in the Baltic Sea posted by Kristian is
>
> 14737935340.098511 m^2
>
> assuming the WGS84 ellipsoid. This is accurate to 1 square-mm and was
> obtained with the MPFR-enabled version of GeographicLib's Planimeter.
>
> This is consistent with the result reported by the PROJ.4
> geod_polygonarea() and it agrees with the result I get with the regular
> Planimeter utility. I don't why someone else is getting a different
> result from Planimeter.
I can confirm this, the results of Planimeter and PROJ 5.2.0 are identical.
Some more confusion:
I created simple boxes for the test polygon, one a bit larger, one a bit
smaller than the test polygon
GRASS: 13,222.778
Planimeter: 13,221.965
geod_polygonarea(): 13,221.965
14.62569 55.36254
14.6256944444444 54.0786
17.1177 54.0786111111111
17.1177777777778 55.3625
GRASS: 22,950.510
Planimeter: 22,946.901
geod_polygonarea(): 22,946.901
GRASS native geodesic area is not too far off from GeographicLib/PROJ
With the test polygon
Planimeter/PROJ: 14,737.935 km^2
GRASS GIS: 14.718.098 km^2
GRASS is much farther off, and here smaller instead of larger than
Planimeter/PROJ.
I will have a look.
Markus M
> It's probably operator error -- but who knows?
>
> --Charles