[GRASS-dev] area calculations in several GIS

Moritz Lennert mlennert at club.worldonline.be
Tue Oct 2 07:57:08 PDT 2018


On 02/10/18 15:01, Markus Metz wrote:
> 
> 
> On Tue, Oct 2, 2018 at 12:08 PM Kristian Evers <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?

I get the same result as PROJ using geosphere in R:

library(geosphere)
library(rgdal)
poly<-readOGR("ellipsoidal_area.geojson")
areaPolygon(poly)
[1] 14737935340

Exporting the coordinates using:

write.table(poly at polygons[[1]]@Polygons[[1]]@coords[,c(2,1)], 
"polygon.csv", row.names=FALSE)

and feeding them to Planimeter on my Debian testing machine 
(GeographicLib version 1.49), I get:

Planimeter -r --input-file polygon.csv
4017 503378.285237 14737935340.1

So, the result is strictly identical to the PROJ result.

I don't know where the below value of 14,722.522 km^2 comes from.


Moritz


> 
> 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>>
>  > Sendt: 2. oktober 2018 11:51
>  > 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 8:47 AM Kristian Evers <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>>
>  > > Sendt: 1. oktober 2018 23:22
>  > > Til: GRASS developers list <grass-dev at lists.osgeo.org 
> <mailto:grass-dev at lists.osgeo.org>>; Kristian Evers <kreve at sdfe.dk 
> <mailto:kreve at sdfe.dk>>
>  > > Cc: Helmut Kudrnovsky <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>> 
> wrote:
>  > >
>  > >
>  > >
>  > > On Tue, Sep 25, 2018 at 7:38 PM Helmut Kudrnovsky <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>
>  > > > https://lists.osgeo.org/mailman/listinfo/grass-dev
> 
> 
> _______________________________________________
> grass-dev mailing list
> grass-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/grass-dev
> 




More information about the grass-dev mailing list