4

I have a question that I believe is similar to Count points in polygon grouping by attribute using PyQGIS but not really the same, as I would like to do it without PyQGIS and also have the information added in one of the original layers, instead of a new one.

In my case, I have a point-layer which contains several points with a lot of data (columns). There is one column that is most significant, as I need the points classified according to that column.

enter image description here enter image description here

I also have a polygon-layer containing a buffer.

I now wish to count the points within each polygon and I want the count and classification to be added (in the original buffer-layer if possible), so that each class gets its own column, which then contains the count for this particular class in the buffer.

If the buffers overlap and the point falls into several buffers, the point should count for each buffer.


EDIT: What I meant with the last statement: When 1 point falls into 3 buffer-polygons for example, I would like each polygon to recognise the point and note down "1", instead of perhaps just couting it once for the first polygon and then ignoring it for the other polygons, or counting it thrice for one polygon and not at all for the other overlapping polygons.

enter image description here

0

3 Answers 3

4

You can achieve this using a virtual layer.

Go to the menu Layer > Add Layer > Add/Edit Virtual Layer... and enter the following query. Modify pt.cat for pt.nutzart (the name of your category field) as well as the name of the layers. It assumes the polygon primary key is called id (maybe change it to oid)

This query will keep all fields from the polygon and count the points for each specified category. You can then persist the layer if you wish.

select p.*, 
    count(case when pt.cat='category A' then 1 else null end) cnt_a,
    count(case when pt.cat='category B' then 1 else null end) cnt_b,
    count(case when pt.cat='category C' then 1 else null end) cnt_c,
    count(case when pt.cat='category D' then 1 else null end) cnt_d
from myPolygonLayer p
join myPointLayer pt
 on st_intersects(p.geometry, pt.geometry)
group by p.id
3
  • One correction on my part. The relevant column isn't 'nutzart' but 'neue_bez'. I have tried your code and ran into two problems. One seems to have to do with the '/' in the class name. But when I remove that, I get a very long error message, stating the layer is invalid because: (...) 1 - near "from": syntax error".
    – gruenbet
    Aug 29, 2022 at 7:23
  • make sure the statement just before the from is properly written (including no comma at the end).
    – JGH
    Aug 29, 2022 at 13:42
  • I have corrected my mistake and the code went through! Sadly, the count is off as well whenever polygons overlap. I have updated my original statement regarding the specifications of counting as it might have been misleading. Any idea how to archieve that way of counting or why the counting is off right now?
    – gruenbet
    Aug 30, 2022 at 6:54
4

Using field calculator with expressions you can create a new integer field on your polygonlayer with the expression:

array_length(
  overlay_intersects(
    'pointlayername',
    $id,
    filter:="fieldcontainingcategory"='nameofcategory'
  )
)

Repeat this as many times as you want the count of different categories.

Instead of overlay_intersects() you can also use overlay_contains() depending on the method you wish.

5
  • I tried your code for Overlay_intersects() and although it goes through without error messages, the count is off when buffers overlay. From looking at some example areas, it seems that when a points falls into three buffer-polygons for example, the point is counted three times for one polygon.
    – gruenbet
    Aug 29, 2022 at 7:43
  • Well you said "If the buffers overlap and the point falls into several buffers, the point should count for each buffer", thats what it does.
    – MrXsquared
    Aug 29, 2022 at 7:47
  • I see now how that is a misleading statement of mine. What I meant is simply that each buffer recognises this point and counts it, instead of only counting it once for the first polygonlayer and then ignoring it in every other overlaying polygon. So, let's say a point falls into 3 overlaying buffer-polygons. Then I would like each polygon to note down "1".
    – gruenbet
    Aug 29, 2022 at 7:56
  • Is there a way to eliminate the double count? I have edited my original question to make my statement/request clearer. I have since tried other overlay-functions (touches, within and crosses) which usually counted 0.
    – gruenbet
    Aug 29, 2022 at 11:56
  • I cannot think of a way fulfilling this condition using expressions. But you can do this with PyQGIS, just as I did in ProcessX Plugin's algorithm Count Points In Polygon With Condition. Using the Plugin you need to run the algorithm as often as you have different categories.
    – MrXsquared
    Sep 3, 2022 at 11:40
1

Here is a pretty quick PyQGIS solution for your updated question. It adds x new fields to the polygon layer naming them count_nameofcategory and adds the count of each category. You can also set that a point is only counted once, if so, it will only be counted for the first polygonfeature ordered by featureid. You can implement an iteration order if you want to change that. Note that it uses QgsSpatialIndexKDBush, so it will only work with 2D-SinglePoints. I am aware the script is not perfect and can be improved, but for now just change the first 4 lines to your needs and run it:

### Settings:
polygonlayer = iface.activeLayer()
pointlayer = QgsProject.instance().mapLayersByName('points')[0]
fieldwithcategories = 'category'
countpointonlyonce = True

### no changes needed ###
attrvaluesdict = {}
categories = []
req = QgsFeatureRequest().setFlags(QgsFeatureRequest.NoGeometry).setSubsetOfAttributes([fieldwithcategories], pointlayer.fields())
for feat in pointlayer.getFeatures(req):
    categories.append(feat[fieldwithcategories])
    attrvaluesdict[feat.id()] = feat[fieldwithcategories]
categories = set(categories)

pointlayer_idx = QgsSpatialIndexKDBush(pointlayer.getFeatures())
pointlayerskip = []

with edit(polygonlayer):
    for cat in categories:
        polygonlayer.dataProvider().addAttributes([QgsField('count_'+str(cat), QVariant.Int)])
        polygonlayer.updateFields()
    
    for feat in polygonlayer.getFeatures():
        for cat in categories:
            feat['count_'+str(cat)] = 0
            polygonlayer.updateFeature(feat)
        featgeom = feat.geometry()
        featgeomengine = QgsGeometry.createGeometryEngine(featgeom.constGet())
        featgeomengine.prepareGeometry()
        intersectingpoints = pointlayer_idx.intersects(feat.geometry().boundingBox())
        for ipoint in intersectingpoints:
            if ipoint.id in pointlayerskip and countpointonlyonce is True:
                continue
            pointgeom = QgsGeometry.fromPointXY(ipoint.point()).constGet()
            if featgeomengine.intersects(pointgeom):
                pointlayerskip.append(ipoint.id)
                category = attrvaluesdict[ipoint.id]
                feat['count_'+str(category)] += 1
                polygonlayer.updateFeature(feat)

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.