3

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

4
  • 2
    Welcome to GIS:SE @JoshuaKidd! Remove self, you don't really need this in a script. Also take a look at Using a QGIS spatial index to speed up your code which may be helpful :)
    – Joseph
    Jan 18, 2017 at 11:42
  • Thank you! If what I've written works, what should I replace self with? I can see in the example you posted that he finds if features in a single vector layer touch another. I'm guessing if I try and adapt what he's written I'll have to find a way to see if a point from one layer overlaps a vector in another. Jan 18, 2017 at 12:01
  • If you have acces to GeoPandas Python module, it could be very straightforward to do that with 'sjoin' method and 'intersects' option (see my answer).
    – xunilk
    Jan 18, 2017 at 12:10
  • @JoshuaKidd -No need to replace self, just remove it. In terms of the example link I posted, you could replace touches with intersects.
    – Joseph
    Jan 18, 2017 at 12:20

2 Answers 2

6

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:

enter image description here

and running the above code, I obtain this output:

enter image description here

13
  • I'm getting an error saying processing not defined. Jan 18, 2017 at 13:26
  • That's weird.. Try typing import processing at the top of the code. However, I need to edit it because at this moment it selects all the polygons for which the bounding box contains a point (it's a minor fix, I'll edit the answer as soon as possible, sorry)
    – mgri
    Jan 18, 2017 at 13:31
  • I just added 'import processing' and it's running now :) Jan 18, 2017 at 13:33
  • WOW, that ran in 40 seconds as opposed to my previous best 16 min! Jan 18, 2017 at 13:37
  • Well! I edited the code, now it should work as expected! Please, let me know if it works for you! =)
    – mgri
    Jan 18, 2017 at 13:37
1

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:

enter image description here

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:

enter image description here

5
  • Alright, I'll give this a try. I'm just installing geopandas now. My main problem isn't that I can't get it to select the features, it's that I'm looking for the quickest way to do it. Does this method use spatial indexing, or something similar? Jan 18, 2017 at 12:23
  • Yes, it does. So, you also need Rtree python module (in my case installed with easy_install) and libspatialindex library (from my Debian Linux repository). If you are Windows user see: gis.stackexchange.com/questions/120527/…
    – xunilk
    Jan 18, 2017 at 12:27
  • hmmm, this geopandas install is more involved than I thought. I'm using spyder through anaconda Jan 18, 2017 at 12:39
  • I use GNU/Debian Linux where installation of these libraries are much easier.
    – xunilk
    Jan 18, 2017 at 12:45
  • Ok, i got it working, and it ran for 16 minutes which is good. A problem though, is that it creates entries for each point in a poly. Jan 18, 2017 at 13:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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