PyQGIS – Function to Calculate Percentage of Each Feature Over the Sum of That Feature’s Field

pyqgisqgis-3

I would like to calculate the percentage of each feature (i.e. each row in "HISTO_1") over the sum of that feature's field ("HISTO_1") as a new column (i.e. "Percent_HISTO_1") in an attribute table by creating a function using PyQGIS but I can't seems to do so.

I can do it easily using the expression in the field calculator as follow:

Result field:

Here is the function that I tried to create following this guide: https://docs.qgis.org/3.28/en/docs/user_manual/expressions/expression.html#function-editor

from qgis.core import *
from qgis.gui import *

@qgsfunction(args='auto', group='Custom')
def percent_feature(field1, feature, parent):
    
    totalSum = aggregate(QgsAggregateCalculator.Sum, field1)
    
    percent = (field1/totalSum)*int(100)
    
    return percent

I can't get it to work with the above function. Could someone please help me with this 🙂

Best Answer

A couple of problems. Firstly, the aggregate method you are using is method of the QgsVectorLayer class, so you need to call it on a QgsVectorLayer object. Then, the second argument needs to be a field name string, however the field1 object is the value in the field you pass to your custom function, evaluated for each feature.

Try the following:

from qgis.core import *
from qgis.gui import *

@qgsfunction(args='auto', group='Custom', referenced_columns=[])
def feature_percent(field1, field_name, feature, parent, context):
    layer_id = context.variable('layer_id')
    layer = QgsProject.instance().mapLayer(layer_id)
    totalSum = layer.aggregate(QgsAggregateCalculator.Sum, field_name)[0]
    percent = (field1/totalSum)*100
    return int(percent)

If you want decimal places, change the last line e.g.

return round(percent, 2)

The first argument of this custom function is the field (double quoted) which will be evaluated for each feature. The second argument is the name of the field as a string (single quoted) which will be passed to the aggregate method.

Example usage:

Result:

Honestly, I don't really see a huge advantage of a custom function over the field calculator expression. Perhaps you would prefer a small script which you could run in the Python console which will add and fill all your percentage fields in one go.

Try:

layer = QgsProject.instance().mapLayersByName('canopyChangeDistrict_2018-2022')[0]

field_names = ['HISTO_0',
                'HISTO_1',
                'HISTO_2',
                'HISTO_3',
                'HISTO_4',
                'HISTO_5']

atts_map = {}
                
for fname in field_names:
    new_field = QgsField(f'Percent_{fname}', QVariant.Double, len=5, prec=2)
    layer.dataProvider().addAttributes([new_field])
    layer.updateFields()
    
for ft in layer.getFeatures():
    atts = {}
    for fname in field_names:
        fld_idx = layer.fields().lookupField(f'Percent_{fname}')
        total_sum = layer.aggregate(QgsAggregateCalculator.Sum, fname)[0]
        pcnt = (ft[fname]/total_sum)*100
        atts[fld_idx] = pcnt
    atts_map[ft.id()] = atts
    
#print(atts_map)
layer.dataProvider().changeAttributeValues(atts_map)

Test result:

Related Solutions

PyQGIS – Iterating Over Layer’s Features and Rendering PNG for Each Feature

Try moving the setExtent() call on the QgsMapSettings object inside your loop and pass the bounding box of the current feature's geometry on each iteration:

def finished():
    img = render.renderedImage()
    # save the image; img.save("C:\QGIS Project\Restructured\Rev1\\test_render.png","png")
    img.save(f'C:\QGIS Project\Restructured\Rev1\\test_render_{feature[0]}.png', "png")

layer = QgsProject.instance().mapLayersByName("Rotated")[0]
settings = QgsMapSettings()
settings.setLayers([layer])
settings.setBackgroundColor(QColor(255, 255, 255))
settings.setOutputSize(QSize(400, 400))

for feature in layer.getFeatures():
    layer.setSubsetString("fid = {}".format(feature.id()))
    settings.setExtent(feature.geometry().boundingBox())
    render = QgsMapRendererParallelJob(settings)
    render.finished.connect(finished)
    # Start the rendering
    render.start()
    render.waitForFinished()

layer.setSubsetString("")

By the way, if you wish to add a small margin around each feature's extent you can use the grow() method on the QgsRectangle object returned by feature.geometry().boundingBox(). So your loop block could look like:

for feature in layer.getFeatures():
    layer.setSubsetString("fid = {}".format(feature.id()))
    extent = feature.geometry().boundingBox()
    extent.grow(0.01) # Play with the value to find a suitable margin
    settings.setExtent(extent)
    render = QgsMapRendererParallelJob(settings)
    render.finished.connect(finished)
    # Start the rendering
    render.start()
    render.waitForFinished()

PyQGIS – Calculate Area in hA by Creating Multiple Fields

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)

Best Answer

Related Solutions

PyQGIS – Iterating Over Layer’s Features and Rendering PNG for Each Feature

PyQGIS – Calculate Area in hA by Creating Multiple Fields

Related Question