1

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

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

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.