I have two layers with 20+ polygon features, layer1 and layer 2. I would like to extract attributes from all of the features in layer 2 that are within 1km of layer 1, for each feature in layer 1. I am doing this by iterating through the features in layer 1, and buffering each feature by 1km. Within this loop, I am iterating through the features in layer 2, and checking whether these features intersect with the 1km buffer:
siteboundary=None
for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
if lyr.name() == "@@@@@":
layer 1 = lyr
break
targetlayer=None
for lyr2 in QgsMapLayerRegistry.instance().mapLayers().values():
if lyr2.name() == "@@@@2":
targetlayer = lyr2
break
siteboundaryshapes = siteboundary.getFeatures()
targetlayershapes = targetlayer.getFeatures()
for count,i in enumerate(siteboundaryshapes):
siteboundarygeom = i.geometry()
siteboundarybuffer = siteboundarygeom.buffer(1000,25)
siteboundaryatt = i.attributes()
print (siteboundaryatt[1])
for j in targetlayershapes:
targetlayergeom = j.geometry()
targetlayeratt = j.attributes()
if targetlayergeom.intersects(siteboundarybuffer):
print (targetlayeratt[1])
else:
continue
I am trying to return the attributes of the features in layer 2 that intersect with the buffer on layer 1 in the following format:
layer 1 feature 1 attribute
layer 1 feature 2 attribute
layer 2 feature 1 attrbute
layer 2 feature 2 attribute
layer 1 feature 3 attribute
layer 2 feature 1 attribute
layer 1 feature 4 attribute
layer 2 feature 1 attribute
etc etc
layer 1 feature 20 attribute
The above currently does not work; it loops through all of the features in layer 1, but only returns attributes for layer 2 in the first instance that layer 2 intersects with the buffer in layer 1. I have tried moving the indent of the if statement, but this just loops through the features in layer 2 for the first feature in layer 1 and then stops looping through layer 2 thereafter.
1 Answer 1
I would make sure you were using an index to cut down on needing to go through the whole 2nd layer each time.
Something like this should do it:
from qgis.core import *
import processing
layer1 = processing.getObject('ne_10m_populated_places')
layer2 = processing.getObject('ne_10m_admin_0_countries')
all_features = {}
index = QgsSpatialIndex() # Spatial index
for ft in layer1.getFeatures():
index.insertFeature(ft)
all_features[ft.id()] = ft
selection = [] # This list stores the features which contains at least one point
features = layer2.getFeatures()
for feat in features:
inGeom = feat.geometry().buffer(1000,25)
idsList = index.intersects(inGeom.boundingBox())
#print "found %d points"%len(idsList)
for id in idsList:
testGeom = all_features[id].geometry()
testGeomP = testGeom.asPoint()
#print testGeomP.x(),testGeomP.y()
if inGeom.intersects(testGeomP):
if feat not in selection:
selection.append( all_features[id])
else:
pass
# Select all the polygon features which contains at least one point
layer1.setSelectedFeatures([k.id() for k in selection])