# Grouping values by type and counting them

I have a set of points in a csv and a shapefile, and I used the buffers of the points to obtain the subset of polygons from the shapefile that intersects the buffers.

Each of the polygons has a specific type, let's say A,B,C,D.

What I need to know is for each point, what are the types of the polygons that intersect its buffer.
Unfortunately, in some cases, a single point has multiple polygons in its buffer that have the same type but I need to count them only once because in the end, I would like to know how many occurrences of each type there are in total.
So, for example, let's say that I have 4 points, each with a subset of polygons

``````Point 1: A,B,C
Point 2: B,C,C,C,C,D
Point 3: A,A,A,B,B
Point 4: C,C,C,D,D
``````

I would like to edit this data to look like this

``````Point 1: A,B,C
Point 2: B,C,D
Point 3: A,B
Point 4: C,D
``````

After doing that I would like to count the total types, in order to know that there are in total 2 occurrences of A, 3 of B, 3 of C and 2 of D.
How can I make something like that?

• Hi, can you please send the data? Feb 16, 2018 at 9:05
• I did a few edits and one bigger update on my answer. With the use of a spatial index you can really speed up the processing time. Feb 20, 2018 at 13:57

The following solution produces this result: Having 1 polygonlayer with 3 types (a:red,b:green,c:blue) and one point layer with 4 points (nr 1-4). Using a virtual layer in QGIS with the following definition produces a virtual layer showing the intersects of each buffered point with the polygon types (buffering and intersecting is done in one operation!). You will get a new feature for each intersect with a polygontype and a column with the number of counts. bufferd point 1 intesects twice the type "a" and also "b" and "c". the Definition for that virtual layer is like that:

``````select st_buffer(punkte.geometry, 1000) as buffer, punkte.nr,
polygon.name, count(polygon.name) as anzahl
from punkte left join polygon on
st_intersects(st_buffer(punkte.geometry, 1000), polygon.geometry)
group by punkte.nr, polygon.name
``````

This is a Python solution. You have to paste this code into the QGIS Python console.

This is my setup: I have two vector layer `GIS_SE_POINT` and `GIS_SE_POLY` each with a field (attribute) called `name` (to describe a point like `Point 1` and a polygon like `A`).

Canvas:

This is the code:

EDITED

``````import collections
from operator import itemgetter
from itertools import groupby

def grouping_intersection(point_layer, point_attribute, polygon_layer, group_by_attribute, buffer, segments, state):
if QgsMapLayerRegistry.instance().mapLayersByName(point_layer):
point = QgsMapLayerRegistry.instance().mapLayersByName(point_layer)[0]
else:
print('{}{}{}'.format('Layer ', point_layer, ' does not exist.'))
return False
index = point.fieldNameIndex(point_attribute)
if index == -1:
print('{}{}{}'.format('Field ', point_attribute, ' does not exist.'))
return False
else:
pass
if QgsMapLayerRegistry.instance().mapLayersByName(polygon_layer):
poly = QgsMapLayerRegistry.instance().mapLayersByName(polygon_layer)[0]
else:
print('{}{}{}'.format('Layer ', polygon_layer, ' does not exist.'))
return False
index = poly.fieldNameIndex(group_by_attribute)
if index == -1:
print('{}{}{}'.format('Field ', group_by_attribute, ' does not exist.'))
return False
else:
pass
list = []
for b in poly.getFeatures():
for a in point.getFeatures():
if buffer >0:
if a.geometry().buffer(buffer, segments).intersects(b.geometry()):
list.append([a[point_attribute], b[group_by_attribute]])
elif buffer == 0:
if a.geometry().intersects(b.geometry()):
list.append([a[point_attribute], b[group_by_attribute]])
for points, polygons in groupby(sorted(list, key=itemgetter(0)), itemgetter(0)):
item_list=[]
for i in polygons:
item_list.append((i)[1])
counter=collections.Counter(item_list)
if state == 0:
print('{}{}{}'.format(points, ': ', ','.join(sorted(item_list))))
elif state == 1:
print('{}{}{}'.format(points, ': ', ','.join(sorted(counter.keys()))))
elif state == 2:
print('{}{}{}'.format(points, ': ', ','.join('{}{}{}{}'.format(key, '(', val, ')') for key, val in sorted(counter.items()))))
else:
print('{}{}{}'.format('State ', state, ' does not exist'))
return False
``````

