PyQGIS – Calculate Area in hA by Creating Multiple Fields

pyqgisqgis-3

I'm writing code that add some attributes. I need create an attribute, 'AREA_HA', and to calculate the area in HA. My code works well, but doesn't create the attribute 'AREA_HA' and calculation produces no result. I looked at some answers here in stack, but none applied to the structure of my code. How I can improve it?

path = 'D:/shps/SHP'
for root, directory, files in os.walk(path):
    for file in files:
        if file.endswith('.shp'):
            full_path = os.path.join(root,file)
            vlayer = QgsVectorLayer(full_path,file[:-4],'ogr')
            feats = vlayer.getFeatures()
            for feat in feats:
                geom = feat.geometry()
                if geom.type() == QgsWkbTypes.PolygonGeometry:
                    vlayer.startEditing()
                    vlayer.addAttribute(QgsField('ID',QVariant.String, 'string',50))
                    vlayer.addAttribute(QgsField('STUDY',QVariant.String,'string',100))
                    vlayer.addAttribute(QgsField('PROJECT',QVariant.String,'string',100))
                    vlayer.addAttribute(QgsField('STAGE',QVariant.String,'string',100))
                    vlayer.addAttribute(QgsField('DATE',QVariant.String,'string',30))
                    vlayer.addAttribute(QgsField('COMPANY',QVariant.String,'string',50))
                    expression = '($area/10000)'
                    field = QgsField('AREA_HA', QVariant.Double, 'double',19)
                    vlayer.addExpressionField(expression, field)
                    vlayer.commitChanges()
                print(vlayer)

Best Answer

You have a number of issues here. Firstly, with the addExpressionField() method, you are trying to add a virtual field to a QgsVectorLayer object which is not loaded into a project, which simply won't work. Virtual fields are added to the end of an attribute table in a project vector layer and exist only in the project settings in which they have been created. They cannot be permanently added to the layer's data provider. Please see the first point under the Note section here.

Other than that you have some logic errors e.g. you are adding fields to the layer for every feature, so you should move that section out of the feature iterator loop. Also you need to call updateFields() after adding the fields.

I strongly recommend studying the PyQGIS Developer Cookbook, especially the following sections:

Modifying Vector Layers

Expressions with features

Try the following script:

path = 'D:\\shps\\SHP'

for root, directory, files in os.walk(path):
    for file in files:
        if file.endswith('.shp'):
            full_path = os.path.join(root,file)
            vlayer = QgsVectorLayer(full_path,file[:-4],'ogr')
            if not vlayer.isValid() or not vlayer.geometryType() == QgsWkbTypes.PolygonGeometry:
                if not vlayer.isValid():
                    print(f'Layer from {file} not valid')
                if vlayer.geometryType() != QgsWkbTypes.PolygonGeometry:
                    print(f'Layer from {file} not a polygon layer')
                continue
            vlayer.dataProvider().addAttributes([
                QgsField('ID', QVariant.String, len=50),
                QgsField('STUDY', QVariant.String, len=100),
                QgsField('PROJECT', QVariant.String, len=100),
                QgsField('STAGE', QVariant.String, len=100),
                QgsField('DATE', QVariant.String, len=30),
                QgsField('COMPANY', QVariant.String, len=50),
                QgsField('AREA_HA', QVariant.Double, len=19, prec=9)
            ])
            vlayer.updateFields()# Need this line!
            
            e = QgsExpression('$area/10000')
            c = QgsExpressionContext()
            c.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(vlayer))

            att_map = {}

            for ft in vlayer.getFeatures():
                c.setFeature(ft)
                att_map[ft.id()] = {vlayer.fields().lookupField('AREA_HA'): e.evaluate(c)}
                
            vlayer.dataProvider().changeAttributeValues(att_map)
Related Question