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

Nicolas Cadieux njacadieux.gitlab at gmail.com
Tue Jan 5 10:04:13 PST 2021


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 
> <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 
> <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
> 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/0aee50d7/attachment.html>


More information about the Qgis-user mailing list