Utilizar cada característica de la capa vectorial como máscara de recorte sobre una capa de puntos separada, almacenando los resultados en una lista de Python o en una capa de raspado en memoria

Question

Tengo una capa de puntos (que representa las ciudades), y un conjunto de capas vectoriales cuyas características son cientos de discos concéntricos, cada uno de los cuales me gustaría utilizar como capa de recorte para seleccionar sólo las ciudades que se encuentran dentro de ese disco, realizar algunas operaciones en esas ciudades, y luego pasar al siguiente disco.

Aquí hay algunas capturas de pantalla del proyecto. Los discos concéntricos de color azul-gris son los polígonos de recorte: topes de 25/50/100/150 millas alrededor de las ubicaciones de las bases militares (representadas como diamantes verdes). Los puntos rojos son las características de la ciudad (mostradas en tamaños graduados por la tasa de criminalidad violenta) que estoy tratando de usar como capa de entrada:

Estaba tratando de usar el QGIS processing módulo para hacerlo. Sin embargo, toda la documentación que he podido encontrar sobre el qgis:clip muestra su salida hacia un shapefile. Por ejemplo las respuestas a la pregunta ¿Existe una forma de llamar a la función clip en pyQGIS desde la consola de python? sugieren:

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

Lo que quiero hacer en su lugar es iterar a través de todos los discos en la(s) capa(s) de recorte, y guardar los resultados de usar cada uno en un en memoria capa vectorial en lugar de un shapefile.

Aquí está el código (que no funciona) que intenté escribir para hacer esto:

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())

... sin embargo, cuando ejecuto esto clippedCities parece estar siempre vacío, por lo que es evidente que algo va mal.

¿Cómo podría recortar así sin crear realmente archivos shape en el disco para cada iteración?

¿Puedo tener processing enviar los resultados de cada clip a una "capa de raspado" en memoria, o (incluso mejor) simplemente obtener una lista de QgsFeature objetos de la capa de ciudades para cada disco?

¿O voy a tener que usar otra cosa además de processing ¿Aquí?

Answer 1

Si entiendo bien tu pregunta, entonces quieres tratar las ciudades que caen dentro de algún polígono como un grupo y manipularlas de alguna manera, sin producir un número de archivos como resultado de algún proceso.

Mis sugerencias:

Cuando cada ciudad deba ser manipulada con respecto al polígono en el que se encuentra, entonces utilice un código como este:

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()

En el ejemplo anterior cada ciudad recibe en su campo name el name del polígono en el que se encuentra.

Cuando necesite acceder a todas las ciudades dentro de cada polígono, puede utilizar un código como éste:

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

En este último caso ten en cuenta que no puedes llamar a funciones fuera de este script repetidamente, ya que se ejecuta en el mismo hilo que el propio Qgis. Así que la aplicación está esperando a que el script termine antes de iniciar otra cosa. Otras lecturas: ¿Seleccionar una característica por una no funciona en pyQGIS? y ¿Cómo puedo evitar que Qgis sea detectado como "no responde" cuando se ejecuta un plugin pesado?