libgis: use GeographicLib from PROJ

[metzm]

• use GeographicLib for geodesic distance and area calculations (see [Feat] Use GeographicLib for computing ellipsoidal areas #1235)

[metzm]

CentOS 7 is out with this PR because it gets proj-4.8.0 from EPEL which is too old...

The big question is, if we accept the fact that this PR does not work with stock CentOS 7, which the GRASS community needs to answer.

[jratike80]

@metzm Markus, you asked me to review but even I am a somewhat experienced GIS user I unfortunately understand very little about coding.

[neteler]

Maybe @kbevers is willing to take a look at this PR?

[kbevers]

Looks about right to me but I am in no way an expert in the geodesic code so take my blessings with a grain of salt

[nilason]

I have tested with testfile linked in QGIS ml:
ellipsoidal_area.geojson.txt
(added .txt file extension as Github doesn't support geojson).

With this PR (proj=7.2.1)
14737935340.099 sq.m

With G 7.8.5 (proj=7.1.1)
14718097678.906 sq.m

QGIS 3.10.11 (proj=5.2.0)
14718097678.949 sq.m

For nc_spm_08 / boundary_county@PERMANENT
feature: NAME:WAKE
PR: 2219365235.93 sq.m
7.8.5: 2219365235.93 sq.m

The results seem to me to be correct, with no regression.
I didn't encounter any other issue either.

[mlennert]

I've tested and can confirm that in current trunk, results between GRASS GIS and the GeographicLib implementation in R (geosphere) are identical (to the sq m):

v.in.ogr ellipsoidal_area.geojson out=ell_test
v.to.db -p ell_test op=area --q
1|14737935340.0989

library(geosphere)
library(rgdal)
poly<-readOGR("ellipsoidal_area.geojson")
areaPolygon(poly)
14737935340

Results in GRASS GIS also seem stable when generalizing and more generalized versions do concur with the geosphere output:

v.generalize ell_test method=douglas out=ell_test_00000001 thresh=0.00000001 --o --q
v.to.db -p ell_test_00000001 op=area --q
1|14737935340.0999

v.generalize ell_test method=douglas out=ell_test_0000001 thresh=0.0000001 --o --q
v.to.db -p ell_test_0000001 op=area --q
1|14737935354.4453

library(geosphere)
library(rgdal)
poly<-readOGR("ell_test_0000001.json")
areaPolygon(poly)
[1] 14737935354

And testing with the data described in this old trac ticket the results are fairly stable:

v.in.ogr GRASS_area_problem/training_single.shp out=poly --o
v.generalize poly method=douglas out=poly_gen_000000001 thresh=0.000000001 --o --q
v.generalize poly method=douglas out=poly_gen_000000003 thresh=0.000000003 --o --q
v.generalize poly method=douglas out=poly_gen_00000001 thresh=0.00000001 --o --q
v.to.db -p poly op=area --q
1|0.126277758915019
v.to.db -p poly_gen_000000001 op=area --q
1|0.126328429
v.to.db -p poly_gen_000000003 op=area --q
1|0.12635081776554964583
v.to.db -p poly_gen_00000001 op=area --q
1|0.126351168928104

But results are not identical to geosphere results:

v.out.ogr poly format=GeoJSON out=poly.json
v.out.ogr poly_gen_000000001 format=GeoJSON out=poly_gen_000000001.json
v.out.ogr poly_gen_000000003 format=GeoJSON out=poly_gen_000000003.json
v.out.ogr poly format=GeoJSON out=poly.json

library(geosphere)
library(rgdal)
poly<-readOGR("poly.json")
areaPolygon(poly)
[1] 0.126354
poly<-readOGR("poly_gen_000000001.json")
areaPolygon(poly)
[1] 0.1263502
poly<-readOGR("poly_gen_000000003.json")
areaPolygon(poly)
[1] 0.1263632
poly<-readOGR("poly_gen_00000001.json")
areaPolygon(poly)
[1] 0.1263439

Any idea why ?

[metzm]

Testing the example of mlennert, one reason for different results with R geosphere is that a and f where not given, however, using the same a and f as in GRASS I still don't get the same results, and my results with GRASS differ from the results of mlennert with GRASS:

mlennert GRASS master + PROJ ?
metzm GRASS master + PROJ 7.2.1
R geosphere without a, f
R geosphere with a = 6.37814e+06, f = 0.00335281

poly

0.126277758915019           mlennert
0.126272547475285           metzm
0.126354                    R geosphere without a, f
0.126373                    R geosphere with a, f

poly_gen_000000001

0.126328429                 mlennert
0.126329002640606           metzm
0.1263502                   R geosphere without a, f
0.1263837                   R geosphere with a, f

poly_gen_000000003

0.12635081776554964583      mlennert
0.126352616683089           metzm
0.1263632                   R geosphere without a, f
0.1264062                   R geosphere with a, f

poly_gen_00000001

0.126351168928104           mlennert
0.126362257568278           metzm
0.1263439                   R geosphere without a, f
0.1263677                   R geosphere with a, f

If mlennert also used PROJ 7.2.1, there must be some numerical instability in the GeographicLib as included in PROJ, otherwise differences between mlennert and metzm might be due to different PROJ versions.

Altogether, the observed differences can not even be explained with single precision floating point limitations. GeographicLib can no longer be used as reference because GeographicLib itself produces different results on different systems.

My results with the R geosphere package are identical to those of mlennert, but according to the geosphere documentation, a and f need to be specified.

[metzm]

For my tests with R geosphere, the precision of a and f was too low. Using more precise values a = 6378137, f = 0.00335281066474746, the results of R geosphere with and without a and f are identical.

[metzm]

There is also the Planimeter tool of GeographicLib. Thus there are three different implementations of the algorithm of Charles Karney: Planimeter, geodesic.[c|h] in PROJ and R geosphere. Unfortunately with four different results:

poly
0.126272547475285           GRASS + PROJ metzm
0.126277758915019           GRASS + PROJ mlennert
0.126354                    R geosphere
0.12638			    Planimeter

Differences begin to occur already at the fourth significant decimal digit, latest at the sixth significant decimal digit. This is IMHO too high.

[agiudiceandrea]

There is also the Planimeter tool of GeographicLib.

Have you also used the "Planimeter" tool with the "-E" option?

Anyway, I think geodesic/ellipsoidal area calculation are not actual useful for polygons of very small extension (like the "training_single" layer which is about 0.1 square meters).

[metzm]

@agiudiceandrea @kbevers

There is also the Planimeter tool of GeographicLib.

Have you also used the "Planimeter" tool with the "-E" option?

Same result. Does not answer the question why Planimeter and PROJ geodesic provide different results.

Anyway, I think geodesic/ellipsoidal area calculation are not actual useful for polygons of very small extension (like the "training_single" layer which is about 0.1 square meters).

But PROJ geodesic.h states that "for the WGS84 ellipsoid, the errors are less than 15 nanometers". Thus I expect accurate and consistent (!) results for our example which is with > 0.1 spm huge relative to the given errors.

I could create a new issue in PROJ for the inconsistency between GeographicLib and geodesic.[c|h] in PROJ.

[kbevers]

I think it's time to call in the cavalry. @cffk, above are showcased inconsistencies between various setups using geodesic calculations from GeographicLib - any insight as to why?

[cffk]

OK, I'll take a look. The intention is that PROJ's geodesic calculations (C code) should mirror very accurately what you get from GeographicLib (C++ code).

[agiudiceandrea]

But PROJ geodesic.h states that "for the WGS84 ellipsoid, the errors are less than 15 nanometers". Thus I expect accurate and consistent (!) results for our example which is with > 0.1 spm huge relative to the given errors.

I would also expect consistent results.

Anyway, the Planimeter online tool description at https://geographiclib.sourceforge.io/cgi-bin/Planimeter states that: "The result for the perimeter is accurate to about 15 nm per vertex. The result for the area is accurate to about 0.1 m² per vertex."

The same is stated in C. F. F. Karney, Algorithms for geodesics, J. Geodesy 87, 43-55 (2013): "The round-off errors in the direct and inverse methods are less than 15 nanometers and the error in the computation of the area S₁₂ is about 0.1 m²."

[mlennert]

If mlennert also used PROJ 7.2.1, there must be some numerical instability in the GeographicLib as included in PROJ, otherwise differences between mlennert and metzm might be due to different PROJ versions.

GRASS ll_wgs85/test_proj:~ > g.version -rbe
GRASS 7.9.dev (2021)

 ./configure  --prefix=/usr/lib --sysconfdir=/etc --sharedstatedir=/var --enable-socket --enable-shared --enable-largefile --with-openmp --with-postgres --with-pthread --with-cxx --with-x --with-gdal --with-freetype --with-motif --with-readline --with-nls --with-odbc --with-sqlite --with-freetype-includes=/usr/include/freetype2 --with-tcltk-includes=/usr/include/tcl --with-postgres-includes=/usr/include/postgresql --with-proj-share=/usr/share/proj --with-python=/usr/bin/python-config --with-cairo --with-geos --with-blas --with-lapack --with-zstd --with-pdal=/usr/bin/pdal-config --with-netcdf=/usr/bin/nc-config
libgis revision: 8b4c4986e
libgis date: 2021-02-08T07:13:27+00:00
PROJ: 7.2.1
GDAL/OGR: 3.2.1
GEOS: 3.9.0
SQLite: 3.34.1

GRASS ll_wgs85/test_proj:~ > v.to.db -p poly op=area --q
1|0.126277758915019

I'm on Debian testing. Is there any env variable I need to set or something ? Is there a way to test that I'm really using the GeographicLib ?

[mlennert]

But PROJ geodesic.h states that "for the WGS84 ellipsoid, the errors are less than 15 nanometers". Thus I expect accurate and consistent (!) results for our example which is with > 0.1 spm huge relative to the given errors.

I would also expect consistent results.

Anyway, the Planimeter online tool description at https://geographiclib.sourceforge.io/cgi-bin/Planimeter states that: "The result for the perimeter is accurate to about 15 nm per vertex. The result for the area is accurate to about 0.1 m2 per vertex."

0.1 m2 per vertex seems quite a lot actually, especially for small polygons. Is this completely independent of the size ?

IIUC, the small test polygon used above has 45 vertices, so that would mean an expected error of 4.5m2 ?

[metzm]

I'm on Debian testing. Is there any env variable I need to set or something ? Is there a way to test that I'm really using the GeographicLib ?

If you use this PR, geodesic.[c|h] (the C version of GeographicLib) from PROJ is automatically used.

[nilason]

I'm on Debian testing. Is there any env variable I need to set or something ? Is there a way to test that I'm really using the GeographicLib ?

If you use this PR, geodesic.[c|h] (the C version of GeographicLib) from PROJ is automatically used.

Judging from libgis revision: 8b4c498 it is not this PR

[cffk]

The "exact" area of the 44-sided training_single polygon is

-0.126347391348170 m^2

if the coordinates in the KML file are treated as exact decimal numbers.
This calculation is performed by the MPFR-enabled version of
GeographicLib for the WGS84 ellipsoid. The plain double version of
Planimeter and a simple planimeter implementation using PROJ both give

-0.1263768025 m^2

so the error is 30 mm^2 or about 1 mm^2 per vertex. The error in the
"GRASS + PROJ metzm" calculation is 75 mm^2 in the other direction. I
don't think you should read too much into the discrepancies between the
various PROJ based calculations. These are probably attributable to
differences in the compiler

My estimate of the error of 0.1 m^2 per vertex tries to capture the
worst case for arbitrary polygons. This reflects the underlying method,
namely applying the trapezium rule to the polygon with the equator as
the baseline. The width of the trapezium is computed by a geodesic
calculation with error 10e-9 m; the height of the trapezium is at most
10e7 m. So the error per side is at most

10e-9 m x 10e7 m = 0.1 m^2

Almost certainly, the calculation of the area benefits from a
cancellation of errors. It might be worth trying to tighten up the
bounds on the errors; e.g., using approximate squares of edges between
10 cm and 100 km. I can generate a set of test data along these lines
together with accurate areas.

As to the question of whether 100 mm^2 = (1 cm)^2 (let's say) is too big
an error for this polygon...

It would be difficult to wring much more accuracy out of a method using
global coordinates (latitude + longitude) and double precision
arithmetic. It would require, for example, reformulating the geodesic
calculations to give high relative accuracy for short lines.
(Currently, they just give high absolute accuracy, which, at 10e-9 m,
is good enough for most applications.) For what it's worth

A / 2^52 = 0.11 m^2

where A is the area of the ellipsoid.

[agiudiceandrea]

The "exact" area of the 44-sided training_single polygon is

-0.126347391348170 m^2

For reference:

• using training_single.shp in GRASS_area_problem.zip from https://trac.osgeo.org/grass/ticket/3356

  • GRASS 7.8.5 [1]
    v.to.db -p training_single_shp op=area --q
    1|0.126347391488286

  • QGIS 3.16.3 [2]
    $area = 0.12634739148828641

• using training_single.kml from https://trac.osgeo.org/grass/ticket/3356

  • GRASS 7.8.5 [1]
    v.to.db -p training_single_kml op=area --q
    1|0.126347391454861

  • QGIS 3.16.3 [2]
    $area = 0.12634739145486146

So, the value of the area of this particular polygon, calculated with the current GRASS/QGIS implementation of the geodesic area calculation algorithm, differs by about 1E-10 square meters from the "exact" value.

[1] GRASS version: 7.8.5 - Code revision: 2b6ab28 - Windows 7 64 bit (OSGeo4W 64 bit)
[2] QGIS 3.16.3 - Code revision: 94ac9f21b8 - Windows 7 64 bit (OSGeo4W 64 bit)

[jratike80]

How about the polygon from the QGIS issue qgis/QGIS#40888? For that the difference in area between current QGIS and GeographicLib was 0.654%.

POLYGON ((
 20.13293641   59.95688345,
 26.94617837   60.47397663,
 29.74782155   62.56499443,
 27.45254202   68.70650340,
 23.75771765   68.24937206,
 25.42698984   65.27444593,
 21.51545237   63.10353609,
 21.40562760   61.12318104,
 19.41123592   60.40477513,
 20.13293641   59.95688345))

[metzm]

This tiny test polygon with a size of about 0.1263 sqm seems to hit the limits of floating point calculations and trigonometric functions. For larger areas, the Planimeter tool of GeographicLib and the corresponding C version in PROJ provide identical results.

Let's assume that GeographicLib provides the most accurate results for geodesic distance and area calculations available in the osgeo world. The different results of the Planimeter tool of GeographicLib and the corresponding C version in PROJ might be explained by different compilers and/or compiler flags. In my case, using -fno-fast-math with PROJ changes the size of the tiny test polygon from 0.126272547475285 to 0.12634739123375 which is very close to 0.126347391348170, the "exact" value calculated by @cffk. Planimeter's 0.12638 from the package manager of my Linux distro is thus a bit off.

Compilers and compiler flags seem to make a difference, explaining any inconsistencies for extreme cases between he Planimeter tool of GeographicLib and the corresponding C version in PROJ.

AFAIK, QGIS uses the same (or very similar) algorithm to calculate area sizes as the GRASS native algorithm and is thus not suitable for comparison of accurate results, only for comparison of consistency.

[cffk]

How about the polygon from the QGIS issue qgis/QGIS#40888? For that the difference in area between current QGIS and GeographicLib was 0.654%.

POLYGON ((
 20.13293641   59.95688345,
 26.94617837   60.47397663,
 29.74782155   62.56499443,
 27.45254202   68.70650340,
 23.75771765   68.24937206,
 25.42698984   65.27444593,
 21.51545237   63.10353609,
 21.40562760   61.12318104,
 19.41123592   60.40477513,
 20.13293641   59.95688345))

The "exact" area here is (assuming coordinates are lon,lat and WGS84):

251199344354.42985178 m^2

[cffk]

One possible source of confusion is that the online planimeter tool is compiled to use long double so you get 11 bits more accuracy compared to the Planimeter tool you get with most distributions of GeographicLib.

[jratike80]

So, the value of the area of this particular polygon, calculated with the current GRASS/QGIS implementation of the geodesic area calculation algorithm, differs by about 1E-10 square meters from the "exact" value.

If you mean that the difference is perhaps too small for having any practical meaning, it may be true with this tiny polygon but have a look at the half-a-Finland polygon that I pasted as a comment. That shows 0.654% difference in areas.

"exact" area (m^2)          = 251199344354.42985178 
QGIS 3.17 with GRASS method = 249566957499.7546

[nilason]

So, the value of the area of this particular polygon, calculated with the current GRASS/QGIS implementation of the geodesic area calculation algorithm, differs by about 1E-10 square meters from the "exact" value.

If you mean that the difference is perhaps too small for having any practical meaning, it may be true with this tiny polygon but have a look at the half-a-Finland polygon that I pasted as a comment. That shows 0.654% difference in areas.

"exact" area (m^2)          = 251199344354.42985178 
QGIS 3.17 with GRASS method = 249566957499.7546

Using this PR (with PROJ 7.2.1) the same polygon gives the result 251199344354.431 which is pretty close to the "exact" area.

[metzm]

"exact" area (m^2)          = 251199344354.42985178 
QGIS 3.17 with GRASS method = 249566957499.7546

Using this PR (with PROJ 7.2.1) the same polygon gives the result 251199344354.431 which is pretty close to the "exact" area.

Close enough for me.

Now that I got my PROJ compiler flags right, I get in all reported cases results very close to these "exact" results of GeographicLib. For me, the inconsistencies are resolved.

[nilason]

The relatively big difference between the "old" grass/qgis and the new calculation (this PR) for the Finland polygon, may perhaps display the limitations of the old method.

Testing a polygon (~ Uganda) along the equator:

equat.geojson.txt

gives ~0.012% difference:

qgis (=old grass)   =  192668499676.993
this PR             =  192692521685.248

[ninsbl]

Bumping up milestone as 8.0.1 is due in two days, while this has not been part of RC1 and there has not been activity for some time.

[cffk]

FYI, I comment on GRASS's current method of calculating areas (which assumes straight edges in a plate carree projection) in

https://geographiclib.sourceforge.io/C++/2.4/rhumb.html#platecarreearea

Bottom line:

• the method is correct;
• it could be improved;
• nevertheless I recommend switching to geodesic edges.

Is this PR likely to make it into GRASS 8.5.0? It seems that it's just a matter of "doing it".

[echoix]

Solved conflicts of the existing code of the PR.

[petrasovaa]

This is causing test failure in v.import (across platforms):

FAIL: test_import_gpkg_wgs84 (__main__.TestVImport)
Import GPKG in same proj, default params
----------------------------------------------------------------------
Traceback (most recent call last):
  File "scripts/v.import/testsuite/test_v_import.py", line 81, in test_import_gpkg_wgs84
    self.assertVectorFitsRegionInfo(
  File "etc/python/grass/gunittest/case.py", line 420, in assertVectorFitsRegionInfo
    self.assertModuleKeyValue(
  File "etc/python/grass/gunittest/case.py", line 283, in assertModuleKeyValue
    self.fail(self._formatMessage(msg, stdMsg))
AssertionError: v.info map=test_v_import_imported layer=1 format=plain -g difference:
mismatch values (key, reference, actual): [('north', 227744.8, 227201.923076887), ('south', 215259.6, 219012.237762206)]
command: v.info map=test_v_import_imported layer=1 format=plain -g {'map': 'test_v_import_imported', 'flags': 'g'}

@metzm any idea?

v.proj uses G_distance to compute segment lengths for densification, so different results probably make sense, but I was surprised that they differ so much.