[GRASS-dev] Fw: Re: SV: area calculations in several GIS
Helmut Kudrnovsky
hellik at web.de
Tue Oct 2 23:12:30 PDT 2018
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
More information about the grass-dev
mailing list