3

I found this code snippet in Selecting features within certain buffer distance of selected feature using PyQGIS?, which is almost what I'm looking for, except I don't have polylines but points and I want to find another layers points in radius.

I tried to edit the code above, but I'm confused when to switch layers and how to select the points in radius properly, so that I can call my "do something"-function.

Here my code so far:

def selectPointsInRadius(bufDist):
  inputlayer = QgsMapLayerRegistry.instance().mapLayersByName('points')[0]
  print ('points-layer found')
  #stores features in a dictionary for fast access
  featDict = {feature.id(): feature for (feature) in inputlayer.getFeatures()}
  #spatial index for faster spatial queries
  featIdx = QgsSpatialIndex(inputlayer.getFeatures())
  try:
    itlayer = QgsMapLayerRegistry.instance().mapLayersByName('itpoint')[0]
    print ('itpoint-layer found')
    itlayer.selectAll()
    for point in itlayer.selectedFeatures():
        xyList = point.geometry().asPoint()
        print(xyList)
        minx = min([xy[0] for xy in xyList])   #minimum x-coordinate
        maxx = max([xy[0] for xy in xyList])   #maximum x-coordinate
        miny = min([xy[1] for xy in xyList])   #minimum y-coordinate
        maxy = max([xy[1] for xy in xyList])   #maximum y-coordinate
        ids = featIdx.intersects(QgsRectangle(minx - bufDist, miny - bufDist, maxx + bufDist, maxy + bufDist))
        #lyr.selectByIds(ids)  #use this if you don't need total accuracy, much faster!!
        tids = []  #true ids of nearby features that we want to select
        for id in ids:
            if featDict[id].geometry().distance(line.geometry()) < bufDist:
                print ("found features in radius")
                tids.append(id)
        print ("selected features in radius")
        lyr.selectByIds(tids, 1)   #add to selection,
        #call doSomething-function here?
  except IndexError:
    print ('itpoint does not exist')
  • 1
    You could also use the Processing algorithm Select by location: processing.runalg("qgis:selectbylocation", inputlayer, itlayer, u'within', bufDist, 0) – Joseph Aug 9 '17 at 10:06
  • 1
    @Joseph, why put this answer as a comment? I would normally search for the answers among the answers and not in the comments - unless the answers doesn't work. I haven't tried your suggestion, but it looks like a good and valid answer :) – Chau Aug 9 '17 at 10:39
  • @Chau - Thanks! But since the question includes code from the API, I'm assuming the OP is looking to see where the problem lies instead of looking to replace it entirely :) – Joseph Aug 9 '17 at 10:48
  • I've seen the processing algorithm several times, but don't know how to use it 😅😬 do I just have to copy and paste it in my editor and execute? – Nico Aug 9 '17 at 11:23
  • @Nico - First you will need to use import processing, then define your layers as you have done in your script (i.e. inputlayer = QgsMapLayerRegistry.instance().mapLayersByName('points')[0]). Then just run the command processing.runalg("qgis:selectbylocation", inputlayer, itlayer, u'within', 50, 0) to create a bounding box precision of 50 units. It might not be what you're looking for as it doesn't create a search buffer, only a bounding box. – Joseph Aug 10 '17 at 11:45
2

I perform things like this with this code snippets:

inputlayer = QgsMapLayerRegistry.instance().mapLayersByName('pieces_sepc')[0]
itlayer = QgsMapLayerRegistry.instance().mapLayersByName('POINT')[0]
dist = 60

def selectPointsInRadius(inputlayer, itlayer, dist):
    #iterate over feature of inputlayer and get geometry
    for feat in inputlayer.getFeatures():
        itlayer.removeSelection()
        geom = feat.geometry()
        #iterate over feature of itlayer and get geometry
        for itfeat in itlayer.getFeatures():
            itgeom = itfeat.geometry()
            #test if geom of itlayer is in the buffer of inputlayer, if yes select the feature 
            if geom.buffer(dist,5).intersects(itgeom):
                itlayer.select(itfeat.id())
        for sfeat in itlayer.selectedFeatures():
            #do something with the selected feature
            pass

selectPointsInRadius(inputlayer, itlayer, dist)

At the end every features from itlayer which are at a distance less than 75 meters from a feature of inputlayer is affect by the things you put where I comment #do something with the selected feature

