Agrupando valores por tipo y contándolos

Question

Tengo un conjunto de puntos en un archivo csv y un shapefile, y utilicé los buffers de los puntos para obtener el subconjunto de polígonos del shapefile que intersectan los buffers.

Cada uno de los polígonos tiene un tipo específico, digamos A, B, C, D.

Lo que necesito saber es para cada punto, cuáles son los tipos de polígonos que intersectan su buffer.
Desafortunadamente, en algunos casos, un solo punto tiene múltiples polígonos en su buffer que tienen el mismo tipo, pero necesito contarlos solo una vez porque al final, me gustaría saber cuántas ocurrencias hay de cada tipo en total.
Entonces, por ejemplo, digamos que tengo 4 puntos, cada uno con un subconjunto de polígonos

Punto 1: A, B, C   
Punto 2: B, C, C, C, C, D   
Punto 3: A, A, A, B, B   
Punto 4: C, C, C, D, D   

Me gustaría editar estos datos para que se vean así

Punto 1: A, B, C  
Punto 2: B, C, D  
Punto 3: A, B   
Punto 4: C, D  

Después de hacer eso, me gustaría contar los tipos totales, para saber que en total hay 2 ocurrencias de A, 3 de B, 3 de C y 2 de D.
¿Cómo puedo hacer algo así?

Answer 1

La siguiente solución produce este resultado: Teniendo 1 capa de polígonos con 3 tipos (a:rojo, b:verde, c:azul) y una capa de puntos con 4 puntos (nr 1-4). Utilizando una capa virtual en QGIS con la siguiente definición se obtiene una capa virtual que muestra las intersecciones de cada punto bufferizado con los tipos de polígonos (¡el bufferizado y la intersección se realizan en una sola operación!). Obtendrás una nueva característica para cada intersección con un tipo de polígono y una columna con el número de conteos. El punto bufferado 1 interseca dos veces el tipo "a" y también "b" y "c". La definición para esa capa virtual es la siguiente:

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

Answer 2

Esta es una solución en Python. Debes pegar este código en la consola de Python de QGIS.

Este es mi setup: Tengo dos capas vectoriales GIS_SE_POINT y GIS_SE_POLY cada una con un campo (atributo) llamado name (para describir un punto como Punto 1 y un polígono como A).

Canvas:

Este es el código:

EDITADO

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

Cómo obtener el resultado:

Debes llamar a la función con 6 parámetros.

1. nombre de tu capa de puntos
2. campo (atributo) con los nombres de los puntos
3. nombre de la capa de polígonos
4. campo (atributo) con los nombres de los polígonos
5. tamaño del buffer en metros (aumenta/disminuye el valor según desees)
6. segmentos del buffer (aumentar la cantidad de segmentos aumenta el tiempo de procesamiento)
7. alternando el estado: 0 sin agrupar ni contar, 1 agrupando, 2 agrupando y contando.

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

---

ACTUALIZACIÓN

Para una pequeña cantidad de datos, el resultado del enfoque anterior es razonable en términos de tiempo. Para conjuntos de datos más grandes, el siguiente script utiliza un índice espacial para los polígonos. -- Atención -- Probé varias capas y noté que si la capa de polígonos no es un shapefile, no se puede establecer el índice espacial.

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

Con el índice espacial realmente puedes acelerar tu tiempo de procesamiento.

Caso de prueba:

Capa de puntos con 2250 elementos y una capa de polígonos con 600 elementos. Pruebas del tiempo de procesamiento sin imprimir.

Sin usar un índice espacial: 65seg

Usando un índice espacial: 0.35seg

¡Más de 180 veces más rápido!