Cluster analysis using buffers and polygon features (NOT layer)

Question

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:
              my_clusters.append(area.id())

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!

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.