The point is that I have never test it on huge datasets, so I have no idea about the performance of my code. In the small dataset I work with it's efficient.

  • thx, also looks great!... but unfortunately can't test it, because I have to fix my layers points (which don't have a geometry defined)... but I'll get back to you soon! thx again! ...@huge datasets: I'm a beginner and my other functions aren't memory optimized as well - therefore I'm always keeping an eye on the memory usage ;-) – Nico Aug 9 '17 at 10:46
  • You just helped me to fix AND improve the other code there (gis.stackexchange.com/questions/251483/…), it should actually be my "inputlayer" in this code here... And yeah, I don't get an error message anymore, but that's all - seems that it doesn't find any matches (that's unfortunately impossible). what does the "5" stand for in "if geom.buffer(dist,5).intersects(itgeom):"? (just trying to understand what's going wrong ;-)) – Nico Aug 9 '17 at 14:57
  • Actually it probably find matches. but at each iteration it remove the selection to work only with features which match with the current feat from inputlayer loop. If your last feat has no matches, at the end of the function, it will seems that nothing happen... add some functions in the last loop and they will be execute on the selected itfeat (the itfeat that are within bufDist from the cureent feat). Please give more information about what you want to do and I would be able to help. – YoLecomte Aug 9 '17 at 15:11
  • the 5 is the number of segments generated to create a pseudo-circle buffer – YoLecomte Aug 9 '17 at 15:12
  • I put some prints in and it seems that it runs all for loops but never gets into the if clause (a print there is never printed); although I put the buffer to 100 meter and there are points everywhere nearby...? What I want to do further on is kind of complex and may be more than one other code snippet away... I want to check what attributes "amenity" the selected points have (there are quite a lot, approximately 10 like "restaurant" or "school" or "museum" etc.) and count each of them. And depending on which and how many are counted, call another function... (is it possible to chat? :-) – Nico Aug 9 '17 at 15:38
3

To select all points in the points layer which is within bufDist of the points in the itpoint layer, you can use the following:

def selectPointsInRadius(bufDist):
    inputlayer = QgsMapLayerRegistry.instance().mapLayersByName('points')[0]
    print ('points-layer found')
    #stores features in a dictionary for fast access
    featDict = {feature.id(): feature for (feature) in inputlayer.getFeatures()}
    #spatial index for faster spatial queries
    featIdx = QgsSpatialIndex(inputlayer.getFeatures())
    try:
        itlayer = QgsMapLayerRegistry.instance().mapLayersByName('itpoint')[0]
        print ('itpoint-layer found')
        tids = []  #true ids of nearby features that we want to select    
        for point in itlayer.getFeatures():
            p = point.geometry().asPoint()
            ids = featIdx.intersects(QgsRectangle(p[0] - bufDist, p[1] - bufDist, p[0] + bufDist, p[1] + bufDist))        
            for id in ids:
                if featDict[id].geometry().distance(point.geometry()) < bufDist:
                    print ("found features in radius")
                    # Add the id of the point in the itpoint layer
                    tids.append(id)
        inputlayer.selectByIds(tids, 1)   #add to selection,    
        #call doSomething-function here?
    except IndexError:
        print ('itpoint does not exist')

Is that what you are looking for?

  • I think so, yes! :D but somehow I get an error message now "p = point.geometry().asPoint() AttributeError: 'QgsGeometry' object has no attribute 'asPoint'" do you know how to fix it? – Nico Aug 9 '17 at 10:03
  • sorry pasted the wrong message: “p = point.geometry().asPoint() AttributeError: 'NoneType' object has no attribute 'asPoint'“ ...does it mean that my itpoint-layers points don't have a point geometry yet? (I created them with this function before: gis.stackexchange.com/questions/251225/… ) – Nico Aug 9 '17 at 10:07
  • @Nico - Yes, it means atleast one of your point features does not have a geometry defined. – Joseph Aug 9 '17 at 10:14
  • @Nico, just be aware that with this code, you may miss some features that are within your bufDist because you do the intersection with a rectangle and not with a circle... – YoLecomte Aug 9 '17 at 10:39
  • @YoLecomte, You propose a method which buffers the point by radius of the circle you want to search within. I would think that the buffer method is more expensive than a normal distance calculation between two points. Further more the rectangle (or square in this case) provides a larger search area which in every case contains the circle you are using - thus ensuring that he'll get at minimum the same results as you. – Chau Aug 9 '17 at 10:45

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.