While most GIS software supports 3D editing capabilities, calculating the surface area of polygons in 3D remains a serious challenge. This issue arose when we were writing the roof analyzer script since it requires QGIS 3D area calculation, namely the calculation of polygon surfaces.
In this article, we present the problem of polygon QGIS 3D area calculation, the script that provides the solution, and its implementation in QGIS.
We present the polygon 3D area calculation in QGIS because it is open-source, making it accessible to anyone. The method we demonstrate can also be used in ArcGIS, but arcpy functions must be applied there to extract geometric data.
Sidenote: QGIS 3D area calculation can also be replaced with Blender Computer Graphics software. In Blender, there is no need to write code, but you need to install the Blender GIS and 3D-Print Toolbox addons if you want to calculate the area of a shapefile.
Basic problem
GIS software was originally designed for map visualization and map-based analysis. Joining the 3D world came later in this technology. The calculation of the area of geometries with map projections is done confidently and with high precision.
In general, polygon 2D area calculation can be done using Heron’s formula. GIS software complements this with the characteristics of the coordinate reference system (CRS). Let’s see what built-in QGIS area calculations are available!
- $area – Returns the area of the current feature. The area calculated with this function takes into account the ellipsoid settings and the unit settings of the area in the current project. For example, if an ellipsoid is set in the current project, the calculated area will be ellipsoidal, and if no ellipsoid is set, the calculated area will be planar.
- area – Returns the area of the polygon object. Calculations are always planar in the geometry coordinate reference system (CRS), and the units of the returned area match the units of the CRS. This differs from the calculation performed by the $area function, which is capable of ellipsoidal calculations based on the ellipsoid and area unit settings of the current project.
So whichever QGIS area calculation we use, in both cases we get the area of the polygon in two dimensions, even though we created a 3D polygon (PolygonZ). Calculating the area of a polygon in QGIS in 3D is a more complex operation and unfortunately not part of QGIS’s built-in functions.
Python 3D polygon area calculation with Stokes’ theorem
To properly apply QGIS 3D area calculation, we reviewed numerous forums, scientific articles, and GitHub repositories until we found the python code by caougolin. He used Stokes’ theorem to write a general python code for 3D polygon area calculation. Using caougolin’s solution, we wrote Python code for QGIS 3D area calculation. Below, we present the original python code followed by the version we used.
Using Stokes’ function, caougolin wrote the following Python code for 3D polygon area calculation:
import math
def compute_3D_polygon_area(points):
if (len(points) < 3): return 0.0
P1X,P1Y,P1Z = points[0][0],points[0][1],points[0][2]
P2X,P2Y,P2Z = points[1][0],points[1][1],points[1][2]
P3X,P3Y,P3Z = points[2][0],points[2][1],points[2][2]
a = pow(((P2Y-P1Y)*(P3Z-P1Z)-(P3Y-P1Y)*(P2Z-P1Z)),2) + pow(((P3X-P1X)*(P2Z-P1Z)-(P2X-P1X)*(P3Z-P1Z)),2) + pow(((P2X-P1X)*(P3Y-P1Y)-(P3X-P1X)*(P2Y-P1Y)),2)
cosnx = ((P2Y-P1Y)*(P3Z-P1Z)-(P3Y-P1Y)*(P2Z-P1Z))/(pow(a,1/2))
cosny = ((P3X-P1X)*(P2Z-P1Z)-(P2X-P1X)*(P3Z-P1Z))/(pow(a,1/2))
cosnz = ((P2X-P1X)*(P3Y-P1Y)-(P3X-P1X)*(P2Y-P1Y))/(pow(a,1/2))
s = cosnz*((points[-1][0])*(P1Y)-(P1X)*(points[-1][1])) + cosnx*((points[-1][1])*(P1Z)-(P1Y)*(points[-1][2])) + cosny*((points[-1][2])*(P1X)-(P1Z)*(points[-1][0]))
for i in range(len(points)-1):
p1 = points[i]
p2 = points[i+1]
ss = cosnz*((p1[0])*(p2[1])-(p2[0])*(p1[1])) + cosnx*((p1[1])*(p2[2])-(p2[1])*(p1[2])) + cosny*((p1[2])*(p2[0])-(p2[2])*(p1[0]))
s += ss
s = abs(s/2.0)
return s
polygon = [[0,0,0],[0,1,0],[1,1,1],[1,0,1]]
result = compute_3D_polygon_area(polygon)
Practical QGIS 3D area calculation
I will now demonstrate the practical aspect of QGIS 3D area calculation, i.e., how to use the Python code based on Stokes’ theorem. Launch QGIS and open the Python console (Ctrl+Alt+P).
Note: The Python code we wrote is suitable for 3D area calculation of PolygonZ type geometries. Furthermore, it can handle more complex polygons with rings!
Click the “Show Editor” button. You can write Python code in the right-hand side section.
Copy and paste the following code. Make sure to replace ‘PATH’ with the file path of the polygon for which you want to calculate the 3D area in the “layer = QgsVectorLayer(‘PATH’)” line.
from qgis.core import *
from PyQt5.QtCore import *
import math
def compute_3D_polygon_area(points):
if (len(points) < 3): return 0.0
P1X,P1Y,P1Z = points[0][0],points[0][1],points[0][2]
P2X,P2Y,P2Z = points[1][0],points[1][1],points[1][2]
P3X,P3Y,P3Z = points[2][0],points[2][1],points[2][2]
a = pow(((P2Y-P1Y)*(P3Z-P1Z)-(P3Y-P1Y)*(P2Z-P1Z)),2) + pow(((P3X-P1X)*(P2Z-P1Z)-(P2X-P1X)*(P3Z-P1Z)),2) + pow(((P2X-P1X)*(P3Y-P1Y)-(P3X-P1X)*(P2Y-P1Y)),2)
cosnx = ((P2Y-P1Y)*(P3Z-P1Z)-(P3Y-P1Y)*(P2Z-P1Z))/(pow(a,1/2))
cosny = ((P3X-P1X)*(P2Z-P1Z)-(P2X-P1X)*(P3Z-P1Z))/(pow(a,1/2))
cosnz = ((P2X-P1X)*(P3Y-P1Y)-(P3X-P1X)*(P2Y-P1Y))/(pow(a,1/2))
s = cosnz*((points[-1][0])*(P1Y)-(P1X)*(points[-1][1])) + cosnx*((points[-1][1])*(P1Z)-(P1Y)*(points[-1][2])) + cosny*((points[-1][2])*(P1X)-(P1Z)*(points[-1][0]))
for i in range(len(points)-1):
p1 = points[i]
p2 = points[i+1]
ss = cosnz*((p1[0])*(p2[1])-(p2[0])*(p1[1])) + cosnx*((p1[1])*(p2[2])-(p2[1])*(p1[2])) + cosny*((p1[2])*(p2[0])-(p2[2])*(p1[0]))
s += ss
s = abs(s/2.0)
return s
# Set the PolygonZ/MultiPolygonZ file path
layer = QgsVectorLayer('PATH')
area3d = layer.dataProvider()
area3d.addAttributes( [ QgsField('Area3D', QVariant.Double, "double", 10, 2)] )
layer.updateFields()
feat = layer.getFeatures()
for feature in feat:
points = []
for v in feature.geometry().vertices():
vertex = [v.x(), v.y(), v.z()]
points.append(vertex)
area = compute_3D_polygon_area(points)
area = round(area, 2)
print(area)
layer.startEditing()
feature['Area3D'] = area
layer.updateFeature(feature)
layer.commitChanges()
QgsProject.instance().reloadAllLayers()
print("Done!")
You can run the code using the “Run script” green triangle button.
The result will appear in the polygon attribute table. It’s important to note that 3D area calculation only works if there isn’t already a column named “Area3D” in the attribute table!
Limitation of QGIS 3D polygon area calculation
It’s important to note that QGIS 3D area calculation also has its limitations. The 3D area of a polygon can only be calculated if its vertices are ordered. Ordering means that the polygon’s breakpoints can start from any point, but they must either be in clockwise or counter-clockwise order. A typical vectorization follows this precisely, but it’s worth mentioning separately. If the vertices are not provided in order, the result of 3D area calculation will be 0.
Furthermore, QGIS 3D polygon calculation cannot be applied to 2D polygons. For 2D polygons, the result of 3D area calculation will be “Not a Number” (NaN). It doesn’t provide accurate results for MultiPolygonZ geometries consisting of more than one part.
Closing thoughts
During performing operations in QGIS, it can be useful in many cases to calculate the 3D area of a given polygon. SurveyTransfer also provides 3D area calculation, which works similarly to the method described here.
The SurveyTransfer team developed the QGIS 3D area calculation script for the roof analyzer solution.
If you really liked what you read than you can share it with your friends. 🙂
Did you like what you read? Do you want to read similar ones?
