1

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

  • 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 '18 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 '18 at 20:54
2

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())
  • 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 '18 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 '18 at 3:29
  • I apologize Ndawson, I did a mistake while constructing the layer. Now, it works perfectly. Many thanks – ennine Jul 30 '18 at 12:13

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.