Opened 3 years ago
Closed 3 years ago
#3356 closed defect (fixed)
v.to.db: incorrect area calculations in lat-long location
Reported by: | mlennert | Owned by: | |
---|---|---|---|
Priority: | normal | Milestone: | 7.4.0 |
Component: | Vector | Version: | svn-trunk |
Keywords: | v.to.db area lat-long | Cc: | |
CPU: | Unspecified | Platform: | Unspecified |
Description
As reported by James Duffy on the ML, using the attached shapefile, raises the following issues:
When I import the polygon in an EPSG 4326 location, I get the same result as reported by James:
v.report training op=area u=me cat|id|area 1|1|38.2256243922775
But when I reproject it to a UTM 43N (EPSG 32743) location, I get a different area, which is close to what QGIS gives as area:
v.proj location=LL_WGS84 mapset=mlennert input=training output=training_reproj_grass v.report training_reproj_grass op=area u=me cat|id|area 1|1|0.126369981910102
Interestingly, when I create a buffer around the area in the lat-long location:
v.buffer training dist=0.0001 out=training_buff_0_0001
The area measurement in GRASS is again very close to the one in QGIS:
v.report training_buff_0_0001 op=area u=me cat|area 1|419.188654810375
$area in field calculator in QGIS: 425.1112801
In addition, in the ESPG 4326 location, it is impossible to zoom to the original polygon as it seems to go below the zoom capacity of the GUI:
"Failed to run command 'd.vect map=training@mlennert type=point,line,boundary,area,face width=1'. Details: GRASS_INFO_ERROR(22724,1): Invalid n-s resolution field: 0.000000 "
Although these two might be unrelated, I do feel that there might be an issue with floating point precision somewhere.
Attachments (3)
Change History (17)
Changed 3 years ago by
Attachment: | GRASS_area_problem.zip added |
---|
comment:1 follow-ups: 2 3 Changed 3 years ago by
Replying to mlennert:
$area in field calculator in QGIS: 425.1112801
just tested it with
QGIS version 2.18.9 (OTF enabled)
1 34.815138348493
and when I do manual measure in QGIS (OTF on)
32,4 m²
Changed 3 years ago by
Attachment: | training_single.kml added |
---|
comment:2 Changed 3 years ago by
Replying to hellik:
Replying to mlennert:
$area in field calculator in QGIS: 425.1112801just tested it with
QGIS version 2.18.9 (OTF enabled)1 34.815138348493and when I do manual measure in QGIS (OTF on)
32,4 m²
added 2 kml files (training_single.kml and bbox of training_single by v.in.region) to see the polygon in real life conditions
comment:3 Changed 3 years ago by
Replying to hellik:
Replying to mlennert:
$area in field calculator in QGIS: 425.1112801just tested it with
QGIS version 2.18.9 (OTF enabled)1 34.815138348493and when I do manual measure in QGIS (OTF on)
32,4 m²
now tested
g.region -p vector=training_single@test projection: 3 (Latitude-Longitude) zone: 0 datum: wgs84 ellipsoid: wgs84 north: 0:13:32.187577N south: 0:13:32.175318N west: 73:12:48.504873E east: 73:12:48.518218E nsres: 0:00:00.012259 ewres: 0:00:00.013345 rows: 1 cols: 1 cells: 1
v.in.region output=reg_training_single
v.report map=reg_training_single@test option=area units=meters cat|area 1|0.155378405506674
then export reg_training_single to a shapefile, reproject this shapefile to EPSG:32743 by ogr2ogr, creating a EPSG:32743-location, import the reprojected shapefile into the new location
v.report map=myreg2_32743@data option=area units=meters cat||area 1|0.155406186822802
now area of the bbox of training_single in wgs84 and EPSG:32743 seems to be similar.
check this bbox vector area also in qgis (OTF enabled).
cat areawgs84 areareproj 1 0.1574797419 0.157479741
QGIS and GRASS area seems similar and reasonable, when you're looking at the kml-file in real life conditions.
these numbers are similar to the number in the original report
v.report training_reproj_grass op=area u=me cat|id|area 1|1|0.126369981910102
but not similar to
v.report training op=area u=me cat|id|area 1|1|38.2256243922775
or
$area in field calculator in QGIS: 425.1112801
any idea?
comment:4 Changed 3 years ago by
g.region -p vector=training_single at test projection: 3 (Latitude-Longitude) zone: 0 datum: wgs84 ellipsoid: wgs84 north: 0:13:32.187577N south: 0:13:32.175318N west: 73:12:48.504873E east: 73:12:48.518218E nsres: 0:00:00.012259 ewres: 0:00:00.013345 rows: 1 cols: 1 cells: 1
looking at resolution and the w-e or n-s borders:
the difference seems to be at the 3rd decimal number of the seconds.
That's maybe the reason the wxgui display can't map the polygon.
And also area calculations may have some issues.
comment:5 follow-up: 6 Changed 3 years ago by
I fixed the rendering issue in r71163 and r71164. But the area computation problem must be somewhere in 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.
comment:6 follow-up: 7 Changed 3 years ago by
Replying to annakrat:
Yes, works great now, thanks !
But the area computation problem must be somewhere in [https://grass.osgeo.org/programming7
/areapoly1_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.
Yes:
v.in.ogr training_single.shp out=poly v.to.db -p poly op=area Reading areas... 100% cat|area 1|38.2256243922775 g.region vect=poly v.in.region out=test v.to.db -p test op=area Reading areas... 100% cat|area 1|0.155378405506674
Another interesting test:
v.generalize poly method=douglas out=poly_gen thresh=0.000000001 --o v.to.db -p poly_gen op=area --q 1|38.2255972503724 v.generalize poly method=douglas out=poly_gen thresh=0.00000001 --o v.to.db -p poly_gen op=area --q 1|2.27719829466825 v.generalize poly method=douglas out=poly_gen thresh=0.00000002 v.to.db -p poly_gen op=area --q 1|0.0634717598906464 v.generalize poly method=douglas out=poly_gen thresh=0.00000003 --o v.to.db -p poly_gen op=area --q 1|0.134256325986436 v.generalize poly method=douglas out=poly_gen thresh=0.0000001 --o v.to.db -p poly_gen op=area --q 1|0.000968180796418213
IOW: just very slight generalization gives much better area calculations, but very small differences in generalization can lead to quite erratic area calculation 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 ?
comment:7 follow-up: 8 Changed 3 years ago by
Replying to mlennert:
Replying to annakrat:
But the area computation problem must be somewhere in [https://grass.osgeo.org/programming7
/areapoly1_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 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.
comment:8 follow-up: 9 Changed 3 years ago by
Replying to mmetz:
Replying to mlennert:
Replying to annakrat:
But the area computation problem must be somewhere in [https://grass.osgeo.org/programming7
/areapoly1_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 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.
comment:9 Changed 3 years ago by
Replying to mmetz:
Replying to mmetz:
Replying to mlennert:
Replying to annakrat:
But the area computation problem must be somewhere in [https://grass.osgeo.org/programming7
/areapoly1_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 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 geosphere package (which apparently uses 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...
comment:11 follow-up: 12 Changed 3 years ago by
Replying to mlennert:
Can we close this ticket ?
Never tested it with other examples, but it seems to work. Is it already backported?
comment:12 follow-up: 13 Changed 3 years ago by
comment:13 Changed 3 years ago by
comment:14 Changed 3 years ago by
Resolution: | → fixed |
---|---|
Status: | new → closed |
Shapefile with polygon causing the problem