[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