[GRASS-dev] [GRASS GIS] #3356: v.to.db: incorrect area calculations in lat-long location

[GRASS GIS]

#3356: v.to.db: incorrect area calculations in lat-long location
--------------------------+-----------------------------------
  Reporter:  mlennert     |      Owner:  grass-dev@…
      Type:  defect       |     Status:  new
  Priority:  normal       |  Milestone:  7.4.0
 Component:  Vector       |    Version:  svn-trunk
Resolution:               |   Keywords:  v.to.db area lat-long
       CPU:  Unspecified  |   Platform:  Unspecified
--------------------------+-----------------------------------

Comment (by mlennert):

 Replying to [comment:8 mmetz]:
 > Replying to [comment:7 mmetz]:
 > > Replying to [comment:6 mlennert]:
 > > > Replying to [comment:5 annakrat]:
 > > >
 > > > > But the area computation problem must be somewhere in
 [https://grass.osgeo.org/programming7
 > > > >/area__poly1_8c.html#af6f1f53bacc34249be98006c95369695
 G_ellipsoid_polygon_area], probably related
 > > > > to very small numbers, but I couldn't pinpoint the problem. There
 is something specific about this
 > > > > polygon, when I draw a similar one, it gives more reasonable
 results.
 > > >
 > > > So, yes, this definitely seems to be a precision issue, but at this
 stage I can't really see where in the source code, and don't have much
 time to delve into it in greater detail. Maybe MarkusM has an idea ?
 > >
 > > The problem seems to be in G_ellipsoid_polygon_area() at
 [https://trac.osgeo.org/grass/browser/grass/trunk/lib/gis/area_poly1.c#L161
 L161]. If dy is much smaller than dx, dx / dy becomes very large and
 erroneus values might sneak in. A simple solution could be to set dy to
 zero if dy is very small, but then how to define "very small"?
 Interestingly, dx must not be snapped to zero, otherwise I get complete
 nonsense results for the test shapefile.
 >
 > Apparently small differences in latitudes of adjacent vertices can cause
 a wrong latitude correction in G_ellipsoid_polygon_area(). Please try
 trunk r71167.
 >
 > Note that this does not affect areas created with v.in.region, but this
 fix also affects larger areas such as country boundaries.

 Thanks for the quick find !

 I tested and compared with the results of the R
 [https://cran.r-project.org/package=geosphere geosphere] package (which
 apparently uses [https://geographiclib.sourceforge.io/ GeographicLib]:

 {{{
 v.to.db poly -p op=area --q
 1|0.126662200952979
 }}}

 In R:

 {{{
 crswgs84=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
 poly=readShapePoly("GRASS_area_problem/training_single.shp",proj4string=crswgs84,verbose=TRUE)
 areaPolygon(poly)
 [1] 0.1262853
 }}}

 Don't know if there are other "authoritative" programs to measure the area
 of this polygon ? At one point, I guess it comes down to floating point
 precision used in the different stages of the algorithm...

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/3356#comment:9>
GRASS GIS <https://grass.osgeo.org>