<div dir="ltr"><br><br>On Wed, Oct 3, 2018 at 8:58 AM Kristian Evers <<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a>> wrote:<br>><br>> Thanks for clearing up the confusion, Charles! Looking back in the thread<br>> it seem to have been me who supplied the rhumb line calculation. I guess<br>> that comes down to a cut'n'paste mishap. I'm glad we got to the bottom<br>> of this finally.<br>><br>> I'll also recommend that GRASS and QGIS switch to using either geodesic<br>> or rhumb line area determination. Or even better, offer the ability to<br>> choose various method of calculating areas. Several of the GIS applications<br>> I ran my test polygon through had such features and it did give me a bit<br>> more confidence in those calculations since it was clear what it was doing<br><div>> (albeit not as accurately as I would like in most cases).</div><div><br></div><div>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.</div><div><br></div><div>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?</div><div><br></div><div>Markus M<br></div>><br>>  /Kristian<br>><br>> -----Oprindelig meddelelse-----<br>> Fra: Charles Karney <<a href="mailto:charles.karney@gmail.com">charles.karney@gmail.com</a>> På vegne af Charles Karney<br>> Sendt: 2. oktober 2018 22:46<br>> Til: Markus Metz <<a href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a>><br>> Cc: Kristian Evers <<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a>>; GRASS developers list <<a href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a>>; Helmut Kudrnovsky <<a href="mailto:hellik@web.de">hellik@web.de</a>><br>> Emne: Re: SV: [GRASS-dev] area calculations in several GIS<br>><br>> In comparisons like this, it's probably a good idea to document what<br>> areas the various packages are attempting to compute.<br>><br>> For Planimeter/geod_polygonarea, it is the area of a polygon with edges<br>> given by geodesics.<br>><br>> For GRASS/QGIS, it is the area of a polygon with edges which are<br>> straight lines in the plate-carree projection.<br>><br>> The documentation for GRASS (in lib/gis/area_poly1.c) suggests that this<br>> calculation is a poor man's attempt to compute the area of a polygon<br>> with edges given by rhumb lines.<br>><br>> In fact, the mathematics for computing the area of rhumb polygons is<br>> straightforward, see<br>><br>>    <a href="https://geographiclib.sourceforge.io/html/rhumb.html#rhumbarea">https://geographiclib.sourceforge.io/html/rhumb.html#rhumbarea</a><br>><br>> and GeographicLib's Planimeter utility can perform such calculations<br>> (with the -R flag).  For Kristian's sample polygon we have<br>><br>>    geodesic area = 14737935340.10 m^2<br>>    rhumb area    = 14722522188.60 m^2<br>><br>> This explains the apparent discrepancy between Planimeter and<br>> geod_polygonarea reported earlier.  The figure given for Planimeter was<br>> the rhumb area.<br>><br>> I recommend that GRASS/QGIS be updated either to compute the area of<br>> either geodesic or rhumb polygons.  Surely plate-carree polygons are not<br>> useful?<br>><br>>    --Charles<br>><br>> On 10/2/18 1:55 PM, Markus Metz wrote:<br>> ><br>> ><br>> > On Tue, Oct 2, 2018 at 7:30 PM Charles Karney <<a href="mailto:charles@karney.com">charles@karney.com</a><br>> > <mailto:<a href="mailto:charles@karney.com">charles@karney.com</a>>> wrote:<br>> >  ><br>> >  > Sorry, I'm coming late to this conversation.  The area of the polygon<br>> >  > in the Baltic Sea posted by Kristian is<br>> >  ><br>> >  >    14737935340.098511 m^2<br>> >  ><br>> >  > assuming the WGS84 ellipsoid.  This is accurate to 1 square-mm and was<br>> >  > obtained with the MPFR-enabled version of GeographicLib's Planimeter.<br>> >  ><br>> >  > This is consistent with the result reported by the PROJ.4<br>> >  > geod_polygonarea() and it agrees with the result I get with the regular<br>> >  > Planimeter utility.  I don't why someone else is getting a different<br>> >  > result from Planimeter.<br>> ><br>> > I can confirm this, the results of Planimeter and PROJ 5.2.0 are identical.<br>> ><br>> > Some more confusion:<br>> ><br>> > I created simple boxes for the test polygon, one a bit larger, one a bit<br>> > smaller than the test polygon<br>> ><br>> > GRASS:                13,222.778<br>> > Planimeter:           13,221.965<br>> > geod_polygonarea():   13,221.965<br>> ><br>> > 14.62569 55.36254<br>> > 14.6256944444444 54.0786<br>> > 17.1177 54.0786111111111<br>> > 17.1177777777778 55.3625<br>> ><br>> > GRASS:                22,950.510<br>> > Planimeter:           22,946.901<br>> > geod_polygonarea():   22,946.901<br>> ><br>> > GRASS native geodesic area is not too far off from GeographicLib/PROJ<br>> ><br>> > With the test polygon<br>> > Planimeter/PROJ: 14,737.935 km^2<br>> > GRASS GIS:       14.718.098 km^2<br>> ><br>> > GRASS is much farther off, and here smaller instead of larger than<br>> > Planimeter/PROJ.<br>> ><br>> > I will have a look.<br>> ><br>> > Markus M<br>> ><br>> >  >  It's probably operator error -- but who knows?<br>> >  ><br>> >  >    --Charles<br>> >  ><br>> >  > On 10/2/18 10:08 AM, Kristian Evers wrote:<br>> >  > > Sorry, I don’t know. It’s possible there’s a bug somewhere.<br>> >  > ><br>> >  > > Charles, do you have any insights here?<br>> >  > ><br>> >  > > /Kristian<br>> >  > ><br>> >  > > *Fra:*Markus Metz <<a href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a><br>> > <mailto:<a href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a>>><br>> >  > > *Sendt:* 2. oktober 2018 15:01<br>> >  > > *Til:* Kristian Evers <<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a> <mailto:<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a>>><br>> >  > > *Cc:* GRASS developers list <<a href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a><br>> > <mailto:<a href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a>>>; Helmut<br>> >  > > Kudrnovsky <<a href="mailto:hellik@web.de">hellik@web.de</a> <mailto:<a href="mailto:hellik@web.de">hellik@web.de</a>>><br>> >  > > *Emne:* Re: [GRASS-dev] area calculations in several GIS<br>> >  > ><br>> >  > ><br>> >  > ><br>> >  > > On Tue, Oct 2, 2018 at 12:08 PM Kristian Evers <<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a><br>> > <mailto:<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a>><br>> >  > > <mailto:<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a> <mailto:<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a>>>> wrote:<br>> >  > >  ><br>> >  > ><br>> >  > >  > Use GeographicLib/PROJ as the reference.<br>> >  > ><br>> >  > > Which one? The results of Planimeter from GeographicLib and<br>> >  > > geod_polygonarea() from PROJ are different. If in doubt, use Planimeter<br>> >  > > as reference?<br>> >  > ><br>> >  > > Markus M<br>> >  > ><br>> >  > >  > I would say that this polygon is large, yes. If you want to get a<br>> >  > > better sense for what constitutes a large polygon in this sense, try<br>> >  > > creating a series of polygons of varying size (e.g. polygons of your<br>> >  > > house, your city, your county, your country and your continent) and<br>> >  > > calculate the area with both GRASS and GeographicLib. At some point the<br>> >  > > calculated areas should start to diverge significantly. I suspect it<br>> >  > > happens somewhere between county and country.<br>> >  > ><br>> >  > >  ><br>> >  > >  ><br>> >  > >  ><br>> >  > >  > /Kristian<br>> >  > >  ><br>> >  > >  ><br>> >  > >  ><br>> >  > >  > Fra: Markus Metz <<a href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a><br>> > <mailto:<a href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a>><br>> >  > > <mailto:<a href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a><br>> > <mailto:<a href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a>>>><br>> >  > >  > Sendt: 2. oktober 2018 11:51<br>> >  > >  > Til: Kristian Evers <<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a> <mailto:<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a>><br>> > <mailto:<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a> <mailto:<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a>>>><br>> >  > >  > Cc: GRASS developers list <<a href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a><br>> > <mailto:<a href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a>><br>> >  > > <mailto:<a href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a><br>> > <mailto:<a href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a>>>>; Helmut Kudrnovsky <<a href="mailto:hellik@web.de">hellik@web.de</a><br>> > <mailto:<a href="mailto:hellik@web.de">hellik@web.de</a>><br>> >  > > <mailto:<a href="mailto:hellik@web.de">hellik@web.de</a> <mailto:<a href="mailto:hellik@web.de">hellik@web.de</a>>>><br>> >  > >  > Emne: Re: [GRASS-dev] area calculations in several GIS<br>> >  > >  ><br>> >  > >  ><br>> >  > >  ><br>> >  > >  ><br>> >  > >  ><br>> >  > >  > On Tue, Oct 2, 2018 at 8:47 AM Kristian Evers <<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a><br>> > <mailto:<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a>><br>> >  > > <mailto:<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a> <mailto:<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a>>>> wrote:<br>> >  > >  > ><br>> >  > >  > > Markus,<br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > > Thanks for extending the list. I do wonder why Planimeter gives<br>> >  > > different results than geod_polygonarea(). I’ll run that by Charles<br>> >  > > Karney when I get a chance.<br>> >  > >  > ><br>> >  > >  ><br>> >  > >  > for completeness, I used geod_polygonarea() from proj-5.2.0<br>> >  > >  ><br>> >  > >  ><br>> >  > >  ><br>> >  > >  > >> Regarding areas in projected space: UTM is not area true, laea is<br>> >  > > but is instead scewing angles. You can check that stuff by running proj<br>> >  > > in very verbose mode.<br>> >  > >  > ><br>> >  > >  ><br>> >  > >  > [...]<br>> >  > >  ><br>> >  > >  > ><br>> >  > >  ><br>> >  > >  > > So yes, the LAEA is the better choice of the two but it is never<br>> >  > > going to represent the true area, especially for large polygons, and I<br>> >  > > would not advice using it as reference for ellipsoidal area<br>> > calculations.<br>> >  > >  ><br>> >  > >  ><br>> >  > >  ><br>> >  > >  > The question is, what can be used as reference? And is the test<br>> >  > > polygon a "large" polygon, causing "large" deviations from the true<br>> > area<br>> >  > > when measured in LAEA?<br>> >  > >  ><br>> >  > >  ><br>> >  > >  ><br>> >  > >  > Markus M<br>> >  > >  ><br>> >  > >  > ><br>> >  > >  ><br>> >  > >  > > /Kristian<br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > > Fra: Markus Metz <<a href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a><br>> > <mailto:<a href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a>><br>> >  > > <mailto:<a href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a><br>> > <mailto:<a href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a>>>><br>> >  > >  > > Sendt: 1. oktober 2018 23:22<br>> >  > >  > > Til: GRASS developers list <<a href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a><br>> > <mailto:<a href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a>><br>> >  > > <mailto:<a href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a><br>> > <mailto:<a href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a>>>>; Kristian Evers <<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a><br>> > <mailto:<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a>><br>> >  > > <mailto:<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a> <mailto:<a href="mailto:kreve@sdfe.dk">kreve@sdfe.dk</a>>>><br>> >  > >  > > Cc: Helmut Kudrnovsky <<a href="mailto:hellik@web.de">hellik@web.de</a> <mailto:<a href="mailto:hellik@web.de">hellik@web.de</a>><br>> > <mailto:<a href="mailto:hellik@web.de">hellik@web.de</a> <mailto:<a href="mailto:hellik@web.de">hellik@web.de</a>>>><br>> >  > >  > > Emne: Re: [GRASS-dev] area calculations in several GIS<br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > > Updated list with area calculations for<br>> >  > >  > ><br>> >  > >  > > <a href="https://gist.github.com/kbevers/207b5bcb9be20e7554abe5f56742ec2c">https://gist.github.com/kbevers/207b5bcb9be20e7554abe5f56742ec2c</a><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > > PROJ [1]:    14,737.935 km^2<br>> >  > >  > > Caris LOTS:  14,737 km^2<br>> >  > >  > > ArcMap:      14,727.446 km^2<br>> >  > >  > > MapInfo:     14,727.352 km^2<br>> >  > >  > > GeoMedia:    14,726.443 km^2<br>> >  > >  > > Planimeter:  14,722.522 km^2<br>> >  > >  > > GRASS GIS:   14.718.098 km^2<br>> >  > >  > ><br>> >  > >  > > EU LAEA [2]: 14,718.098 km^2<br>> >  > >  > > UTM 33 N:    14,707.742 km^2<br>> >  > >  > ><br>> >  > >  > > QGIS 3.2:    14,652.181 km^2<br>> >  > >  > > QGIS 2.8:    14,652.181 km^2<br>> >  > >  > ><br>> >  > >  > > [1] geodesic.h:geod_polygonarea()<br>> >  > >  > > [2] EPSG:3035<br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > > In this case, GRASS GIS provides the best match of geodesic<br>> > area to<br>> >  > > metric area.<br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > > @Kristian: are the metric area measurements in "EU LAEA" and "UTM<br>> >  > > 33 N" suitable as reference?<br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > > Discussion started on<br>> >  > >  > ><br>> >  > >  > ><br>> >  > ><br>> > <a href="http://osgeo-org.1560.x6.nabble.com/Re-Qgis-user-New-Features-in-Shape-Tools-3-2-0-td5378898.html">http://osgeo-org.1560.x6.nabble.com/Re-Qgis-user-New-Features-in-Shape-Tools-3-2-0-td5378898.html</a><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > > Markus M<br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > > On Mon, Oct 1, 2018 at 2:56 PM Markus Metz<br>> >  > > <<a href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a><br>> > <mailto:<a href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a>><br>> > <mailto:<a href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a><br>> > <mailto:<a href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a>>>><br>> >  > > wrote:<br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > > On Tue, Sep 25, 2018 at 7:38 PM Helmut Kudrnovsky<br>> > <<a href="mailto:hellik@web.de">hellik@web.de</a> <mailto:<a href="mailto:hellik@web.de">hellik@web.de</a>><br>> >  > > <mailto:<a href="mailto:hellik@web.de">hellik@web.de</a> <mailto:<a href="mailto:hellik@web.de">hellik@web.de</a>>>> wrote:<br>> >  > >  > > ><br>> >  > >  > > > fyi see<br>> >  > >  > > ><br>> >  > ><br>> > <a href="https://lists.osgeo.org/pipermail/qgis-developer/2018-September/054644.html">https://lists.osgeo.org/pipermail/qgis-developer/2018-September/054644.html</a><br>> >  > >  > > ><br>> >  > >  > > > with GRASS mentioned<br>> >  > >  > > > ------------------<br>> >  > >  > > > Kristian Evers:<br>> >  > >  > > ><br>> >  > >  > > > Right, here are the calculated areas as returned by a number of<br>> >  > > different<br>> >  > >  > > > GIS applications and the planimeter app of GeographicLib for<br>> >  > > reference:<br>> >  > >  > > ><br>> >  > >  > > > Caris LOTS: 14.737 km^2<br>> >  > >  > > > ArcMap:     14.727,446 km^2<br>> >  > >  > > > MapInfo:    14.727,352 km^2<br>> >  > >  > > > GeoMedia:   14.726,443 km^2<br>> >  > >  > > > Planimeter: 14.722,522 km^2<br>> >  > >  > > > QGIS 3.2:   14.652,181 km^2<br>> >  > >  > > > QGIS 2.8:   14.652,181 km^2<br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > > adding to the confusion:<br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > > I used the geographiclib API as included in PROJ 5.2.0 following<br>> >  > > the example for geod_polygonarea() in geodesic.h and get<br>> >  > >  > ><br>> >  > >  > > geographiclib: 14,737.935 km^2<br>> >  > >  > ><br>> >  > >  > > quite different from<br>> >  > >  > ><br>> >  > >  > > Planimeter: 14,722.522 km^2<br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > > GRASS native gives 14,718.097679<br>> >  > >  > ><br>> >  > >  > > as also reported by Helmut and Stefan<br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > ><br>> >  > >  > > Markus M<br>> >  > >  > ><br>> >  > >  > > ><br>> >  > >  > ><br>> >  > >  > > > The polygon that I have used to get the numbers above can be<br>> >  > > found here:<br>> >  > >  > > > <a href="https://gist.github.com/kbevers/207b5bcb9be20e7554abe5f56742ec2c">https://gist.github.com/kbevers/207b5bcb9be20e7554abe5f56742ec2c</a><br>> >  > >  > > ><br>> >  > >  > > > I am quite confident that GeographicLib delivers the most<br>> >  > > accurate result<br>> >  > >  > > > (if you have doubts, this reference [0] should convince you). As<br>> >  > > can be seen<br>> >  > >  > > > from the table above all but QGIS come fairly close. I expect<br>> >  > > some variation<br>> >  > >  > > > in the results as these are numerical approximations, although I<br>> >  > > think QGIS<br>> >  > >  > > > is too far of the mark. My suspicion is that the geodesic<br>> >  > > algorithm used by<br>> >  > >  > > > QGIS (and apparently GRASS) is to blame here.<br>> >  > >  > > ><br>> >  > >  > > > /Kristian<br>> >  > >  > > ><br>> >  > >  > > > [0] <a href="https://arxiv.org/pdf/1102.1215.pdf">https://arxiv.org/pdf/1102.1215.pdf</a><br>> >  > >  > > > -----------------<br>> >  > >  > > ><br>> >  > >  > > ><br>> >  > >  > > ><br>> >  > >  > > ><br>> >  > >  > > ><br>> >  > >  > > > -----<br>> >  > >  > > > best regards<br>> >  > >  > > > Helmut<br>> >  > >  > > > --<br>> >  > >  > > > Sent from:<br>> >  > > <a href="http://osgeo-org.1560.x6.nabble.com/Grass-Dev-f3991897.html">http://osgeo-org.1560.x6.nabble.com/Grass-Dev-f3991897.html</a><br>> >  > >  > > > _______________________________________________<br>> >  > >  > > > grass-dev mailing list<br>> >  > >  > > > <a href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a> <mailto:<a href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a>><br>> > <mailto:<a href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a> <mailto:<a href="mailto:grass-dev@lists.osgeo.org">grass-dev@lists.osgeo.org</a>>><br>> >  > >  > > > <a href="https://lists.osgeo.org/mailman/listinfo/grass-dev">https://lists.osgeo.org/mailman/listinfo/grass-dev</a><br>> >  > ></div>