Selección de polígonos que contienen al menos un punto con índice espacial utilizando PyQGIS

Question

Tengo un shapefile de puntos y otro de polígonos. Me gustaría seleccionar todos los polígonos que tienen al menos un punto en ellos.

El problema con el que me encuentro es el tiempo que lleva esto. Tengo 1 millón de puntos y unos 320.000 polígonos, por lo que utilizar la consulta espacial lleva demasiado tiempo. He oído que tendría que escribir un script en Python con indexación espacial para obtener un resultado factiblemente rápido, pero no tengo ni idea de cómo abordar esto.

No me importa que sea un poco inexacto si corre mucho más rápido.

Lo que he tratado de improvisar a partir de otras preguntas de stack overflow es:

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

Esto arroja un NameError:

el nombre 'self' no está definido

Answer 1

He tomado un rollo de acuerdo a las necesidades de @Joshua Kidd que buscaba la solución más rápida (y no la más correcta).

El código publicado a continuación seleccionará una característica poligonal cuando su caja delimitadora contenga al menos un punto. En cambio, la solución más correcta (es decir, la selección de una característica poligonal cuando su geometría real contiene al menos un punto) se informa en este revisión previa de mi respuesta.

La solución más rápida que podría aplicar en PyQGIS sería esta:

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

A partir de esta muestra de capas:

y ejecutando el código anterior, obtengo esta salida:

Answer 2

Si tiene acceso al módulo GeoPandas de Python, podría ser muy sencillo hacerlo con sjoin método y "se cruza opción. Probé mi enfoque con la siguiente situación:

El código era:

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

Después de ejecutar el código en la consola de Python de QGIS, cargué el shapefile resultante (polyWithPoints.shp) y era el esperado: todos los polígonos (color azul) que tienen al menos un punto en ellos: