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.

• Alex: having problems installing SciPy (Windows...) but hopefully I'll be able to give this a try. does look far more efficient than my approach!
• dmh126: Scipy is avilable in QGIS Python console.
• Alex: 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?
• dmh126: Maybe 'C:\\QGIS_testing\\x07ccidents.csv' ? Remember that your csv should contains only X,Y values.
• Alex: 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?