2

I have written a python script in QGIS which estimates the amount of residents for each building polygon by using census block data and building polygons (about 500 000 buildings in an area of about 900 km2).

The relevant data from the census blocks was merged to the building layer, so only the building layer is used in the script to collect the data. It has the following fields:

| A_building | id_census | pop_census_block |

The steps of the script are:

  1. Sort buildings by id_census
  2. Sum all building footprint areas for each census block
  3. Calculate the areal share and from that the amount of residents for each footprint.
  4. Update 2 attributes of the feature (layer.updateFeature(feat[share, residents]))

The 3rd and 4th step happen in a loop, that iterates over each building polygon. This loop takes two days, which seems to me too long for such a simple operation. Is there a faster way to update the features of a layer, like sending the calculated data as a list?

Here is the relevant section in the script:

# feats is a python list:
feats = [feat for feat in layer.getFeatures(QgsFeatureRequest().setFlags(QgsFeatureRequest.NoGeometry).setSubsetOfAttributes(limit_cols))
]

# sort feats
feats.sort(key=itemgetter(col_census_id))

...

# Get unique census ids
census_ids = []
areas = dict()
for feat in feats:
    census_id = feat[col_census_id]
    if census_id not in census_ids:
        census_ids.append(census_id)
        areas[census_id] = 0

# Sum building area per census block
for i in range(len(feats)):
    # Sum area
    census_id = feats[i][col_census_id]
    areas[census_id] += feats[i][col_bldg_area]

# Start an undo block
layer.beginEditCommand('Calculate share of area')
i, j, k = 0, 0, 0
layer.startEditing()
for feat in feats:
    # Share of area
    share = feat[col_bldg_area] / float(areas[feat[col_census_id]])
    # Add share to feature
    feat[col_share] = share
    # Building population
    feat[col_bldg_pop] = int(feat[col_pop] * feat[col_share])

    layer.updateFeature(feat)

    ...

# End the undo block
layer.endEditCommand()

# Comit changes
layer.commitChanges()
return layer
  • How big of an area are we talkin' here? Summing polygons for a county could take awhile. One approach I would look at would be to assign all polygons an ID based on the census block - then do your calculations? Looks like your for loop is doing a calculation for each feature which may take longer. – GISKid Jul 5 '16 at 13:19
  • The area is about 900 km2, however summing the polygons only takes a few seconds. My impression is that calling layer.updateFeature() for each polygon takes a lot of time. Do you know how to do some kind of a multi update? – Alexander Gogl Jul 5 '16 at 14:07
  • The loop takes about 7 minutes for each 5000 polygons to calculate and update. However after a while the progress slows down and takes hours for each iteration. – Alexander Gogl Jul 5 '16 at 14:40
1

Right, there is a way to build a dictionary of changed attributes and pass it directly to the provider's changeAttributeValues method after the whole for loop. The script would look like this:

# feats is a python list:
feats = [feat for feat in layer.getFeatures(QgsFeatureRequest().setFlags(QgsFeatureRequest.NoGeometry).setSubsetOfAttributes(limit_cols))
]

# sort feats
feats.sort(key=itemgetter(col_census_id))

...

# Get unique census ids
census_ids = []
areas = dict()
for feat in feats:
    census_id = feat[col_census_id]
    if census_id not in census_ids:
        census_ids.append(census_id)
        areas[census_id] = 0

# Sum building area per census block
for i in range(len(feats)):
    # Sum area
    census_id = feats[i][col_census_id]
    areas[census_id] += feats[i][col_bldg_area]

# Create the changed attributes dict
attrsMap = {}
col_share_idx = layer.fieldNameIndex( col_share )
col_bldg_pop_idx = layer.fieldNameIndex( col_bldg_pop )
i, j, k = 0, 0, 0
for feat in feats:
    # Share of area
    share = feat[col_bldg_area] / float(areas[feat[col_census_id]])    
    attrsMap[feat.id()] = { 
        col_share_idx: share, 
        col_bldg_pop_idx: int(feat[col_pop] * share) }

    ...

# Save changes
layer.dataProvider().changeAttributeValues( attrsMap )
return layer

It would be interesting to know how long the script takes with this provider-based approach.

  • 1
    Hi Carrillo! Thank you for pointing out the hint. I will try it during the christmas holidays and report how it works out. – Alexander Gogl Dec 30 '16 at 18: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.