Iterating over point layer and select another layer's points in radius using PyQGIS?

Question

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

Answer 1

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?

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