Performing a Spatial Query in a loop in PyQGIS

Question

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 1

You don't need a special function (as "Ray Casting"), everything is in PyQGIS (contains() in PyQGIS Geometry Handling)

polygons = [feature for feature in polygons.getFeatures()]
points = [feature for feature in points.getFeatures()]
for pt in points: 
     point = pt.geometry() # only and not pt.geometry().asPolygon() 
     for pol in polygons:
        poly = pol.geometry()
        if poly.contains(point):
             print "ok" 

or in one line

 polygons = [feature for feature in polygons.getFeatures()]
 points = [feature for feature in points.getFeatures()]
 resulting = [pt for pt in points for poly in polygons if poly.geometry().contains(pt.geometry())]
 print len(resulting)
 ...

You can also directly use

[pt.geometry().asPoint() for pt in points for poly in polygons if poly.geometry().contains(pt.geometry())]

The problem here is that you must iterate through all the geometries (polygons and points). It is more interesting to use a bounding spatial index : you iterate only through the geometries which have a chance to intersect with your current geometry ('filter', look at How to efficiently access the features returned by QgsSpatialIndex?)

Answer 2

You can use the "Ray Casting" algorithm that I have slightly adapted for using with PyQGIS:

def point_in_poly(point,poly):
    x = point.x()
    y = point.y()

    n = len(poly)
    inside = False

    p1x,p1y = poly[0]
    for i in range(n+1):
        p2x,p2y = poly[i % n]
        if y > min(p1y,p2y):
            if y <= max(p1y,p2y):
                if x <= max(p1x,p2x):
                    if p1y != p2y:
                        xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                    if p1x == p2x or x <= xints:
                        inside = not inside
        p1x,p1y = p2x,p2y

    return inside

## Test
mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

#For polygon 
polygon = [feature.geometry().asPolygon() 
            for feature in layers[1].getFeatures()]

points = [feat.geometry().asPoint() 
           for feat in layers[0].getFeatures()]

## Call the function with the points and the polygon
count = [0]*(layers[1].featureCount())

for point in points:
    i = 0
    for feat in polygon:
        if point_in_poly(point, feat[0]) == True:
            count[i] += 1
        i += 1

print count

Applied to this situation:

the result, at the Python Console, was:

[2, 2]

It worked.

Editing Note:

Code with the more concise proposal of gene:

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

count = [0]*(layers[1].featureCount())

polygon = [feature
           for feature in layers[1].getFeatures()]

points = [feature
          for feature in layers[0].getFeatures()]

for point in points:

    i = 0

    geo_point = point.geometry()

    for pol in polygon:
        geo_pol = pol.geometry()

        if geo_pol.contains(geo_point):
            count[i] += 1
        i += 1

print count

Answer 3

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