I'd like to loop though each feature in a layer and calculate its ellipsoidal area using QgsDistanceArea() in PyQGIS 3.
Everything works fine on unprojected layers, but I get strange results when input layer is in projected CRS.
So far I have written this:
from qgis.core import * def calculate_area(in_lyr_name, ellipsoid, units): input_layer = QgsProject.instance().mapLayersByName(input_lyr_name) area = QgsDistanceArea() area.setEllipsoid(ellipsoid) features = input_layer.getFeatures() feature_area = 0 error_features_counter = 0 for i, feat in enumerate(features): geom = feat.geometry() polygon_area = 0 try: if geom.isMultipart(): polygons = geom.asMultiPolygon() for polygon in polygons: polygon_area += area.measurePolygon(polygon) else: polygon = geom.asPolygon() polygon_area = area.measurePolygon(polygon) feature_area += polygon_area except AttributeError: # catches NoneType geometries (broken etc.) error_features_counter += 1 pass # calculated area is in sq. metres (see the "else" case) if units == 'km²': final_area = feature_area / 1e6 elif units == 'Ha': final_area = feature_area / 10000 else: final_area = feature_area return final area, error_features_counter
How do I get a constant result that does not depend on the presence of projection?