Search code examples
polygonpointqgisspatial-index

QGIS Select polygons which intersect points with python


I'm very new to using QGIS, what I have is a points shapefile and a polygon shapefile. I would like to select all the polygons which have at least one point in them. The problem I'm running into is how long this takes. I have 1 million points and about 320,000 polygons, so using spatial query takes far too long. I've heard that I'd need to write a python script with spatial indexing to get a feasibly quick result, but I have no idea how to approach this. Any help would be greatly appreciated.

What I've tried to cobble together from other stack overflow questions is:

pointProvider = self.pointLayer.dataProvider()
all_point = pointProvider.getFeatures()
delta = 0.1

for point in all_point:

    searchRectangle = QgsRectangle(point.x() - delta, point.y()  - delta, point.x() + delta, point.y() + delta)

    candidateIDs = line_index.intesects(searchRectangle)

    for candidateID in candidateIDs:
        candFeature == rotateProvider.getFeatures(QgsFeatureRequest(candidateID)).next()
        if candFeature.geometry().contains(point):

            break

This throws up a NameError: name 'self' is not defined


Solution

  • I found an answer over on GIS Stack Exchange, which you can find here

    The code I used was:

    from qgis.core import *
    import processing
    
    layer1 = processing.getObject('MyPointsLayer')
    layer2 = processing.getObject('MyPolygonsLayer')
    
    index = QgsSpatialIndex() # Spatial index
    for ft in layer1.getFeatures():
        index.insertFeature(ft)
    
    selection = [] # This list stores the features which contains at least one point
    for feat in layer2.getFeatures():
        inGeom = feat.geometry()
        idsList = index.intersects(inGeom.boundingBox())
        if idsList:
            selection.append(feat)
    
    # Select all the polygon features which contains at least one point
    layer2.setSelectedFeatures([k.id() for k in selection])