v.to.db: incorrect area calculations in lat-long location

[mlennert]

​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.

[mlennert]

Shapefile with polygon causing the problem

[hellik]

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²

[hellik]

v.in.region of training_single

[hellik]

Replying to hellik:

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²

added 2 kml files (training_single.kml and bbox of training_single by v.in.region) to see the polygon in real life conditions

[hellik]

Replying to hellik:

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²

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, reprojet this shapefile to EPSG:32743, creating a EPSG:32743-location by ogr2ogr, 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?

[hellik]

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.

[annakrat]

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.

[mlennert]

Replying to annakrat:

I fixed the rendering issue in r71163 and r71164.

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 ?

[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.

[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.

[mlennert]

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...

[mlennert]

Can we close this ticket ?

[hellik]

Replying to mlennert:

Can we close this ticket ?

Never tested it with other examples, but it seems to work. Is it already backported?

[mlennert]

Replying to hellik:

Replying to mlennert:

Can we close this ticket ?

Never tested it with other examples, but it seems to work. Is it already backported?

Yes, in r71260.

[hellik]

Replying to mlennert:

Replying to hellik:

Replying to mlennert:

Can we close this ticket ?

Never tested it with other examples, but it seems to work. Is it already backported?

Yes, in r71260.

Then lets close it