6

For each point in my layer, I would like to calculate the number of other points within x distance, and for that information to be in the attribute table.

I have a method but it is a few steps. As it seems like a basic but common issue, I was just wondering if there was a one-step method/tool I am missing in QGIS.

My current steps are:

  1. Create circular buffers around all points at x distance
  2. Count points in polygon
  3. Spatial join resulting counts polygons back to the original points layer
  4. extra step - field calculator to get counts/buffer area to get points/km2
Babel
  • 71,072
  • 14
  • 78
  • 208
AdSad
  • 163
  • 6
  • 2
    I think you way is the best – BERA Oct 19 '23 at 06:29
  • 2
    The only "one step" approach I can see would be a Python script, which in this case especially would only be a few lines using QgsFeatureRequest.setDistanceWithin – Kalak Oct 19 '23 at 06:44
  • 1
    I would rather join the count using the ID, otherwise you'd join the number of points within a polygon to all points within said polygon. Otherwise this works fine. Alternatively you could use heatmap kernel density estimation and add the density data back to your points. – Erik Oct 19 '23 at 07:29

3 Answers3

4

The best one-step method I can think of is using the Field Calculator to directly create a new attribute with this expression:

array_length(
    overlay_nearest(
        @layer, 
        $id,
        max_distance:=20, -- set distance here
        limit:=count($id)
    )
)

Remark: Unfortunately, the limit argument doesn't work using -1 when applied to the same layer (returning zero instead of all matching features). To address this, simply use the number of points contained in the layer with count($id).

Taras
  • 32,823
  • 4
  • 66
  • 137
Babel
  • 71,072
  • 14
  • 78
  • 208
  • 1
    This seems a bit of a weak best one-step solution if either it needs 2 steps or it forces you to set a parameter that depending on your data could give you an inaccurate result. – Kalak Oct 19 '23 at 08:27
  • 4
    See edited answer that solves this issue – Babel Oct 19 '23 at 08:53
3

There is a possibility of using a "Virtual Layer" through Layer > Add Layer > Add/Edit Virtual Layer....

Let's assume there is a point layer called 'points_in_polygon', see the image below.

input

With the following query, it is possible to calculate the number of other points within x distance, and for that information to be in the attribute table.

SELECT p1.*, COUNT(p2.id) AS WITHIN50
FROM "points_in_polygon" AS p1
LEFT JOIN "points_in_polygon" AS p2 ON
      NOT ST_Equals(p1.geometry, p2.geometry)
      AND ST_Distance(p1.geometry, p2.geometry) <= 50 -- this value can be changed
GROUP BY p1.id

It encapsulates two functions i.e. ST_Distance and ST_Equals.

Unfortunatelly the ST_DWithin function is not implemented in Virtual Layers, however, another solution will be available through the SpatialLite PtDistWithin function, as was suggested by @BERA.

SELECT p1.*, COUNT(p2.id) AS WITHIN50
FROM "points_in_polygon" AS p1
LEFT JOIN "points_in_polygon" AS p2 ON
      NOT ST_Equals(p1.geometry, p2.geometry)
      AND PtDistWithin(p1.geometry, p2.geometry, 50.0) -- this value can be changed
GROUP BY p1.id

In both cases, the output point layer with its attribute table will look like:

result

Note: More sophisticated solutions can be performed using the PostGIS, see the first reference.

References:

Taras
  • 32,823
  • 4
  • 66
  • 137
1

One-step solution using PyQGIS script:

from qgis.core import QgsProject, QgsFeatureRequest

"points" is the name of the layer in legend


layer_name = "points"


Distance within which to look for neighbors


distance = 100


Field name (already in layer) into which the result will be saved


field_name = "field_x"


point_layer = QgsProject.instance().mapLayersByName(layer_name)[0]
idx_field_to_change = point_layer.fields().indexFromName(field_name)


dict_attribute_change = {}


for point_feature in point_layer.getFeatures():
    point_geometry = point_feature.geometry()
    request = QgsFeatureRequest().setDistanceWithin(point_geometry, distance)
    dict_attribute_change[point_feature.id()] = {idx_field_to_change: 0}


for other_point_feature in point_layer.getFeatures(request):

    # comment the next two lines if you want to count the origin point 
    # as being a neighbor of itself
    if other_point_feature.id() == point_feature.id():
        continue

    dict_attribute_change[point_feature.id()][idx_field_to_change] += 1




point_layer.dataProvider().changeAttributeValues(dict_attribute_change)
Kalak
  • 3,848
  • 8
  • 26