3

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?

2
  • 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.
    – Stefan
    Feb 20, 2018 at 13:57

2 Answers 2

1

The following solution produces this result: enter image description here 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
1

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

enter image description here

enter image description here

Canvas:

enter image description here

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!

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.