# Calculating point density as attribute in point layer in QGIS

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
• I think you way is the best
– Bera
Commented Oct 19, 2023 at 6:29
• 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` Commented Oct 19, 2023 at 6:44
• 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
Commented Oct 19, 2023 at 7:29

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

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

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.

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:

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

References:

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