Buscar Polígono en "Capa A" que contenga Polígono en "Capa B"

Question

Todavía no he encontrado a nadie con el mismo problema, así que publico una nueva pregunta. El tema parece bastante estándar, así que estoy seguro de que hay una manera de hacerlo en QGIS.

Tengo dos capas:

Capa de regiones Capa de códigos postales

ambos con características multipoligonales.

Me gustaría recorrer todos los códigos postales, encontrar cuál es la región que contiene ese código postal y, finalmente, escribir el ID de la región dentro del código postal.

Hasta ahora he intentado escribir algo en la consola de Python, que puedo compartir aquí:

from qgis.utils import iface
from PyQt4.QtCore import QVariant
_ID_FIELD = 'RS'
_NAME_FIELD = 'GEN'
_PLZ_FIELD = 'plz'
_NEW_ID_FIELD = 'PARENT_ID'
_NEW_NAME_FIELD = 'PARENT_NAME'
Land_layer = iface.mapCanvas().layer(1)
PLZ_layer = iface.mapCanvas().layer(0)
# Create 2 new fields in the layer that will hold the list of neighbors and sum
# of the chosen field.

#PLZ_layer.startEditing()
#PLZ_layer.dataProvider().addAttributes(
#        [QgsField(_NEW_ID_FIELD, QVariant.String),
#         QgsField(_NEW_NAME_FIELD, QVariant.Int)])
#PLZ_layer.updateFields()

## Create a dictionary of all features
PLZ_feature_dict = {f.id(): f for f in PLZ_layer.getFeatures()}
Land_feature_dict = {f.id(): f for f in Land_layer.getFeatures()}
# Build a spatial index
index = QgsSpatialIndex()
print 'PLZ Layer: %s' % PLZ_layer.name()
print 'Landkreis Layer: %s' % Land_layer.name()

for f in Land_feature_dict.values():
    index.insertFeature(f)
    print 'Inserting %s in Landkreis index' % f[_NAME_FIELD]

for f in PLZ_feature_dict.values():
    print 'Working on %s' % f[_PLZ_FIELD]
    geom = f.geometry()
    intersect_ids = index.intersects(geom.boundingBox())
    for id in intersect_ids:
        tempf = Land_feature_dict[id]
        print 'got parent Landkreis %s' % tempf[_NAME_FIELD]

El problema es que obtengo el siguiente error en la línea

intersect_ids = index.intersects(geom.boundingBox())

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "//VAIMUCO9/G114305$/Anwendungsdaten/Desktop/PLZ/ContainerFinder.py", line 34, in  <module>
    intersect_ids = index.intersects(geom.boundingBox())
AttributeError: 'NoneType' object has no attribute 'boundingBox'

Además, no estoy muy seguro de que deba utilizar estas funciones, e incluso este enfoque. ¿Puede alguien sugerirme alguna forma más inteligente?

Answer 1

Unas cuantas cosas. En primer lugar, .getFeatures() devuelve un iterador, y no las características en sí, por lo que ponerlas en un diccionario sólo sirve para consumir memoria.

Rehice tu trabajo usando un par de capas que tengo para el área de Phoenix, AZ y funcionó bien. Véase más abajo.

El código original tampoco garantiza que la característica esté contenida, sino que se intersecte. He incluido una prueba para asegurarse de que se cruza.

No me preocupé por los nombres, ya que el mío sería diferente al tuyo, pero el código funciona y puedes añadirle los mensajes.

PLZ_layer = iface.mapCanvas().layer(0)
Land_layer = iface.mapCanvas().layer(1)

# Build a spatial index
index = QgsSpatialIndex()
for f in Land_layer.getFeatures():
    a=index.insertFeature(f)

print 'PLZ Layer: '+ PLZ_layer.name()
print 'Landkreis Layer: '+ Land_layer.name()

for f in PLZ_layer.getFeatures():
    geom = f.geometry()
    intersect_ids = index.intersects(geom.boundingBox())

    #Now we filter only the IDS that intersected our f feature
    request = QgsFeatureRequest()
    a=request.setFilterFids(intersect_ids)

    #And we will loop through them in order to see which actually contain our f feature
    for feat in Land_layer.getFeatures(request):
        if feat.geometry().contains(geom):
            print 'FOUND A FATHER'