Coordinate reference systems and projections are at the heart of any GIS. Unless something went very wrong you usually don't see much of it, once you have settled for one system. However, I stumbled upon a CRS-related issue in QGIS, which is quite hard to spot on first glance and for which I currently don't find a solution.
1) If you use the Field Calculator to calculate certain properties of your geometries, you may have noticed, that you have two functions for most calculations, e.g. $area and area($geometry). On first glance the results appear to be the same, but if you compare the results of two calculations directly you will notice a difference, especially for larger polygons. For small features differ only by a few square centimeters. The same is the case for $length and length($geometry) etc.
The issue becomes clear, if you look into the description of the two functions:
$area: "Returns the area of the current feature. The area calculated by this function respects both the current project's ellipsoid setting and area unit settings. For example, if an ellipsoid has been set for the project then the calculated area will be ellipsoidal, and if no ellipsoid is set then the calculated area will be planimetric."
area: "Returns the area of a geometry polygon object. Calculations are always planimetric in the Spatial Reference System (SRS) of this geometry, and the units of the returned area will match the units for the SRS. This differs from the calculations performed by the $area function, which will perform ellipsoidal calculations based on the project's ellipsoid and area unit settings."
This isn't a problem per se, since it only means that you have to choose the right function for your calculation. For most problems $area should be the right option. Note that the results differ according to the CRS you have set (compare the results in the image above (EPSG 3857) to the ones below (EPSG 25833)).
2) However, this shows clearly how much the ellipsoidal and the planimetric calculations differ. This is an issue, if you want to construct features using the Advanced Digitizing Toolbar > CAD-Tool. As you may know, this tool does not work for geographic coordinates, so you will have to use a projected coordinate system. This would not be a problem if the CAD-Tool would be able to use ellipsoidal distances, but this is not the case.
The polygon in this image has been constructed using the CAD-Tool and with a fixed distance of exactly 100 m and a fixed angle of 90° in EPSG 25833 (note that the results differ a lot if you use e.g. EPSG 3857). The results of the calculations show that the construction was done planimetric and therefore differs from the real result on the curved surface of the earth.
Now my questions:
1) Is there a way to construct features in QGIS that take the ellipsoid into acount?
2) If not, which of the results can be taken to be "correct"?
3) Is the difference of in this case 0.05% a real problem for the daily archaeological practice?
3 Kommentare
Benjamin Ducke sagt:
Hi Lukas, your example above is very extreme. You use EPSG 3857 (Web Mercator), which is a messed up system that pretends to be "WGS84" but in fact disregards the WGS84 ellipsoid (it scales both geographic axes equally, which is great for 3D toys such as GoogleEarth but results in wrong distance and area calculations compared to a true geodetic system), and you compare to EPSG 25833, which is ETRS89-based and uses the GRS80 ellipsoid (slightly different from the WGS84 ellipsoid). So you have two sources of error here which combined could easily account for the difference in area calculations that you see. Try staying with "compatible" systems that use the same WGS84 datum and ellipsoid (e.g. modern UTM with WGS84 for your area) and see how that works out.
Lukas Goldmann sagt: Autor
Hi Benjamin, thanks for the hint. Of course Pseudo-Mercator is totally unsuitable for tasks like this, however, the effect is still there even if you use different systems. I just tested it again with a square of 1000 m side length constructed in WGS84/UTM. Result: $area = 1000530,603 m², area($geometry) = 1000000,00 m².
Benjamin Ducke sagt:
I think this could be one for the QGIS mailing list. I vaguely remember a discussion there, about the different area calculation modes, a few months ago.