Cannot calculate new values for an empty attribute using existing attributes in QGIS with python

by Oski Turanoglu   Last Updated January 13, 2018 22:22 PM

I'm having difficulty in calculating the values for features in an empty attribute column by using existing values from other attribute columns in a joined table. For example, the feature 1 in attribute column 1 + feature 1 in attribute column 2 = new feature 1 value in attribute column 3. Specifically, that the column names being used are hard coded and that is definitely not ideal (see code below).

expression = QgsExpression("to_int(\"Primary\" * \"POP06\")")
expression.prepare(huff_model.pendingFields())

for feature in huff_model.getFeatures():
    value = expression.evaluate(feature)
    feature[idx_pri_pop] = value  #idx_pri_pop is the index of the column where the new values go
    huff_model.updateFeature(feature)

"Primary" and "POP06" fields are both hard coded which will cause the script to fail if the user changes the field names or uses different data. I want the code to accept a field name no matter what it is, so it's flexible. However, it appears that QgsExpression does not accept variables, even if they hold strings (or it can and I can't seem to make it work). As such, I've also tried using the changeAttributeValue method, however that does not work at all, as the new field only has null values.

  • I have looked at other answer with joins, however I could not find one that explicitly covers how to calculate a value for a feature using features in other fields as variables for the calculation.
  • The attribute table that I am using is from the target layer after a join. I'm not sure if that might cause any issues or not.
  • It might also be the case that the rest of my code (see below) may not be in the right order or some other parts of it might be incorrect, despite running with no errors. That could be why the changeAttributeValue method might not be working. I'm not entirely sure.

FULL CODE

Processing Script Parameters and Imports

##Huff_Model_Layer=vector
##Mall=field Huff_Model_Layer
##Census_Layer=vector
##Population_Data=field Census_Layer
##Mall_Layer_Name= string test
##PDF_file=output file pdf

from qgis.core import *
from qgis.gui import *
from qgis.utils import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
import os

Creating the layer

# Load the huff model layer
ori_huff_model = processing.getObject(Huff_Model_Layer)
# Load the census layer
census_layer = processing.getObject(Census_Layer)

# If script was run previous, remove the previously created layer to keep 
# layer panel with only one newly created layer at a time
try:
    last_layer = QgsMapLayerRegistry.instance().mapLayersByName(Mall_Layer_Name)[0]
    QgsMapLayerRegistry.instance().removeMapLayers([last_layer])
    del last_layer
except IndexError:
    child1 = root.children()[1]
    if child1.name() == census_layer.name():
        ids = root.findLayerIds()
        last_layer = root.findLayer(ids[0])
        print last_layer.name()
        root.removeChildNode(last_layer)

# Make a copy of the huff model layer in memory to work with so that all 
# changes made to copy only keeping original data intact
huff_model = QgsVectorLayer("Polygon?crs=epsg:4326", Mall_Layer_Name, "memory")
huff_model_data = huff_model.dataProvider()

huff_model.startEditing()

attr = ori_huff_model.dataProvider().fields().toList()
huff_model_data.addAttributes(attr)
huff_model.updateFields()

feat = QgsFeature()
for elem in ori_huff_model.getFeatures():
    feat.setGeometry(elem.geometry())
    feat.setAttributes(elem.attributes())
    huff_model.addFeatures([feat])
    huff_model.updateExtents()

# Choose the mall (field in table) to find the primary and secondary market 
# areas 
mall = ori_huff_model.fieldNameIndex(Mall)
for field in huff_model.fields():
    field_id = huff_model.fieldNameIndex(field.name())
    if mall == field.name():
        mall = huff_model.fieldNameIndex(field.name())
        break

col_num = int(mall)

# Create two new fields that will only show the primary and secondary market 
# areas based on what is chosen as primary and secondary market 
# probabilities
primary = QgsField('Primary', QVariant.Double)
huff_model.addAttribute(primary)
index_pri = huff_model.fieldNameIndex('Primary')

secondary = QgsField('Secondary', QVariant.Double)
huff_model.addAttribute(secondary)
index_sec = huff_model.fieldNameIndex('Secondary')

# Will hold population data based on primary and secondary markets 
# multiplied by the population for each census tract.
pri_pop = QgsField('PRI_POP', QVariant.Int)
huff_model.addAttribute(pri_pop)
idx_pri_pop = huff_model.fieldNameIndex('PRI_POP')

sec_pop = QgsField('SEC_POP', QVariant.Int)
huff_model.addAttribute(sec_pop)
idx_sec_pop = huff_model.fieldNameIndex('SEC_POP')

# Find the relevant features for the chosen mall
for feature in huff_model.getFeatures():
    probability = feature.attributes()
    ctuid = feature["CTUID"]
    if not probability[col_num]:
        print ("This mall does not exist")
    else:
        # Primary market values will go into the "Primary" field
        if probability[col_num] >= .6:
            huff_model.changeAttributeValue(feature.id(), index_pri, probability[col_num])
        # Secondary market values will go into the "Secondary" field
        if .4 <= probability[col_num] < .6:
            huff_model.changeAttributeValue(feature.id(), index_sec, probability[col_num])

