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
[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
_______________________________________________
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
