5

I have a 'POINT' layer named export_bdd with a field called denom. Many of the points overlap because they are saved at the same coordinates.

For each record where denom = 'Teuc pseu', I want to search within a 10-meter buffer zone for the number of entities where denom = 'Teuc pseu'.

To do this, I tried the following expression before applying any filters:

CASE WHEN denom = 'Teuc pseu' THEN array_length( overlay_intersects(@layer, buffer($geometry, 10))) ELSE 0 END

However, the results are always incorrect. The number of entities found is not right, and it doesn't change regardless of the buffer size (whether 0.0001 or 1,000,000 meters, the result is the same).

Do you know what might be going wrong?

Example:

I choose a random point :

random_point

You can see in the last field on the right, the result of my expression above returned 2. On the picture, we can see that more features are located within the 10m. And here is the attribute table :

Attribute table - selected features

My expression should return 30 instead of 2 (before I apply a filter on denom).

Here is also a link to download the layer: Test layer

I tried with a new project, with correct CRS.

Improve this question
3
  • 1
    You can't use overlay_intersects in the way you want. Your expression evaluates intersection (without buffer), then only for the results creates a buffer as output. For what you want to do, use overlay_nearestwith max_distance Commented Jan 22 at 19:00
  • Why do you use a CRS that is not valid in your region of interest? Commented Jan 22 at 20:55
  • 1
    I relocated the entities to another country. The data isn’t confidential but is still private, so I modified both the content and the location. However, I didn’t consider the CRS while making these changes. Commented Jan 23 at 2:14

3 Answers 3

Reset to default
2

Try aggregate

Returns an aggregate value calculated using features from another layer.

(It can be the same layer)

I have a point layer named building with an integer field named kkod. For each point I want to count other points within 1000 m sharing the same kkod, but with a different id ogc_fid (to not count the point itself). Replace ogc_fid with the name of your unique id column

aggregate( layer:='building', aggregate:='count', expression:=1, filter:= "kkod"=attribute(@parent, 'kkod') AND distance(geometry(@parent), $geometry)<1000 AND "ogc_fid"<>attribute(@parent, 'ogc_fid') )

I have labeled the points using the expression above. The circle is a 1000 m buffer: enter image description here

Improve this answer
1
  • 1
    I've tried this solution, and it worked. I like it because in my case, this is a one time problem, and this expression is solving it easily. Commented Jan 23 at 2:56
2

There seem to be several problems.

  • You can't use function overlay_intersects in the way you want. Your expression evaluates intersection (without buffer) with the current feature's geometry, so here returning points that overlap 100% with the current point, and then only for these results creates a buffer as output. You in fact count point clusters (overlaps). To state it otherwise: this function does not use buffers (or distance from current feature) to evaluate what intersects. The buffer is the result, created around the geometries (points) that intersect the current feature's geometry. For what you want to do, however, use overlay_nearest with max_distance.

  • Overlay functions seem not to work properly when applied to the current layer (@layer). To avoid problems, create a duplicate of your layer and use the expression referring to this layer.

  • You use a CRS that is not valid in your region of interest. Reproject the layer (do not set CRS).

Use this expression. Be aware that the current point is included in the count, so you might want to subtract it from the result (-1):

array_length( overlay_nearest ( 'duplicate', -- name of duplicated layer @id, limit:=-1, max_distance:=10, filter:=denom = 'Teuc pseu' ) )

Point clusters (showing how many points are in the same place), including the denom values as label in black; the expression above is used to create the red label. For the point in the middle, the black circle shows the 10 m distance (buffer): enter image description here

Improve this answer
1
  • A bit similar to the one proposed by @Bera, it worked with little adjustment to meet specificities of my project. Commented Jan 23 at 2:57
2

Imho it's a bad approach to use buffer+intersection for pure point-distance-queries, because the buffer-algorithm is very expensive and doesn't produce mathematical circles but circle-approximated polygons. Below a script for your needs, which does this calculation via distance and writes the results to your layer if in edit-mode. In an extra string-field it additionally writes the fid's of the features in 10m distance. Copy&paste to Python-Console, adapt to your needs, set your layer active and editable, and run the script.

## settings begin # field containing the aggregate-information aggregate_field = 'denom' # narrowing down the scope, if empty, the scripts checks all features filter_value = 'Teuc pseu' # find neighbours having the same aggregate-attribute within this distance check_distance = 10 # target-field which is filled with the number of the neighbours, type int num_neighbours_field = 'test' # target-field filled with the concatenated neighbour-feature-ids # type string (stringified list of integers, sufficient size required!) # optional, disregarded if empty neighbour_ids_field = 'sub_fids' ## settings end vl = iface.activeLayer() request = QgsFeatureRequest() # only filter if required if filter_value: request.setFilterExpression(f"{aggregate_field}='{filter_value}'") for feat in vl.getFeatures(request): sub_request = QgsFeatureRequest() # sub-request using the queried denom from parent iteration, the parent feature-id is excluded sub_request.setFilterExpression(f"{aggregate_field}='{feat['denom']}' and $id!={feat.id()}") sub_request.setDistanceWithin(feat.geometry(), check_distance) sfc = 0 sub_fids = [] for sub_feat in vl.getFeatures(sub_request): sfc += 1 sub_fids.append(sub_feat.id()) # dist = feat.geometry().distance(sub_feat.geometry()) # print("sub_feat", sub_feat.id(), sub_feat['denom'], dist) # print("feat", feat.id(), feat['denom'], sfc, sub_fids) if vl.isEditable(): feat[num_neighbours_field] = sfc if neighbour_ids_field: feat[neighbour_ids_field] = str(sub_fids) vl.updateFeature(feat)
Improve this answer
3
  • I don’t typically work with Python, but I decided to give it a try in the Python console. The layer updates correctly, behaving the same way as expressions. However, I encountered an error message about a wrong data type (QString). Since I’m using integers, I don’t quite understand the issue. Additionally, there’s a strange behavior with the fid (or sub_fids) field. Initially, it contains unique values, but after execution, I end up with array values like '[83, 84, 154, 156, 329, 332, 339, 342, 346, 367]'. Commented Jan 23 at 3:15
  • Edit: It does what you said actually. This is not en issue. Commented Jan 23 at 7:14
  • Glad I could help. I expended and commented the script. Commented Jan 24 at 12:41

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.