# Delete all other features (rows) that do not have a probability in the  
# primary and the secondary fields
expr = QgsExpression("\"Primary\" is NULL and \"Secondary\" is NULL")
for f in huff_model.getFeatures(QgsFeatureRequest(expr)):
    huff_model.deleteFeature(f.id())

huff_model.updateExtents()
huff_model.commitChanges()

Setting up the Join

# Add newly created layer to map to prepare for the join
reg = QgsMapLayerRegistry.instance()
reg.addMapLayer(huff_model)

# Set the appropriate layers for the join. Mall_Layer_Name is the name for 
# the huff_model layer 
targetLyr = QgsMapLayerRegistry.instance().mapLayersByName(Mall_Layer_Name)[0]
censusLyr = QgsMapLayerRegistry.instance().mapLayersByName(census_layer.name())[0]

# I probably should move the block of code (not part of join). It moves the 
# census layer to the second position in the layer panel as a conditional to 
# remove the newly created layer on a subsequent run of the script.
mycensusLyr = root.findLayer(censusLyr.id())
censusClone = mycensusLyr.clone()
parent = mycensusLyr.parent()
parent.insertChildNode(1, censusClone)
parent.removeChildNode(mycensusLyr)

# Set properties for the join
targetField = 'CTUID'
inField = 'CTUID'
joinObject = QgsVectorJoinInfo()
joinObject.joinLayerId = censusLyr.id()
joinObject.joinFieldName = inField
joinObject.targetFieldName = targetField
joinObject.memoryCache = True
joinObject.prefix = ''
print(targetLyr.addJoin(joinObject))  # You should get True as response.
targetLyr.addJoin(joinObject)

# The population data used to calculate the the values for the new fields 
# and populate the features(see code above). However, it can't be used in  
# QgsExpression!
ct_pop = Population_Data
ct_pop = "{}\\".format(ct_pop)

huff_model.startEditing()

# Expression used to populate the "PRI_POP" field (see above code). Both 
# fields are hard coded. Primary is fine, since it's created in the code but 
# "POP06" should not be hard coded. 
expression = QgsExpression("to_int(\"Primary\" * \"POP06\")")
expression.prepare(huff_model.pendingFields())

for feature in huff_model.getFeatures():
    value = expression.evaluate(feature)
    feature[idx_pri_pop] = value
    huff_model.updateFeature(feature)

# Expression used to populate the "PRI_POP" field (see above code). Both 
# fields are hard coded. Primary is fine, since it's created in the code but 
# "POP06" should not be hard coded. 
expression2 = QgsExpression("to_int(\"Secondary\" * \"POP06\")")
expression2.prepare(huff_model.pendingFields())

for feature in huff_model.getFeatures():
    value = expression2.evaluate(feature)
    feature[idx_sec_pop] = value
    huff_model.updateFeature(feature)

# Delete all other columns except for the CTUID, Primary, Secondary columns. 
# For some reason, this has to be here AFTER the join, otherwise the newly  
# calculated features don't show up in "PRI_POP" and "SEC_POP.
# It also doesn't delete fields NOT included in the set in the newly joined 
# table. It does however delete the fields in target layer (huff_model)  
# before the join that are not in the set.
fields = []

fieldnames = {'CTUID', 'Primary', 'Secondary', 'PRI_POP', 'SEC_POP'}
for field in huff_model.fields():
    if field.name() not in fieldnames:
        fields.append(huff_model.fieldNameIndex(field.name()))
        # print(fields.append(huff_model.fieldNameIndex(field.name())))

huff_model.deleteAttributes(fields)
huff_model.commitChanges()

Creating the PDF output

# Make sure that the newly created and joined layer is the active layer.
huff_model = iface.activeLayer()

# Create the attribute using the map composer in QGIS and show the output in 
# PDF format.
mapRenderer = iface.mapCanvas().mapRenderer()
c = QgsComposition(mapRenderer)
c.setPlotStyle(QgsComposition.Print)
x, y = 0, 0
w, h = c.paperWidth(), c.paperHeight()
composerMap = QgsComposerMap(c, x, y, w, h)

table = QgsComposerAttributeTable(c)
table.setComposerMap(composerMap)
table.setScale(.85)
table.setVectorLayer(huff_model)
table.setMaximumNumberOfFeatures(huff_model.featureCount())
c.addItem(table)

# Checks for the existence of a file with the same name and removes it.
if os.path.isfile(PDF_file):
    os.remove(PDF_file)

printer = QPrinter()
printer.setOutputFormat(QPrinter.PdfFormat)
printer.setOutputFileName(PDF_file)
printer.setPaperSize(QSizeF(c.paperWidth(), c.paperHeight()), QPrinter.Millimeter)
printer.setFullPage(True)

pdfPainter = QPainter(printer)
paperRectMM = printer.pageRect(QPrinter.Millimeter)
paperRectPixel = printer.pageRect(QPrinter.DevicePixel)
c.render(pdfPainter, paperRectPixel, paperRectMM)
pdfPainter.end()

The image below shows which fields I require for the calculation of the features in the "PRI_POP" and "SEC_POP" fields for clarification purposes. Please let me know if I need to edit anything if I didn't explain a process well!

Joined Attribute Table After Script is Run

enter image description here



Related Questions