I'm looking for a tool or method to aggregate a vector fishnet based on a minimum value. Specifically, I have a fishnet containing a population value and I would like to join together cells that contain fewer people than X in an intelligent way.
Is there a straight forward way to do this in QGIS? Preferably without too much command line, but I'll get my hands dirty if that's what it takes.
I took a shot at it using python scripting. Below image shows an input grid with random population values. The script tries to merge the neighboring cells so that the value is > 50. Hope this is useful.
This is a pyqgis script, so you will need to open Plugins -> Python Console -> Show Editor and paste it there.. Then select your grid layer and click 'Run script'. This script assumes you have only 1 field named 'population' in your input grid. Change the threshold value to your needs. The script can be surely improved, but this gives you a general idea how to approach the problem.
from PyQt4.QtCore import QVariant
from qgis.utils import iface
THRESHOLD = 50
POPULATION_FIELD = 'population'
def aggregate(input_feature, threshold=THRESHOLD):
"""Recursive function to merge cells."""
if input_feature[POPULATION_FIELD] > threshold:
print 'Iteration complete. Returning feature %d with %d' % (
input_feature.id(), input_feature[POPULATION_FIELD])
if input_feature.id() not in feature_dict:
new_layer.startEditing()
new_layer.addFeatures([input_feature])
new_layer.commitChanges()
return
else:
print 'Population %d is less than threshold %d' % (
input_feature[POPULATION_FIELD], threshold)
geom = input_feature.geometry()
intersecting_ids = index.intersects(geom.boundingBox())
candidate_features = []
for intersecting_id in intersecting_ids:
# Look up the feature from the dictionary
intersecting_f = feature_dict[intersecting_id]
if (intersecting_f != input_feature and
not intersecting_f.geometry().disjoint(geom)):
candidate_features.append((
intersecting_id, intersecting_f[POPULATION_FIELD]))
# Sort by POPULATION FIELD value
sorted_candidate_features = sorted(candidate_features, key=lambda x: x[1])
if sorted_candidate_features:
# Choose the feature with the lowest population and merge it.
candidate_id, candidate_population = sorted_candidate_features[0]
print 'Found a candidate_feature %d with population %d' % (
candidate_id, candidate_population)
candidate_feature = feature_dict[candidate_id]
new_population = input_feature[POPULATION_FIELD] + candidate_population
new_feature = QgsFeature()
new_feature.setFields(iface.activeLayer().dataProvider().fields())
new_geom = geom.combine(candidate_feature.geometry())
new_feature.setGeometry(new_geom)
new_feature.setAttribute(POPULATION_FIELD, new_population)
new_layer.startEditing()
new_layer.deleteFeature(input_feature.id())
new_layer.deleteFeature(candidate_feature.id())
new_layer.commitChanges()
print ('Feature %d with population %d merged with '
'feature %d with population %d') % (
input_feature.id(), input_feature[POPULATION_FIELD],
candidate_id, candidate_population)
print 'Population is now at %d' % new_population
deleted_features.append(input_feature.id())
deleted_features.append(candidate_feature.id())
print 'Deleted features %d and %d' % (
input_feature.id(), candidate_feature.id())
index.deleteFeature(input_feature)
index.deleteFeature(candidate_feature)
# Make the recursive call with the new feature as input
aggregate(new_feature)
else:
print ('No more features to merge with. '
'Returning with %d') % input_feature[POPULATION_FIELD]
new_layer.startEditing()
new_layer.addFeatures([input_feature])
new_layer.commitChanges()
return
original_layer = iface.activeLayer()
# Copy to new layer
new_layer = QgsVectorLayer('Polygon?crs=EPSG:4326', 'combined', 'memory')
provider = new_layer.dataProvider()
provider.addAttributes([QgsField(POPULATION_FIELD, QVariant.Int)])
# First copy all features to our new layer
new_layer.startEditing()
for f in original_layer.getFeatures():
new_f = QgsFeature(f)
new_f[POPULATION_FIELD] = f[POPULATION_FIELD]
provider.addFeatures([new_f])
new_layer.updateExtents()
new_layer.commitChanges()
del original_layer
# Build a dictionary of features so we can access them easily
feature_dict = {f.id(): f for f in new_layer.getFeatures()}
# Build a spatial index so we can find neighoring features fast
index = QgsSpatialIndex()
for f in feature_dict.values():
index.insertFeature(f)
# Keep track of deleted features so they are not processed
deleted_features = []
for f in feature_dict.values():
if f.id() not in deleted_features:
aggregate(f)
new_layer.updateExtents()
new_layer.commitChanges()
QgsMapLayerRegistry.instance().addMapLayer(new_layer)
