4

I have a shapefile with about 2000 polygons. I am using QGIS and I need to write a standalone Python script which opens the shapefile and calculates the area of all polygons separately and outputs the calculated area in a column "Area" in square kilometers. I tried to do it using the example in Calculate area in square meters from degrees?

eligibleLayer = QgsVectorLayer("/home/usr/Desktop/eligible_areas.shp", "eligible_areas", "ogr")

if not eligibleLayer.isValid():
    print "Eligible areas layer failed to load!" 
lProvider = eligibleLayer.dataProvider()
lProvider.addAttributes( [ QgsField("Area",QVariant.Double) ] )
area = 0 
for gFeat in eligibleLayer.getFeatures():
    calculator = QgsDistanceArea()
    calculator.setEllipsoid('WGS84')
    calculator.setEllipsoidalMode(True)
    geom = gFeat.geometry()
    landArea = gFeat.attributes()[gProvider.fieldNameIndex('Area')].toDouble()[0]
    if geom.isMultipart() is False: # if only simple polygon, calculate only for this
        polyg = geom.asPolygon() #transform to list of points
        if len(polyg) > 0:
            area = calculator.measurePolygon(polyg[0])

    else: #is Multipart
        multi = geom.asMultiPolygon()
        for polyg in multi:
            area = area + calculator.measurePolygon(polyg[0])

The script calculates the variable area but I am unsure if it is doing the right thing since I checked with the geometry tools and I get something completely different. Furthermore, I am unsure about how to update the above code such that I save the information in a newly created attribute column called "Area".

1

You should use the QgsDistanceArea().computeAreaInit() before calculate any surfaces. It's mentionned in the second answer of the post you've mentioned When you update an field you can use the edit() class of qgis.core to edit your field area and use a with (context manager of python) to edit so the code looks like this :

from qgis.core import *

eligibleLayer = QgsVectorLayer("/home/usr/Desktop/eligible_areas.shp", "eligible_areas", "ogr")

if not eligibleLayer.isValid():
    print "Eligible areas layer failed to load!" 

lProvider = eligibleLayer.dataProvider()
lProvider.addAttributes( [ QgsField("Area",QVariant.Double) ] )

area = 0
with edit(eligibleLayer):
    for gFeat in eligibleLayer.getFeatures():
        calculator = QgsDistanceArea()
        calculator.setEllipsoid('WGS84')
        calculator.setEllipsoidalMode(True)

        # here you need to put the code
        calculator.computeAreaInit()

        geom = gFeat.geometry()
        landArea = gFeat['Area'] # it's so simple to call the attribute with edit()
            if geom.isMultipart() is False: # if only simple polygon, calculate only for this
                polyg = geom.asPolygon() #transform to list of points
                if len(polyg) > 0:
                    area = calculator.measurePolygon(polyg[0])
                    landArea = area
        # this section can be simple if you perform the alg single to multiple from the processing toolbox.
        # Doing so you will not have any multiple Polygon 
        else: #is Multipart
            multi = geom.asMultiPolygon()
            for polyg in multi:
                area = area + calculator.measurePolygon(polyg[0])
                landArea = area
  • Thank you for your answer. Unfortunately, when I add the line calculator.computeAreaInit() , I get an error 'QgsDistanceArea' object has no attribute 'computeAreaInit' ...I have no idea why this is happening – ELD Feb 15 '16 at 9:35
  • Which version of Qgis you are using ? – SIGIS Feb 15 '16 at 12:29
  • I am using QGIS 2.0.1 – ELD Feb 15 '16 at 13:48
  • I suggest you to switch to Qgis 2.8 which is a long term revised realease of Qgis (and the first LTR one). Because using the Qgis 2.0.1 you are using an old API version of PyQGIS. The new version implemented the computeAreaInit method of QgsDistanceArea class. Also I recommended you to use osgeo intall if on windows because it's more customizable. – SIGIS Feb 15 '16 at 14:00
  • @SIGIS can you help me on this gis.stackexchange.com/questions/323006/… – ps1 May 16 at 15:25

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.