9

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:

enter image description here

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.

3 Answers 3

11

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

2
  • 1
    Also see nathanw.net/2013/01/04/…
    – Nathan W
    Commented Nov 3, 2015 at 0:07
  • @gene, can i ask a question, is there anyway to avoid using two for loops, in large datasets this is not good for complexity reasons? Commented Oct 4, 2021 at 15:31
5

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:

enter image description here

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
2
  • Great reference and great answer! I will mark the one I just posted, however, as the solution, as it is a little easier to implement. You should be rewarded with a lot of upvotes though. +1 from me, for sure. Commented Nov 2, 2015 at 13:02
  • You don't need to specify if geo_pol.contains(geo_point) == True: because it is implicit in if geo_pol.contains(geo_point) (always True)
    – gene
    Commented Nov 2, 2015 at 17:55
4

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

enter image description here

1
  • I tried your solution. It works only if you have one polygons, otherwise only the points within the first polygon will be selected
    – ilFonta
    Commented Apr 14, 2019 at 11:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.