[GIS] Using each feature of vector layer as clip mask over separate point layer, storing results in either a Python list or in-memory scratch layer

bufferclipqgisqgis-processingvector-layer

I have a point layer (representing cities), and a set of vector layers whose features are hundreds of concentric discs, each of which I'd like to use as a clipping layer to select only the cities that lie within that disc, perform some operations on those cities, and then move to the next disc.

Here are some screenshots of the project. The concentric blue-grey discs are the clipping polygons – 25/50/100/150 mile buffers around military base locations(represented as green diamonds). The red dots are the city features (displayed in graduated sizes by violent crime rate) that I'm trying to use as the input layer :

Screenshot of QGIS mapping project, concentric buffers
Zoomed in

I was trying to use the QGIS processing module to do this. However all of the documentation I've been able to find on the qgis:clip function shows its output going to a shapefile. For instance the answers to the question Is there a way to call the clip function in pyQGIS from the python console? suggest:

processing.runalg("qgis:clip",inputlayer,overlaylayer,"output_file.shp")
processing.runandload("qgis:clip",inputlayer,overlaylayer,"output_file.shp")

What I want to do instead is iterate through all of the discs in the clipping layer(s), and save the results from using each into a temporary in-memory vector layer instead of a shapefile.

Here's the (non-working) code that I tried to write to do this:

for disc in troopBuffers:
    # Clipping requires the cities layer and troop buffer disc to both be selected
    buffLayer.select(disc.id())
    cityLayer.selectAll()

    clipLayer = QgsVectorLayer("Point", "temporary_points", "memory:tmpClipLayer")
    QgsMapLayerRegistry.instance().addMapLayer(clipLayer)
    processing.runandload("qgis:clip", cityLayer, disc, clipLayer)

    clippedCities = [f for f in clipLayer.getFeatures()]

    for c in clippedCities:
        # do something with each city that lies within this disc ...

    buffLayer.removeSelection()
    QgsMapLayerRegistry.instance().removeMapLayer(clipLayer.id())

… however, when I run this clippedCities appears to always be empty, so something is obviously wrong.

How would I clip like this without actually creating shapefiles on disk for each iteration?

Can I just have processing output the results of each clip to an in-memory "scratch layer", or (even better) just get a list of QgsFeature objects from the cities layer for each disc?

Or am I going to have to use something else besides processing here?

Best Answer

When I get your question right, then you want treat the cities falling inside some polygon as a group and manipulate them in some way, without producing a number of files as a result of some process.

My suggestions:

When each city should be manipulated with respect to the polygon it is within, then use a code like this:

buff_layer = QgsMapLayerRegistry.instance().mapLayersByName('buffLayer')[0]
city_layer = QgsMapLayerRegistry.instance().mapLayersByName('cityLayer')[0]

# switch to edit mode
city_layer.startEditing()

# eg get index of a field to manipulate 
fni = city.fieldNameIndex('name')

# loop over polygons
for buff in buff_layer.getFeatures():
    for city in city_layer.getFeatures():
        if buff.geometry().contains(city.geometry()):
            # do something with a single city
            # access to buff[attribute]
            city_layer.changeAttributeValue(city.id(), fni, buff['name'])
city_layer.commitChanges()

In the example above each city receives in its field name the name of the polygon it is within.

When you need access to all cities within each polygon, you may use a code like this:

buff_layer = QgsMapLayerRegistry.instance().mapLayersByName('buffLayer')[0]
city_layer = QgsMapLayerRegistry.instance().mapLayersByName('cityLayer')[0]

# loop over polygons
for buff in buff_layer.getFeatures():
    selection = []
    for city in city_layer.getFeatures():
        if buff.geometry().contains(city.geometry()):
            selection.append(city.id())

    # select all cities in a single buffer    
    city_layer.setSelectedFeatures(selection)

    # do something with the selected cities

In this latter case please be aware that you cannot call functions outside this script repeatedly, since it is running in the same thread as Qgis itself. So the application is waiting for the script to finish before something else is started. Further readings: Select feature one by one not working in pyQGIS? and How do I prevent Qgis from being detected as "not responding" when running a heavy plugin?

Related Question