[Qgis-user] Confusion with ellipsoidal method of $area
Nicolas Cadieux
njacadieux.gitlab at gmail.com
Tue Jan 5 12:52:27 PST 2021
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
> <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
> <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 <mailto:Qgis-user at lists.osgeo.org>
>
> List info:https://lists.osgeo.org/mailman/listinfo/qgis-user <https://lists.osgeo.org/mailman/listinfo/qgis-user>
>
> Unsubscribe:https://lists.osgeo.org/mailman/listinfo/qgis-user <https://lists.osgeo.org/mailman/listinfo/qgis-user>
>
> --
> Nicolas Cadieux
> https://gitlab.com/njacadieux <https://gitlab.com/njacadieux>
--
Nicolas Cadieux
https://gitlab.com/njacadieux
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/qgis-user/attachments/20210105/64c25ca5/attachment.html>
More information about the Qgis-user
mailing list