Count points per polygon, classify point count in original layer in QGIS

Question

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.

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.

Answer 1

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

Answer 2

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.

Answer 3

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)