[GRASS-dev] area calculations in several GIS
Markus Metz
markus.metz.giswork at gmail.com
Thu Oct 4 12:15:56 PDT 2018
On Wed, Oct 3, 2018 at 8:58 AM Kristian Evers <kreve at sdfe.dk> wrote:
>
> Thanks for clearing up the confusion, Charles! Looking back in the thread
> it seem to have been me who supplied the rhumb line calculation. I guess
> that comes down to a cut'n'paste mishap. I'm glad we got to the bottom
> of this finally.
>
> I'll also recommend that GRASS and QGIS switch to using either geodesic
> or rhumb line area determination. Or even better, offer the ability to
> choose various method of calculating areas. Several of the GIS
applications
> I ran my test polygon through had such features and it did give me a bit
> more confidence in those calculations since it was clear what it was doing
> (albeit not as accurately as I would like in most cases).
GRASS is already using PROJ, thus it would be easy for GRASS to make use of
geodesic.[h|c] as included in PROJ. IMHO, we should not reinvent the wheel
for GRASS.
Out of curiosity, I found that the polygon size does not matter much with
simple polygons. The number of vertices has a much larger influence:
results diverge with an increasing number of vertices, results of
geod_polygonarea() become larger than GRASS native area size which stays
relatively constant. Is there an easy explanation for this divergence
increasing with the number of vertices?
Markus M
>
> /Kristian
>
> -----Oprindelig meddelelse-----
> Fra: Charles Karney <charles.karney at gmail.com> På vegne af Charles Karney
> Sendt: 2. oktober 2018 22:46
> Til: Markus Metz <markus.metz.giswork at gmail.com>
> Cc: Kristian Evers <kreve at sdfe.dk>; GRASS developers list <
grass-dev at lists.osgeo.org>; Helmut Kudrnovsky <hellik at web.de>
> Emne: 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 <charles at karney.com
> > <mailto:charles at karney.com>> 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
> > >
> > > On 10/2/18 10:08 AM, Kristian Evers wrote:
> > > > Sorry, I don’t know. It’s possible there’s a bug somewhere.
> > > >
> > > > Charles, do you have any insights here?
> > > >
> > > > /Kristian
> > > >
> > > > *Fra:*Markus Metz <markus.metz.giswork at gmail.com
> > <mailto:markus.metz.giswork at gmail.com>>
> > > > *Sendt:* 2. oktober 2018 15:01
> > > > *Til:* Kristian Evers <kreve at sdfe.dk <mailto:kreve at sdfe.dk>>
> > > > *Cc:* GRASS developers list <grass-dev at lists.osgeo.org
> > <mailto:grass-dev at lists.osgeo.org>>; Helmut
> > > > Kudrnovsky <hellik at web.de <mailto:hellik at web.de>>
> > > > *Emne:* Re: [GRASS-dev] area calculations in several GIS
> > > >
> > > >
> > > >
> > > > On Tue, Oct 2, 2018 at 12:08 PM Kristian Evers <kreve at sdfe.dk
> > <mailto:kreve at sdfe.dk>
> > > > <mailto:kreve at sdfe.dk <mailto:kreve at sdfe.dk>>> wrote:
> > > > >
> > > >
> > > > > Use GeographicLib/PROJ as the reference.
> > > >
> > > > Which one? The results of Planimeter from GeographicLib and
> > > > geod_polygonarea() from PROJ are different. If in doubt, use
Planimeter
> > > > as reference?
> > > >
> > > > Markus M
> > > >
> > > > > I would say that this polygon is large, yes. If you want to get
a
> > > > better sense for what constitutes a large polygon in this sense,
try
> > > > creating a series of polygons of varying size (e.g. polygons of
your
> > > > house, your city, your county, your country and your continent) and
> > > > calculate the area with both GRASS and GeographicLib. At some
point the
> > > > calculated areas should start to diverge significantly. I suspect
it
> > > > happens somewhere between county and country.
> > > >
> > > > >
> > > > >
> > > > >
> > > > > /Kristian
> > > > >
> > > > >
> > > > >
> > > > > Fra: Markus Metz <markus.metz.giswork at gmail.com
> > <mailto:markus.metz.giswork at gmail.com>
> > > > <mailto:markus.metz.giswork at gmail.com
> > <mailto:markus.metz.giswork at gmail.com>>>
> > > > > Sendt: 2. oktober 2018 11:51
> > > > > Til: Kristian Evers <kreve at sdfe.dk <mailto:kreve at sdfe.dk>
> > <mailto:kreve at sdfe.dk <mailto:kreve at sdfe.dk>>>
> > > > > Cc: GRASS developers list <grass-dev at lists.osgeo.org
> > <mailto:grass-dev at lists.osgeo.org>
> > > > <mailto:grass-dev at lists.osgeo.org
> > <mailto:grass-dev at lists.osgeo.org>>>; Helmut Kudrnovsky <hellik at web.de
> > <mailto:hellik at web.de>
> > > > <mailto:hellik at web.de <mailto:hellik at web.de>>>
> > > > > Emne: Re: [GRASS-dev] area calculations in several GIS
> > > > >
> > > > >
> > > > >
> > > > >
> > > > >
> > > > > On Tue, Oct 2, 2018 at 8:47 AM Kristian Evers <kreve at sdfe.dk
> > <mailto:kreve at sdfe.dk>
> > > > <mailto:kreve at sdfe.dk <mailto:kreve at sdfe.dk>>> wrote:
> > > > > >
> > > > > > Markus,
> > > > > >
> > > > > >
> > > > > > Thanks for extending the list. I do wonder why Planimeter
gives
> > > > different results than geod_polygonarea(). I’ll run that by Charles
> > > > Karney when I get a chance.
> > > > > >
> > > > >
> > > > > for completeness, I used geod_polygonarea() from proj-5.2.0
> > > > >
> > > > >
> > > > >
> > > > > >> Regarding areas in projected space: UTM is not area true,
laea is
> > > > but is instead scewing angles. You can check that stuff by running
proj
> > > > in very verbose mode.
> > > > > >
> > > > >
> > > > > [...]
> > > > >
> > > > > >
> > > > >
> > > > > > So yes, the LAEA is the better choice of the two but it is
never
> > > > going to represent the true area, especially for large polygons,
and I
> > > > would not advice using it as reference for ellipsoidal area
> > calculations.
> > > > >
> > > > >
> > > > >
> > > > > The question is, what can be used as reference? And is the test
> > > > polygon a "large" polygon, causing "large" deviations from the true
> > area
> > > > when measured in LAEA?
> > > > >
> > > > >
> > > > >
> > > > > Markus M
> > > > >
> > > > > >
> > > > >
> > > > > > /Kristian
> > > > > >
> > > > > >
> > > > > >
> > > > > >
> > > > > >
> > > > > > Fra: Markus Metz <markus.metz.giswork at gmail.com
> > <mailto:markus.metz.giswork at gmail.com>
> > > > <mailto:markus.metz.giswork at gmail.com
> > <mailto:markus.metz.giswork at gmail.com>>>
> > > > > > Sendt: 1. oktober 2018 23:22
> > > > > > Til: GRASS developers list <grass-dev at lists.osgeo.org
> > <mailto:grass-dev at lists.osgeo.org>
> > > > <mailto:grass-dev at lists.osgeo.org
> > <mailto:grass-dev at lists.osgeo.org>>>; Kristian Evers <kreve at sdfe.dk
> > <mailto:kreve at sdfe.dk>
> > > > <mailto:kreve at sdfe.dk <mailto:kreve at sdfe.dk>>>
> > > > > > Cc: Helmut Kudrnovsky <hellik at web.de <mailto:hellik at web.de>
> > <mailto:hellik at web.de <mailto:hellik at web.de>>>
> > > > > > Emne: Re: [GRASS-dev] area calculations in several GIS
> > > > > >
> > > > > >
> > > > > >
> > > > > > Updated list with area calculations for
> > > > > >
> > > > > >
https://gist.github.com/kbevers/207b5bcb9be20e7554abe5f56742ec2c
> > > > > >
> > > > > >
> > > > > >
> > > > > > PROJ [1]: 14,737.935 km^2
> > > > > > Caris LOTS: 14,737 km^2
> > > > > > ArcMap: 14,727.446 km^2
> > > > > > MapInfo: 14,727.352 km^2
> > > > > > GeoMedia: 14,726.443 km^2
> > > > > > Planimeter: 14,722.522 km^2
> > > > > > GRASS GIS: 14.718.098 km^2
> > > > > >
> > > > > > EU LAEA [2]: 14,718.098 km^2
> > > > > > UTM 33 N: 14,707.742 km^2
> > > > > >
> > > > > > QGIS 3.2: 14,652.181 km^2
> > > > > > QGIS 2.8: 14,652.181 km^2
> > > > > >
> > > > > > [1] geodesic.h:geod_polygonarea()
> > > > > > [2] EPSG:3035
> > > > > >
> > > > > >
> > > > > >
> > > > > > In this case, GRASS GIS provides the best match of geodesic
> > area to
> > > > metric area.
> > > > > >
> > > > > >
> > > > > >
> > > > > > @Kristian: are the metric area measurements in "EU LAEA" and
"UTM
> > > > 33 N" suitable as reference?
> > > > > >
> > > > > >
> > > > > >
> > > > > > Discussion started on
> > > > > >
> > > > > >
> > > >
> >
http://osgeo-org.1560.x6.nabble.com/Re-Qgis-user-New-Features-in-Shape-Tools-3-2-0-td5378898.html
> > > > > >
> > > > > >
> > > > > >
> > > > > > Markus M
> > > > > >
> > > > > >
> > > > > >
> > > > > >
> > > > > >
> > > > > > On Mon, Oct 1, 2018 at 2:56 PM Markus Metz
> > > > <markus.metz.giswork at gmail.com
> > <mailto:markus.metz.giswork at gmail.com>
> > <mailto:markus.metz.giswork at gmail.com
> > <mailto:markus.metz.giswork at gmail.com>>>
> > > > wrote:
> > > > > >
> > > > > >
> > > > > >
> > > > > > On Tue, Sep 25, 2018 at 7:38 PM Helmut Kudrnovsky
> > <hellik at web.de <mailto:hellik at web.de>
> > > > <mailto:hellik at web.de <mailto:hellik at web.de>>> wrote:
> > > > > > >
> > > > > > > fyi see
> > > > > > >
> > > >
> >
https://lists.osgeo.org/pipermail/qgis-developer/2018-September/054644.html
> > > > > > >
> > > > > > > with GRASS mentioned
> > > > > > > ------------------
> > > > > > > Kristian Evers:
> > > > > > >
> > > > > > > Right, here are the calculated areas as returned by a
number of
> > > > different
> > > > > > > GIS applications and the planimeter app of GeographicLib for
> > > > reference:
> > > > > > >
> > > > > > > Caris LOTS: 14.737 km^2
> > > > > > > ArcMap: 14.727,446 km^2
> > > > > > > MapInfo: 14.727,352 km^2
> > > > > > > GeoMedia: 14.726,443 km^2
> > > > > > > Planimeter: 14.722,522 km^2
> > > > > > > QGIS 3.2: 14.652,181 km^2
> > > > > > > QGIS 2.8: 14.652,181 km^2
> > > > > >
> > > > > >
> > > > > >
> > > > > > adding to the confusion:
> > > > > >
> > > > > >
> > > > > >
> > > > > > I used the geographiclib API as included in PROJ 5.2.0
following
> > > > the example for geod_polygonarea() in geodesic.h and get
> > > > > >
> > > > > > geographiclib: 14,737.935 km^2
> > > > > >
> > > > > > quite different from
> > > > > >
> > > > > > Planimeter: 14,722.522 km^2
> > > > > >
> > > > > >
> > > > > >
> > > > > > GRASS native gives 14,718.097679
> > > > > >
> > > > > > as also reported by Helmut and Stefan
> > > > > >
> > > > > >
> > > > > >
> > > > > > Markus M
> > > > > >
> > > > > > >
> > > > > >
> > > > > > > The polygon that I have used to get the numbers above can be
> > > > found here:
> > > > > > >
https://gist.github.com/kbevers/207b5bcb9be20e7554abe5f56742ec2c
> > > > > > >
> > > > > > > I am quite confident that GeographicLib delivers the most
> > > > accurate result
> > > > > > > (if you have doubts, this reference [0] should convince
you). As
> > > > can be seen
> > > > > > > from the table above all but QGIS come fairly close. I
expect
> > > > some variation
> > > > > > > in the results as these are numerical approximations,
although I
> > > > think QGIS
> > > > > > > is too far of the mark. My suspicion is that the geodesic
> > > > algorithm used by
> > > > > > > QGIS (and apparently GRASS) is to blame here.
> > > > > > >
> > > > > > > /Kristian
> > > > > > >
> > > > > > > [0] https://arxiv.org/pdf/1102.1215.pdf
> > > > > > > -----------------
> > > > > > >
> > > > > > >
> > > > > > >
> > > > > > >
> > > > > > >
> > > > > > > -----
> > > > > > > best regards
> > > > > > > Helmut
> > > > > > > --
> > > > > > > Sent from:
> > > > http://osgeo-org.1560.x6.nabble.com/Grass-Dev-f3991897.html
> > > > > > > _______________________________________________
> > > > > > > grass-dev mailing list
> > > > > > > grass-dev at lists.osgeo.org <mailto:grass-dev at lists.osgeo.org>
> > <mailto:grass-dev at lists.osgeo.org <mailto:grass-dev at lists.osgeo.org>>
> > > > > > > https://lists.osgeo.org/mailman/listinfo/grass-dev
> > > >
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-dev/attachments/20181004/1ad30c03/attachment-0001.html>
More information about the grass-dev
mailing list