#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>
_______________________________________________ grass-dev mailing list grass-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-dev