1

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)[0]

    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[0])
            else:
                polygon = geom.asPolygon()
                polygon_area = area.measurePolygon(polygon[0])

            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?

3

If I´m not mistaken, you will need to provide a transformation context and set the layers source CRS to the QgsDistanceArea() object accordingly. Try:

...

def calculate_area(in_lyr_name, ellipsoid, units):
    input_layer = QgsProject.instance().mapLayersByName(input_lyr_name)[0]
    lyr_crs = input_layer.crs()

    elps_crs = QgsCoordinateReferenceSystem()
    elps_crs.createFromUserInput(ellipsoid)

    trans_context = QgsCoordinateTransformContext()
    trans_context.calculateDatumTransforms(lyr_crs, elps_crs)

    area = QgsDistanceArea()
    area.setEllipsoid(ellipsoid)
    area.setSourceCrs(lyr_crs, trans_context)

...

Seems to work for me (uses ellipsoid, returns sqm), but pyQGIS is not my expertise.

  • Thank you! Exact method name in QGIS3 turned to be .createFromUserInput() – Basile Jun 7 '18 at 19:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.