[GRASS-dev] Fw: Re: SV: area calculations in several GIS

[Helmut Kudrnovsky]

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 at 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