Polígonos de Thiessen con el algoritmo de Fortune en PyQGIS

En el post pasado se consideró un posible bug de QGIS que no hacía posible obtener los polígonos de Thiessen (o Voronoi como también se les conoce) cuando el número de puntos de control era bajo. En la búsqueda de soluciones se corroboró que ArcMap no presentaba el problema y también se encontró una opción programada para Python que empleaba el algoritmo de Fortune.

Tal como se observa en la imagen siguiente, la opción en Python, ejecutada a través de la consola de bash, imprime los resultados en una ventana programada con el módulo ‘Tkinter’ después de hacer doble click sobre puntos arbitrarios del Canvas (para “marcarlos”) y click sencillo sobre el botón de ‘Calculate’.

vorono1

En la imagen anterior también se observa que en la consola de bash se imprimen todas las líneas que conforman los “polígonos”; que más bien constituyen figuras abiertas que se extienden hacia el infinito. Generalmente, el software SIG se encarga de delimitar estos polígonos exteriores a través de una Bounding Box arbitraria.

Como a mi me interesa que los resultados se impriman en el Map Canvas de QGIS tuve que encontrar primero el nombre de las variables a usar y como pasarle los parámetros a la función (‘Voronoi’). Como todo estaba muy bien estructurado no fue difícil. La variable era ‘points’ (una lista de puntos que puede obtenerse de los features de una capa) que imprimía una lista de segmentos de recta (sólo 2 puntos) denominada ‘lines’. Ahora sólo es cuestión de adaptar el formato de las líneas que aquí se imprime al de QgsPoint de PyQGIS para conformar la respectiva memory layer. El código completo se encuentra a continuación:

from qgis.core import QgsPoint, QgsVectorLayer, QgsFeature, QgsGeometry, QgsMapLayerRegistry
from Voronoi import Voronoi

layer = iface.activeLayer()

points = [ feat.geometry().asPoint() for feat in layer.getFeatures() ]

vp = Voronoi(points)
vp.process()
lines = vp.get_output()

qgis_lines = []

for line in lines:
    qgis_lines.append([QgsPoint(line[0], line[1]), QgsPoint(line[2], line[3])]) 

print qgis_lines

crs = layer.crs()
epsg = crs.postgisSrid()

uri = "LineString?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                          'line',
                          'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(qgis_lines)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPolyline(qgis_lines[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

Al ejecutarlo en la Consola de Python, para la misma capa de tres puntos ya considerada, se obtiene lo siguiente:

voronoi1

Incorporando 5 puntos más al shapefile y ejecutando nuevamente:

voronoi2

el resultado se obtuvo de manera rápida y eficiente. Sólo faltaría la ‘Bounding Box’ que se deja como tarea pendiente.

Esta entrada fue publicada en PyQGIS, QGIS, SIG, Software Libre. Guarda el enlace permanente.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s