Selecting polygons which contain at least one point with spatial index using PyQGIS

Question

I have 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.

I don't mind if it's a little inaccurate if it runs significantly quicker.

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

Answer 1

I have taken a rollback according to the needs of @Joshua Kidd, which was looking for the fastest (and not the most correct) solution.

The code posted in the following will select a polygon feature when its bounding box contains at least one point. Instead, the most correct solution (i.e. the selection of a polygon feature when its real geometry contains at least one point) is reported at this previous revision of my answer.

The fastest solution I could apply in PyQGIS would be this one:

##Points=vector point
##Polygons=vector polygon

from qgis.core import *
import processing

layer1 = processing.getObject(Points)
layer2 = processing.getObject(Polygons)

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

Starting from this sample layers:

and running the above code, I obtain this output:

Answer 2

If you have access to the GeoPandas Python module, it could be very straightforward to do that with sjoin method and 'intersects' option. I tried out my approach with the next situation:

The code was:

from geopandas import gpd 

points = gpd.GeoDataFrame.from_file('/home/zeito/pyqgis_data/random_points.shp')  
polys = gpd.GeoDataFrame.from_file('/home/zeito/pyqgis_data/polygon8.shp') 
polyWithPoints = gpd.sjoin(polys, points, op='intersects') 

polyWithPoints.to_file('/home/zeito/pyqgis_data/polyWithPoints.shp')

After running the code at the Python Console of QGIS, I loaded the resulting shapefile (polyWithPoints.shp) and it was as expected: all the polygons (blue color) that have at least one point in them: