qgis 3 - Calculate ellipsoidal area for a projected layer in PyQGIS 3

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?

Answer

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.