3

I'm very new to using QGIS, what I have is 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

  • 1
    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 '17 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. – Joshua Kidd Jan 18 '17 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 '17 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 '17 at 12:20
4

EDIT 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

  • I'm getting an error saying processing not defined. – Joshua Kidd Jan 18 '17 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 '17 at 13:31
  • I just added 'import processing' and it's running now :) – Joshua Kidd Jan 18 '17 at 13:33
  • WOW, that ran in 40 seconds as opposed to my previous best 16 min! – Joshua Kidd Jan 18 '17 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 '17 at 13:37
1

If you have acces to GeoPandas Python module, it could be very straightforward to do that with 'sjoin' method and 'intersects' option. I tried out my approach with 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 resulting shapefile (polyWithPoints.shp) and it was as expected: all the polygons (blue color) which have at least one point in them:

enter image description here

  • 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? – Joshua Kidd Jan 18 '17 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 '17 at 12:27
  • hmmm, this geopandas install is more involved than I thought. I'm using spyder through anaconda – Joshua Kidd Jan 18 '17 at 12:39
  • I use GNU/Debian Linux where installation of these libraries are much easier. – xunilk Jan 18 '17 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. – Joshua Kidd Jan 18 '17 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.