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 <[email protected]> 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 < > [email protected]> 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 <[email protected]> >> <[email protected]> >> *Lähetetty:* tiistai 5. tammikuuta 2021 20.04 >> *Vastaanottaja:* Rahkonen Jukka (MML) >> <[email protected]> >> <[email protected]>; [email protected] >> *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 >> >> [email protected] >> >> 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 >> [email protected] >> List info: https://lists.osgeo.org/mailman/listinfo/qgis-user >> Unsubscribe: https://lists.osgeo.org/mailman/listinfo/qgis-user >> >
_______________________________________________ Qgis-user mailing list [email protected] List info: https://lists.osgeo.org/mailman/listinfo/qgis-user Unsubscribe: https://lists.osgeo.org/mailman/listinfo/qgis-user
