[Qgis-user] Confusion with ellipsoidal method of $area

Rahkonen Jukka (MML) jukka.rahkonen at maanmittauslaitos.fi
Tue Jan 5 10:33:24 PST 2021


Hi,

I suppose that PostGIS is using the WGS 84 ellipsoid but I am not sure where from the documentation I could find that information. The ST_Area document https://postgis.net/docs/ST_Area.html says only "For geography types by default area is determined on a spheroid with units in square meters". Same thing with Spatialite, documentation suggests just that it is "the" spheroid http://www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html.

I did not notice this paragraph in the ST_Area document earlier "Enhanced: 2.2.0 - measurement on spheroid performed with GeographicLib for improved accuracy and robustness. Requires Proj >= 4.9.0 to take advantage of the new feature." So no wonder that the web app and PostGIS give the same results because they both use GeographicLib. And SpatiaLite 5.0 is using RTTopo that is a library that is based on LWGeom so close connection in there too.

So perhaps the way to get identical areas from QGIS would be to make it to use GeographicLib as well. I have no idea if it is a realistic approach and worth making a feature request.

-Jukka Rahkonen-

Lähettäjä: Nicolas Cadieux <njacadieux.gitlab at gmail.com>
Lähetetty: tiistai 5. tammikuuta 2021 20.04
Vastaanottaja: Rahkonen Jukka (MML) <jukka.rahkonen at maanmittauslaitos.fi>; qgis-user at lists.osgeo.org
Aihe: Re: [Qgis-user] Confusion with ellipsoidal method of $area


Hi,

Your method in QGIS is sound.  Area is calculated using the wgs84 ellipsoid EPSG 7030. If you reproject to wgs84 zone 35N (I think this is close), area goes from 249566957499.7546m2  to 249566957499.721m2 or a difference of 0.0336 m2.  I don't think densification would change things much.

My question is the following: You know that QGIS uses the WGS84 Ellipsoid.  What Ellipsoid are using used in the other software???  If you don't know, I would proceed until you figure that out.

This is the currently used ellipsoid in the proj database in QGIS 3.16.  You can get this by typing proj -le in the OSFeo4W Shell.

    MERIT a=6378137.0      rf=298.257       MERIT 1983
    SGS85 a=6378136.0      rf=298.257       Soviet Geodetic System 85
    GRS80 a=6378137.0      rf=298.257222101 GRS 1980(IUGG, 1980)
    IAU76 a=6378140.0      rf=298.257       IAU 1976
     airy a=6377563.396    rf=299.3249646   Airy 1830
   APL4.9 a=6378137.0      rf=298.25        Appl. Physics. 1965
    NWL9D a=6378145.0      rf=298.25        Naval Weapons Lab., 1965
 mod_airy a=6377340.189    b=6356034.446    Modified Airy
   andrae a=6377104.43     rf=300.0         Andrae 1876 (Den., Iclnd.)
   danish a=6377019.2563   rf=300.0         Andrae 1876 (Denmark, Iceland)
  aust_SA a=6378160.0      rf=298.25        Australian Natl & S. Amer. 1969
    GRS67 a=6378160.0      rf=298.2471674270 GRS 67(IUGG 1967)
  GSK2011 a=6378136.5      rf=298.2564151   GSK-2011
   bessel a=6377397.155    rf=299.1528128   Bessel 1841
 bess_nam a=6377483.865    rf=299.1528128   Bessel 1841 (Namibia)
   clrk66 a=6378206.4      b=6356583.8      Clarke 1866
   clrk80 a=6378249.145    rf=293.4663      Clarke 1880 mod.
