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

Sebastian Gutwein bas at rdgland.com
Wed Jan 6 08:11:16 PST 2021


Hello,
Sorry if this gets posted twice I previously attached a .gpkg and it
bounced due to the 100kb limit.
I am enjoying this fascinating discussion. This is all fairly new to me so
I looked some things up and want to make sure that I am getting things
right,

   1. In QGIS the $area function uses the sphere/ellipsoid settings of the
   project not of the CRS of the layer.
   2. QGIS uses the PROJ library to calculate area
   3. According to the PROJ documentation
   even.rouault.free.fr/proj_cpp_api/rst_generated/html/geodesic.html PROJ
   uses GeographicLib to calculate area.

If this is all true then if as  Jukka Rahkonen describes the layer and the
project were both EPSG:4326 and the ellipsoid in the Project-Settings was
WGS 84 (EPSG:7030) then the QGIS measurements should be consistent (if they
are correct is open to discussion).with this GeographicLib web
calculator that says it uses the WGS84 ellipsoid
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

Using Jukka Rahkonen's points and processes I ended up with the same non
consistent results as they reported. Attached is a geopackage with the
polygon that I used,

On Wed, Jan 6, 2021 at 10:29 AM Sebastian Gutwein <bas at rdgland.com> wrote:

> Hello,
> I am enjoying this fascinating discussion. This is all fairly new to me so
> I looked some things up and want to make sure that I am getting things
> right,
>
>    1. In QGIS the $area function uses the sphere/ellipsoid settings of
>    the project not of the CRS of the layer.
>    2. QGIS uses the PROJ library to calculate area
>    3. According to the PROJ documentation
>    even.rouault.free.fr/proj_cpp_api/rst_generated/html/geodesic.html
>    PROJ uses GeographicLib to calculate area.
>
> If this is all true then if as  Jukka Rahkonen describes the layer and the
> project were both EPSG:4326 and the ellipsoid in the Project-Settings was
> WGS 84 (EPSG:7030) then the QGIS measurements should be consistent (if they
> are correct is open to discussion).with this GeographicLib web
> calculator that says it uses the WGS84 ellipsoid
> 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
>
> Using Jukka Rahkonen's points and processes I ended up with the same non
> consistent results as they reported. Attached is a geopackage with the
> polygon that I used,
> -Sebastian
>
> On Tue, Jan 5, 2021 at 3:52 PM Nicolas Cadieux <
> njacadieux.gitlab at gmail.com> wrote:
>
>> Hi,
>>
>> QGIS is currently built on Proj version 6.3.2-1.  If the other libraries
>> are using a spheroid by default, then they use a sphere for speed and not
>> an ellipsoid for precision.  You can probably force this measurement in
>> QGIS by creating a custom CRS with a spheroid rather than an ellipsoid. The
>> other option is to use the python in QGIS and to force a geoid.
>>
>> Nicolas
>> On 2021-01-05 1:33 p.m., Rahkonen Jukka (MML) wrote:
>>
>> 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>
>> <njacadieux.gitlab at gmail.com>
>> *Lähetetty:* tiistai 5. tammikuuta 2021 20.04
>> *Vastaanottaja:* Rahkonen Jukka (MML)
>> <jukka.rahkonen at maanmittauslaitos.fi>
>> <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
>>
>> 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
>>
>> --
>> Nicolas Cadieuxhttps://gitlab.com/njacadieux
>>
>> _______________________________________________
>> 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
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/qgis-user/attachments/20210106/590c3557/attachment-0001.html>


More information about the Qgis-user mailing list