1

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 don't think there's a ready-made tool. But depending on your use-case a simple pyqgis srcript can help. Can you tell more about what is the goal of the analysis? How would you like to aggregate cells? (merge neighbor cells till the sum is > X? in what order? etc.) – spatialthoughts Mar 9 '15 at 10:02
  • I don't have specific demands on how the aggregation is done, but as you said, I'd like to merge cells until sum > X. The goal is to preserve personal integrity by not showing cells with too few data points (ie people) in them. It would be nice if it could done in such a way that, given the same data set, you get the same result. But anything that gets me joined cells with sufficient population is fine. – hexamon Mar 9 '15 at 10:08
  • Related to / possible duplicate of: gis.stackexchange.com/questions/31324 (Also note other linked questions from there). – Chris W Mar 9 '15 at 18:19
4

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.

before

after

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)
  • Thanks a lot for this! It answers my question, but I realised reality (as usual) is a bit more complicated than I considered. Your script works for fairly evenly distributed grids, but my data has a lot of 0 value cells in the countryside, and as such the recursive iterations on isolated cells take a very long time, leading to memory errors. But it's a great start! Thanks again! (Ps. The script uses the first column in the source file, not the column named 'population') – hexamon Mar 12 '15 at 14:59
  • Cool. Depending on your use case, you can either stop the iteration after N tries or skip aggregation for a cell if the population is 0. – spatialthoughts Mar 13 '15 at 3:30

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.