Make ellipsoidal area ($area) consistent with Proj and PostGIS

[jratike80]

See mailing list discussion https://lists.osgeo.org/pipermail/qgis-user/2021-January/047638.html.

Main points about the test environment:

• QGIS 3.17-dev
• QGIS project and layer are in EPSG:4326
• The ellipsoid for the measurements set into WGS 84 (EPSG:7030)

SQL query that is used as PostGIS test. The same WKT geometry was used in QGIS.

select st_area(st_geomfromtext('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))',4326)::geography , true)

QGIS gives slightly different ellipsoidal areas than programs which are using geographiclib https://geographiclib.sourceforge.io/ for computing geodesic lengths and areas. For example PostGIS and Proj are using geographiclib https://postgis.net/docs/ST_Area.html https://proj.org/geodesic.html?highlight=geodesic%20calcu

It would be nice if QGIS would give same ellipsoidal areas than those two other fine open source programs under the same OSGeo umbrella. It does already give consistent results for geodesic lengths.

This link can be used as a test case https://geographiclib.sourceforge.io/cgi-bin/Planimeter?type=polygon&rhumb=geodesic&input=59.95688345+20.13293641%0D%0A60.47397663+26.94617837%0D%0A62.56499443+29.74782155%0D%0A68.70650340+27.45254202%0D%0A68.24937206+23.75771765%0D%0A65.27444593+25.42698984%0D%0A63.10353609+21.51545237%0D%0A61.12318104+21.40562760%0D%0A60.40477513+19.41123592%0D%0A59.95688345+20.13293641&option=Submit

Results from geographiclib

Perimeter (m)      = 2578086.202363
area (m^2)         = 251199344354.4

Results from QGIS 3.17-Master (with GDAL/OGR 3.3.0dev, PROJ 8.0.0)

Perimeter (m)      = 2578086,202361983
area (m^2)         = 249566957499,7546

Notice that the difference is in area. Perimeters are matching.

Implementation details do not matter but because Proj is always available and geographiclib as well then perhaps it would give some benefit to utilize https://geographiclib.sourceforge.io/1.49/classGeographicLib_1_1PolygonAreaT.html.

[nyalldawson]

As requested on the mailing list thread, can you please check and confirm that both tests are using identical ellipsoid parameters? Can you post the semi major and semi minor axis used by both please?

[jratike80]

Sorry, I thought that the axis were known by the name of the Spheroid.

Semi major and semi minor axis are the same. Verified with command proj -le from the OSGeo4W shell (my QGIS is installed with the same OSGeo4W).

WGS84 a=6378137.0 rf=298.257223563 WGS 84

From PostGIS:
select srtext from spatial_ref_sys where srid=4326;
"GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]"

[nyalldawson]

Can you check the parameters shown in QGIS project properties though?

[jratike80]

Screen capture added to the issue text. I also included the main points from the mail thread.

[agiudiceandrea]

Hi @jratike80, PostGIS uses GeographicLib for geodesic area calculation on ellipsoid/spheroid. SpatiaLite 5.0.0 uses the rttopo library, while SpatiaLite 4.3.0 uses lwgeom lib (that AFAIK is the same of GeographicLib for geodesic area calc). In the mailing list you reported a slight difference between PostGIS and SpatiaLite values: PotGIS: 251199344354.4308 SpatiaLite 5: 251195856999.549927

but using SpatiaLite 4.3.0 you will get the same result of PostGIS SpatiaLite 4.3.0: 251199344354.4308

QGIS instead uses the same algorithm used by GRASS for geodesic area calculation. The relevant part of the code is here and here.

In fact, using GRAS 7.8.4 to perform the area calculation you will get a value almost identical to that obtained using QGIS (>= 3.4) one:

QGIS >= 3.4: 249566957499.7546 GRASS 7.8.4: 249566957499.688.

While a bug in the squared eccentricity calc in the algorithm was fixed since QGIS 3.4 https://github.com/qgis/QGIS/pull/8090 by @nyalldawson, it seems the other parts of the code are in place almost unchanged since very long time.

[jratike80]

Thanks a lot @agiudiceandrea! Do you think that this https://docs.qgis.org/3.16/en/docs/user_manual/working_with_vector/functions_list.html#area is the best place to mention that the GRASS algorithm is used?

I noticed this comment in the sources:

// IMPORTANT // don't change anything here without reporting the changes to upstream (GRASS) // let's all be good opensource citizens and share the improvements!

I will join the grass-devel list and ask what they think about the GeographicLib vs. the GRASS algorithm.

[agiudiceandrea]

@jratike80 Please see also these thread in grass-dev ml http://osgeo-org.1560.x6.nabble.com/area-calculations-in-several-GIS-td5379541.html and in qgis-dev ml http://osgeo-org.1560.x6.nabble.com/area-calculations-in-several-GIS-td5379541.html

[jratike80]

So nothing new on this area. I jumped in just because my now retired colleague from the National Land Survey of Finland asked me why he can't get accurate results from QGIS. He counts on GeographicLib and his own formulas as reference.

There is one thing that I do not understand. The difference is said to be because

For GRASS/QGIS, it is the area of a polygon with edges which are straight lines in the plate-carree projection.

But GeographicLib and QGIS compute equal perimeters. Does it mean that GRASS/QGIS handles the lines differently for length and area, as geodetic lines for lenght but straight in plate-carree for areas?

[roya0045]

5% difference is quite considerable. Good move on raising this issue!

[jratike80]

The difference is not that big, only 0.654% with my test geometry.

[roya0045]

My bad, just had a quick glance, I just checked the first two digit. That's 10 times less.

[nyalldawson]

For reference: while there's clearly justification for the values currently calculated by QGIS, I still think this is still a valid feature request and agree that it would be preferable at some stage to switch from the grass algorithm to instead of Geographic lib for these calculations.

[elpaso]

@manisandro do I remember correctly that one of your derived QGIS applications is using Geographic lib? I kind of remember having seen that library in the dependencies. I'm asking just in case you have already integrated that library into QGIS.

[agiudiceandrea]

Anyway, in the meantime it seems GRASS devs could switch the geodesic area calculation algorithm to use the same C implementation of the GeographicLib algorithm used by PROJ: https://github.com/OSGeo/grass/pull/1283.

[manisandro]

@elpaso Yes we use GeographicLib for KADAS

[elpaso]

@elpaso Yes we use GeographicLib for KADAS

Can you give us some more background about that choice and the advantages/disadvantages of using that library instead of the current approach? (If it is related to the issue in this ticket of course)

[nyalldawson]

For reference: https://github.com/qgis/QGIS/pull/41726

[agiudiceandrea]

Is this feature request resolved now that #41726 was merged in master?

[nyalldawson]

Yes, it's resolved!