clrk80ign a=6378249.2      rf=293.4660212936269 Clarke 1880 (IGN).
      CPM a=6375738.7      rf=334.29        Comm. des Poids et Mesures 1799
   delmbr a=6376428.       rf=311.5         Delambre 1810 (Belgium)
  engelis a=6378136.05     rf=298.2566      Engelis 1985
  evrst30 a=6377276.345    rf=300.8017      Everest 1830
  evrst48 a=6377304.063    rf=300.8017      Everest 1948
  evrst56 a=6377301.243    rf=300.8017      Everest 1956
  evrst69 a=6377295.664    rf=300.8017      Everest 1969
  evrstSS a=6377298.556    rf=300.8017      Everest (Sabah & Sarawak)
  fschr60 a=6378166.       rf=298.3         Fischer (Mercury Datum) 1960
 fschr60m a=6378155.       rf=298.3         Modified Fischer 1960
  fschr68 a=6378150.       rf=298.3         Fischer 1968
  helmert a=6378200.       rf=298.3         Helmert 1906
    hough a=6378270.0      rf=297.          Hough
     intl a=6378388.0      rf=297.          International 1909 (Hayford)
    krass a=6378245.0      rf=298.3         Krassovsky, 1942
    kaula a=6378163.       rf=298.24        Kaula 1961
    lerch a=6378139.       rf=298.257       Lerch 1979
    mprts a=6397300.       rf=191.          Maupertius 1738
 new_intl a=6378157.5      b=6356772.2      New International 1967
  plessis a=6376523.       b=6355863.       Plessis 1817 (France)
     PZ90 a=6378136.0      rf=298.25784     PZ-90
   SEasia a=6378155.0      b=6356773.3205   Southeast Asia
  walbeck a=6376896.0      b=6355834.8467   Walbeck
    WGS60 a=6378165.0      rf=298.3         WGS 60
    WGS66 a=6378145.0      rf=298.25        WGS 66
    WGS72 a=6378135.0      rf=298.26        WGS 72
    WGS84 a=6378137.0      rf=298.257223563 WGS 84
   sphere a=6370997.0      b=6370997.0      Normal Sphere (r=6370997)

Nicolas


On 2021-01-05 11:43 a.m., Rahkonen Jukka (MML) wrote:
Hi,

I wonder what method QGIS is using when it computes ellipsoidal area with $area function. I made a test with this EPSG:4326 polygon:
POLYGON (( 20.13293641 59.95688345, 26.94617837 60.47397663, 29.74782155 62.56499443, 27.45254202 68.7065034, 23.75771765 68.24937206, 25.42698984 65.27444593, 21.51545237 63.10353609, 21.4056276 61.12318104, 19.41123592 60.40477513, 20.13293641 59.95688345 ))

I checked that the CRS of the project and layer were both EPSG:4326. The ellipsoid in the Project-Settings was WGS 84 (EPSG:7030). The area that $area returns is
249566957499.7546

As a references I used PostGIS and web site https://geographiclib.sourceforge.io/scripts/geod-calc.html#area. With my test polygon they both give this result:
251199344354.4308

PostGIS returns bigger area. The difference is 0.654% so not huge but not negligible either.

A third test with SpatiaLite 5.0 gives a result that is very close to PostGIS (difference 0.001%).
251195856999.549927

As a conclusion I have decided to trust in PostGIS because it gives the same results than the web site https://geographiclib.sourceforge.io/scripts/geod-calc.html that feels scientifically sound. However, I wonder if I have used QGIS in a correct way or if there is anything I could do for getting areas to match better with PostGIS for example by densifying the long segments in my test polygon.

For getting testing easy I attach directly runnable  SQL for PostGIS and Spatialite (version 5.0 with RTTopo is needed) as well as a coordinate list for the web app.

PostGIS
select st_area(st_geogfromtext('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))') , true);

Spatialite
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))') ,true);

Web site (notice lat-lon order)
59.95688345 20.13293641
60.47397663 26.94617837
62.56499443 29.74782155
68.70650340 27.45254202
68.24937206 23.75771765
65.27444593 25.42698984
63.10353609 21.51545237
61.12318104 21.40562760
60.40477513 19.41123592
59.95688345 20.13293641

-Jukka Rahkonen-



_______________________________________________

Qgis-user mailing list

Qgis-user at lists.osgeo.org<mailto:Qgis-user at lists.osgeo.org>

List info: https://lists.osgeo.org/mailman/listinfo/qgis-user

Unsubscribe: https://lists.osgeo.org/mailman/listinfo/qgis-user

--

Nicolas Cadieux

https://gitlab.com/njacadieux
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/qgis-user/attachments/20210105/8c20e5ed/attachment-0001.html>


More information about the Qgis-user mailing list