Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
jratike80 opened this issue Jan 7, 2021 · 19 comments
Closed

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

jratike80 opened this issue Jan 7, 2021 · 19 comments
Labels
Feature Request Feedback

Comments

@jratike80
Copy link

@jratike80 jratike80 commented Jan 7, 2021

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)

image

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
Copy link
Collaborator

@nyalldawson nyalldawson commented Jan 7, 2021

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?

@nyalldawson nyalldawson added the Feedback label Jan 7, 2021
@jratike80
Copy link
Author

@jratike80 jratike80 commented Jan 7, 2021

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
Copy link
Collaborator

@nyalldawson nyalldawson commented Jan 7, 2021

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

@jratike80
Copy link
Author

@jratike80 jratike80 commented Jan 7, 2021

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

@gioman gioman removed the Feedback label Jan 7, 2021
@agiudiceandrea
Copy link
Contributor

@agiudiceandrea agiudiceandrea commented Jan 7, 2021

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 #8090 by @nyalldawson, it seems the other parts of the code are in place almost unchanged since very long time.

@jratike80
Copy link
Author

@jratike80 jratike80 commented Jan 7, 2021

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.

@jratike80
Copy link
Author

@jratike80 jratike80 commented Jan 7, 2021

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
Copy link
Contributor

@roya0045 roya0045 commented Jan 7, 2021

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

@jratike80
Copy link
Author

@jratike80 jratike80 commented Jan 7, 2021

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

@roya0045
Copy link
Contributor

@roya0045 roya0045 commented Jan 7, 2021

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

@nyalldawson
Copy link
Collaborator

@nyalldawson nyalldawson commented Jan 7, 2021

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
Copy link
Contributor

@elpaso elpaso commented Feb 19, 2021

@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
Copy link
Contributor

@agiudiceandrea agiudiceandrea commented Feb 19, 2021

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: OSGeo/grass#1283.

@manisandro
Copy link
Member

@manisandro manisandro commented Feb 19, 2021

@elpaso Yes we use GeographicLib for KADAS

@elpaso
Copy link
Contributor

@elpaso elpaso commented Feb 20, 2021

@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
Copy link
Collaborator

@nyalldawson nyalldawson commented Feb 22, 2021

For reference: #41726

@agiudiceandrea
Copy link
Contributor

@agiudiceandrea agiudiceandrea commented Apr 23, 2021

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

@gioman gioman added the Feedback label Apr 23, 2021
@nyalldawson
Copy link
Collaborator

@nyalldawson nyalldawson commented Apr 23, 2021

Yes, it's resolved!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Feature Request Feedback
Projects
None yet
Development

No branches or pull requests

7 participants