Hi, Sorry, I am amateur on this area and I do not know what is the difference between spheroids and ellipsoids. Do you mean that QGIS is using more accurate method than PostGIS for calculating the lengths and areas over the WGS 84 ellipsoid?
The article that is used as a reference on the GeographicLib page (Charles F. F. Karney, Algorithms for geodesics, J. Geodesy 87(1), 43-55 (Jan. 2013) does not use word spheroid in the text at all but ellipsoid appears there 29 times. So perhaps the library deals actually with ellipsoid but PostGIS and Spatialite documentation talks about spheroid? Could you give one linestring and one polygon as reference geometries and the most accurate lengths and areas that you know for some of your favorite ellipsoids to be used as true values? -Jukka Rahkonen- Lähettäjä: Nicolas Cadieux <njacadieux.git...@gmail.com> Lähetetty: tiistai 5. tammikuuta 2021 22.52 Vastaanottaja: Rahkonen Jukka (MML) <jukka.rahko...@maanmittauslaitos.fi>; qgis-user@lists.osgeo.org Aihe: Re: [Qgis-user] Confusion with ellipsoidal method of $area 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.git...@gmail.com><mailto:njacadieux.git...@gmail.com> Lähetetty: tiistai 5. tammikuuta 2021 20.04 Vastaanottaja: Rahkonen Jukka (MML) <jukka.rahko...@maanmittauslaitos.fi><mailto:jukka.rahko...@maanmittauslaitos.fi>; qgis-user@lists.osgeo.org<mailto:qgis-user@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@lists.osgeo.org<mailto:Qgis-user@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 Cadieux https://gitlab.com/njacadieux
_______________________________________________ Qgis-user mailing list Qgis-user@lists.osgeo.org List info: https://lists.osgeo.org/mailman/listinfo/qgis-user Unsubscribe: https://lists.osgeo.org/mailman/listinfo/qgis-user
