3

I use the following code to calculate the areas of a multi-polygon shapefile. However, to validate these results, I calculate the polygon’ areas directly in table attribute by creating a new field and applying the formula: $area. I found discrepancies in calculated areas by both methods.

Does the QgsDistanceArea class consider the layer CRS (in my case I use EPSG:32632) and $area does not, that means the geographic WSG84 EPSG:4326 is considered by default?

def area_poly(layer):
  L = []
  features = layer.getFeatures()
  d = QgsDistanceArea()
  d.setSourceCrs(layer.crs(), QgsProject.instance().transformContext())    
  d.setEllipsoid(QgsProject.instance().ellipsoid())
  for i, element in enumerate(features):
    for poly in element.geometry().asMultiPolygon():
        m = d.measurePolygon(poly[0])
        L.append((i, round(m,1))
  return L

enter image description here

2
  • It would be helpful if you'd show some output results. We don't know if you're talking about a 0.1% difference or a 100% difference. If it's a large difference, it is most likely due to the CRS interpretation.
    – Jon
    Jul 28, 2018 at 20:18
  • The results are shown in the attached figure. The differences are not very significant 0.2%. Thank you for more clarifications.
    – ennine
    Jul 28, 2018 at 20:54

1 Answer 1

6

It's the opposite - $area will consider both the layer CRS and project ellipsoid setting. But your code does not, as creating a default constructed QgsDistanceArea is not sufficient to get ellipsoid based measurements, and will currently be returning euclidean areas only.

You should revise your code to (for QGIS 3.x):

d = QgsDistanceArea()
d.setSourceCrs(layer.crs(), QgsProject.instance().transformContext())
d.setEllipsoid(QgsProject.instance().ellipsoid())
3
  • Many thanks ndawson for your clarifications, it works perfectly. However, I am confused for another issue: when I save the line layer (EPSG 32632) in another CSR (for example EPSG: 22391), the object’s length changes considerably: L_32632=68.7m and L_22391=82.5m. Why couldn't I have the same length? Many thanks in advance for more clarifications.
    – ennine
    Jul 29, 2018 at 8:06
  • How are you measuring this? If it's euclidean (non-ellipsoidal) distances then different measurements are expected across different projections. If you're using ellipsoidal measurements, double check that the same ellipsoid is being used in both cases.
    – ndawson
    Jul 30, 2018 at 3:29
  • I apologize Ndawson, I did a mistake while constructing the layer. Now, it works perfectly. Many thanks
    – ennine
    Jul 30, 2018 at 12:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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