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
3
  • 2
    I think you way is the best
    – BERA
    Oct 19, 2023 at 6: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, 2023 at 6: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, 2023 at 7:29

3 Answers 3

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).

2
  • 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, 2023 at 8:27
  • 4
    See edited answer that solves this issue
    – Babel
    Oct 19, 2023 at 8:53
4

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:

0
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)
2
  • what is the other_point_feature ? maybe you meant the point_feature ?
    – Taras
    Oct 21, 2023 at 11:53
  • @Taras Thanks, there was a typo I edited.
    – Kalak
    Oct 21, 2023 at 15:54

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.