I have a point-based vector layer, where each point represents a road accident. I am trying to write a function that:

  1. Loops through my accident layer
  2. creates a buffer, of a given radius, around each point/accident
  3. counts the number of accidents, from the accident layer, within each buffer
  4. if the number of accidents within (i.e. that intersect) the buffer is above a given threshold, then add this buffer as a polygon to a separate, pre-defined layer

What I have so far:

def get_clusters(layer, out_layer, buff_size, acc_threshold):

    segs = 100   #buffer segments

    source_layer = core.QgsMapLayerRegistry.instance().mapLayersByName(layer)[0]
    # get CRS
    CRS = source_layer.crs().authid()
    # make the output layer
    cluster_layer = QgsVectorLayer("Polygon?crs={0}".format(CRS),out_layer,"memory")

    # set up loop and empty list for areas of interest
    my_clusters = []
    for accident_feature in source_layer.getFeatures():
        acc_buffer = accident_feature.geometry().buffer(buff_size,segs)

        # At this point I'm not sure how to proceed. I think I need to use something like

        request = QgsFeatureRequest().setFilterRect(acc_buffer.geometry().boundingBox())
        areas_of_interest = source_layer.getFeatures(request)

        for area in areas_of_interest:
        ## pseudo-code - haven't got far enough to test this    
           if area.featureCount() >= acc_threshold:

but that's about as far as I can get. I need the intersection-object to be a circle (i.e. just the buffer polygon), but I can't find something like setFilterArea that I could pass the acc_buffer.geometry() in to.

I hope the above is clear, if someone could point me in the right direction I'm quite content working out the rest for myself but currently I'm a bit stuck!

1 Answer 1


It's not the most efficient way to do this. Instead of creating buffers and looking for intersections, you can do this with scipy.spatial library.

First get all of your points (accidents) as numpy array. From csv file it will be:

from scipy import loadtxt
data = loadtxt("/path/to/file.csv", delimiter=",")

Then you can use the most efficient way to find nearest neighbour for each point: KDTree. It's a structure that allows you to find nearest neighbours very fast. So create 2D Tree with your data:

from scipy.spatial import KDTree
tree = KDTree(data)

KDTree class has a method query_ball_point(x, r[, p, eps]) that find all points within distance r of point(s) x. So now for each point in dataset find NN in given radius (your buffer radius), and check if number of NN are greather than your threshold, if true, save this point:

for point in data:
    if len(tree.query_ball_point(point, radius)) > threshold:
        # add this point to the new layer

Now take this layer and create buffers.

I hope it will be helpful.

  • having problems installing SciPy (Windows...) but hopefully I'll be able to give this a try. does look far more efficient than my approach!
    – Alex
    Commented Jul 7, 2016 at 8:37
  • Scipy is avilable in QGIS Python console.
    – dmh126
    Commented Jul 7, 2016 at 9:28
  • finally got QGis running with SciPy, I didn't choose that package when I installed it. I used this writeAsVectorFormat(layer, r'c:\QGIS_testing\accidents.csv', "utf-8", None, "CSV", layerOptions='GEOMETRY=AS_XYZ') to export the file as CSV, this appears to have worked. Unable to open using loadtxt though IOError: [Errno 22] invalid mode ('U') or filename: 'C:\\QGIS_testing\x07ccidents.csv'. Filename is different to the input?
    – Alex
    Commented Jul 7, 2016 at 14:52
  • Maybe 'C:\\QGIS_testing\\x07ccidents.csv' ? Remember that your csv should contains only X,Y values.
    – dmh126
    Commented Jul 7, 2016 at 14:57
  • OK the csv file (when viewed in excel) had all the attributes so perhaps I didn't export correctly? I used this post for the syntax, is this not correct?
    – Alex
    Commented Jul 7, 2016 at 16:06

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.