python - Performing a Spatial Query in a loop in PyQGIS

What I am trying to do: loop through a point shapefile and select each point that falls into a polygon.

The following code is inspired by a spatial query example I found in a book:

mitte_path = r"D:\PythonTesting\SelectByLocation\mitte.shp"
punkte_path = r"D:\PythonTesting\SelectByLocation\punkte.shp"

polygon = QgsVectorLayer(mitte_path, 'Mitte', 'ogr')

points = QgsVectorLayer(punkte_path, 'Berlin Punkte', 'ogr')

QgsMapLayerRegistry.instance().addMapLayer(polygon)
QgsMapLayerRegistry.instance().addMapLayer(points)

polyFeatures = polygon.getFeatures()

pointsCount = 0

for poly_feat in polyFeatures:

    polyGeom = poly_feat.geometry()
    pointFeatures = points.getFeatures(QgsFeatureRequest().setFilterRect(polyGeom.boundingBox()))
    for point_feat in pointFeatures:
        points.select(point_feat.id())
        pointsCount += 1

print 'Total:',pointsCount

This works, and it does select datasets, but the problem is that it selects by bounding box, hence obviously returning points I am not interested in:

How could I go about only returning points within the polygon without using qgis:selectbylocation?

I have tried using the within() and intersects() methods, but as I was not getting them to work, I resorted to the code above. But perhaps they're the key after all.

Answer

With some advice from a workmate I finally got it to work using within().

General logic

1. get features of polygon(s)


2. get features of points


3. loop through each feature from polygon file, and for each:
   
   
   
   
   • get geometry
   
   
   • loop through all points
     
     
     
     • get geometry of single point
     
     
     • test if geometry is within geometry of polygon
     
     
     
     
   
   
   
   

Here is the code:

mitte_path = r"D:\PythonTesting\SelectByLocation\mitte.shp"
punkte_path = r"D:\PythonTesting\SelectByLocation\punkte.shp"

poly = QgsVectorLayer(mitte_path, 'Mitte', 'ogr')
points = QgsVectorLayer(punkte_path, 'Berlin Punkte', 'ogr')

QgsMapLayerRegistry.instance().addMapLayer(poly)
QgsMapLayerRegistry.instance().addMapLayer(points)


polyFeatures = poly.getFeatures()
pointFeatures = points.getFeatures()

pointCounter = 0

for polyfeat in polyFeatures:
    polyGeom = polyfeat.geometry()
    for pointFeat in pointFeatures:
        pointGeom = pointFeat.geometry()
        if pointGeom.within(polyGeom):

            pointCounter += 1
            points.select(pointFeat.id())

print 'Total',pointCounter

This would also work with intersects() instead of within(). When using points it does not matter which one you would use, as they will both return the same result. When checking for lines/polygons, however, it can make an important difference: within() returns objects that are entirely within, whereas intersects() retuns objects that are entirely within and that are partly within (i.e. that intersect with the feature, as the name indicates).