2

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!

Improve this question

1 Answer 1

Reset to default
3

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.

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

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.