How to get the result:

You have to call the function with 6 parameters.

1. name of your point layer
2. field (attribute) with the names of the points
3. name of the polygon layer
4. field (attribute) with the names of the polygons
5. buffer in meter (increase/decrease the value as you like)
6. buffer segments (increasing the amount of segments increases the processing time)
7. toggling the state: `0` no grouping and no counting, `1` grouping, `2` grouping and counting.
``````>>>grouping_intersection('GIS_SE_POINT', 'name', 'GIS_SE_POLY', 'name', 1000, 20, 0)
Point 1: A,A,A,B,C
Point 2: A,A,B,C,D
Point 3: B,B
Point 4: A,A,B,C,C,D,D

>>>grouping_intersection('GIS_SE_POINT', 'name', 'GIS_SE_POLY', 'name', 1000, 20, 1)
Point 1: A,B,C
Point 2: A,B,C,D
Point 3: B
Point 4: A,B,C,D

>>>grouping_intersection('GIS_SE_POINT', 'name', 'GIS_SE_POLY', 'name', 1000, 20, 2)
Point 1: A(3),B(1),C(1)
Point 2: A(2),B(1),C(1),D(1)
Point 3: B(2)
Point 4: A(2),B(1),C(2),D(2)
``````

UPDATE

For a small amount of data the result of the approach above yields in reasonable time. For larger data sets the following script uses a spatial index for the polygons. -- Attention -- I tried a few layer and recognized that if the polygon layer is not a shapefile the spatial index cannot be set.

``````import collections
from operator import itemgetter
from itertools import groupby

def grouping_intersection_spIndex(point_layer, point_attribute, polygon_layer, group_by_attribute, buffer, segments, state):
if QgsMapLayerRegistry.instance().mapLayersByName(point_layer):
point = QgsMapLayerRegistry.instance().mapLayersByName(point_layer)[0]
else:
print('{}{}{}'.format('Layer ', point_layer, ' does not exist.'))
return False
index = point.fieldNameIndex(point_attribute)
if index == -1:
print('{}{}{}'.format('Field ', point_attribute, ' does not exist.'))
return False
else:
pass
if QgsMapLayerRegistry.instance().mapLayersByName(polygon_layer):
poly = QgsMapLayerRegistry.instance().mapLayersByName(polygon_layer)[0]
else:
print('{}{}{}'.format('Layer ', polygon_layer, ' does not exist.'))
return False
index = poly.fieldNameIndex(group_by_attribute)
if index == -1:
print('{}{}{}'.format('Field ', group_by_attribute, ' does not exist.'))
return False
else:
pass
points = [feature for feature in point.getFeatures()]
polygons = [feature for feature in poly.getFeatures()]
polygon_spIndex = QgsSpatialIndex()
for feat in polygons:
polygon_spIndex.insertFeature(feat)
list = []
for a in points:
if buffer >0:
pt = a.geometry().buffer(buffer,segments)
if buffer == 0:
pt = a.geometry()
for id in polygon_spIndex.intersects(pt.boundingBox()):
if pt.intersects(polygons[id].geometry()):
list.append([a[point_attribute], polygons[id][group_by_attribute]])
for points, polygons in groupby(sorted(list, key=itemgetter(0)), itemgetter(0)):
item_list=[]
for i in polygons:
item_list.append((i)[1])
counter=collections.Counter(item_list)
if state == 0:
print('{}{}{}'.format(points, ': ', ','.join(sorted(item_list))))
elif state == 1:
print('{}{}{}'.format(points, ': ', ','.join(sorted(counter.keys()))))
elif state == 2:
print('{}{}{}'.format(points, ': ', ','.join('{}{}{}{}'.format(key, '(', val, ')') for key, val in sorted(counter.items()))))
else:
print('{}{}{}'.format('State ', state, ' does not exist'))
return False
``````

With the spatial index you can really speed up your processing time.

Test case:

Point layer with 2250 features and a polygon layer with 600 features. Testing the processing time without printing.

Without using a spatial index: 65sec

Using a spatial index: 0.35sec

More than 180 times faster!