Teniendo en cuenta los agujeros / restricciones en la creación del polígono Voronoi en QGIS?

12

Estoy tratando de crear polígonos voronoi en QGIS que consideren "agujeros" en el dominio general. Un ejemplo sería:

ingrese la descripción de la imagen aquí

De hecho, creé el Voronois en esta imagen usando QGIS a través del comando GRASS, luego usando la herramienta "Diferencia" para crear los agujeros. Se usó un archivo de forma de polígono separado, que contiene la extensión de los agujeros, como la capa "Diferencia". Una aplicación de ejemplo sería crear polígonos alrededor de puntos de muestreo que se recopilaron entre estructuras que deberían excluirse del análisis.

Aquí surgen dos problemas:

  1. La función "diferencia" no parece funcionar al 100% correctamente, con algunos límites de polígono que se extienden hacia los "agujeros". Esto se puede solucionar encontrando una fila en la tabla de atributos que no tenga un número de identificación de polígono (o una identificación de "0").

  2. Este tipo de "perforación" después del hecho puede dar lugar a polígonos discontinuos, como lo muestra la flecha roja en la imagen.

Mi pregunta es: ¿existe una herramienta o complemento Voronoi que pueda considerar la presencia de "agujeros" en el centro del dominio, como un proceso de un solo paso, y también eliminar la generación de polígonos discontinuos? Me imagino que una herramienta así extendería un límite de polígono a la intersección más cercana con otro límite, a menos que ese otro límite golpee primero un límite de "agujero".

Cactus Inclinado
fuente
Esto sería similar pero opuesto a (creo) usar una máscara de entorno en ArcGIS . Eso le permitiría restringir los polígonos creados dentro de un límite particular. Sin embargo, no conozco ninguna herramienta que haga uso de límites / agujeros complejos (aunque tal vez en ArcGIS la máscara podría ser tan compleja; no la he probado y podría intentarlo más tarde si tengo tiempo).
Chris W
He probado la teoría ArcGIS y no funcionará. Según la pregunta vinculada, puede restringir los resultados a una forma externa. Sin embargo, los agujeros resultantes ignoran un agujero cortado en la forma. Además, si ese agujero también tiene algunos puntos, la herramienta falla y no se ejecuta. No puedo explicar su primer problema con la diferencia, pero el segundo resultado en astillas no es del todo inesperado; después de todo, esa área aún se asignaría al mismo punto, incluso si el agujero está presente. Puede usar ese método y luego incorporar las astillas en sus vecinos con un método de limpieza.
Chris W
2
Potencialmente podría resolver esto yendo a la trama. Con una máscara ráster, la distancia euclidiana saldría de sus puntos hasta que golpeara las celdas que salían de otro punto o su ráster de máscara (su descripción de golpe de límite). Luego haces una limpieza zonal y vectorizas el resultado para obtener polígonos.
Chris W
1
Me aseguraría de que la geometría voronoi sea válida ejecutando v.clean y luego verifique la geometría. Finalmente, ejecute Diferencia para crear los agujeros.
klewis
¿Qué es Voronoi acerca de estos agujeros? ¿No quieres perforar agujeros limpiamente? ¿Por qué no haría ninguna capa de polígono?
mdsumner

Respuestas:

3

Esto podría ser posible usando rásteres. Primero convierta sus puntos y polígonos de límite en un ráster de alta resolución. Establezca una máscara para sus límites usando r.mask. Luego, corre r.grow.distanceen GRASS y usa el Value= output. Esto te dará por cada píxel, que es el punto más cercano. Convierta esto nuevamente en polígonos vectoriales. Es posible que se necesiten pasos adicionales para deshacerse de los polígonos de astilla.

Guillaume
fuente
2

Esto es ciertamente posible con rásteres.

Con suerte, esta captura de pantalla muestra el problema más claramente. La parte B del voronoi está más cerca 'mientras el cuervo vuela' al centro original del voronoi, pero esto no tiene en cuenta el hecho de que tomaría más tiempo caminar por el edificio. Entiendo que la pregunta del OP es que el voronoi debe tener en cuenta esta distancia adicional para caminar por el edificio.

ingrese la descripción de la imagen aquí

Me gusta la sugerencia de @Guillaume. Sin embargo, cuando lo probé tuve problemas r.grow.distancepara honrar la máscara (ver más abajo. Las ondas no deberían pasar por los edificios).

Mi conocimiento de Grass no es tan fuerte como podría ser, así que tal vez estoy haciendo algo estúpido. Definitivamente, primero revisa esa sugerencia ya que será mucho menos trabajo que la mía ;-)

ingrese la descripción de la imagen aquí

Paso 1: crear una superficie de costo

El primer paso es crear una superficie de costo. Esto solo debe hacerse una vez.

  • crear una capa editable, agujeros y todo.
  • agregue un campo llamado 'unidad', configúrelo en 1.
  • usando polígono a ráster en su capa vectorial "perforada" (la que tiene los agujeros), usando el campo 'unidad'. Ahora tiene una capa "máscara", donde 1 es espacio libre y 0 se está construyendo.
  • use la calculadora ráster para convertir esto en una superficie de costo. Estableceré 'afuera' en 1 y 'adentro' en 9999. Esto hará que moverse por los edificios sea extremadamente difícil.

    (("máscara @ 1" = 1) * 1) + (("máscara @ 1" = 0) * 9999)

Puede obtener resultados más 'orgánicos' agregando un poco de ruido a la superficie de costo (por ejemplo, use un número aleatorio del 1 al 3, en lugar de solo 1 para pxiels al aire libre).

Paso 2. Crear rásteres de costos acumulativos para cada centro voronoi

Ahora podemos ejecutar (para una celda voronoi a la vez) el algoritmo GRASS r.cost.coordinatescontra nuestra capa de superficie de costo.

Para la coordenada de inicio, use el centro vornoi. Para la coordenada final, elija una de las esquinas de su área. Sugiero usar 'Knights Tour' ya que esto da resultados más suaves.

El resultado muestra líneas de igual tiempo de viaje desde un centro voronoi. Tenga en cuenta cómo las bandas se envuelven alrededor de los edificios.

ingrese la descripción de la imagen aquí

No estoy seguro de la mejor manera de automatizar esto. Tal vez el modo de procesamiento por lotes, o hecho en pyqgis.

Paso 3. Combina los rásteres

Esto probablemente necesitará código. El algoritmo sería

create a raster 'A' to match the size of your cumulative cost images
fill raster 'A' with a suitably high number e.g. 9999
create an array of the same size as the raster.
for each cumulative cost raster number 1..N
    for each cell in image
        if cell < value in raster 'A'
            set value in raster 'A' to cell value
            set corresponding cell in array to cum. cost image number
write out array as a raster

Ese enfoque debería generar un ráster donde cada celda está categorizada por el centro voronoi más cercano, teniendo en cuenta los obstáculos.

Entonces podría usar raster-to-polygon. A continuación, puede usar el complemento Generalizar para eliminar los artefactos de efecto "paso" del ráster.

Disculpas por la vaguedad en los pasos 2 y 3 ... Espero que alguien intervenga con una solución más elegante :)

Steven Kay
fuente
1
Gracias Steven, tengo un trabajo raster GRASS que funciona, pero esperaba una solución más elegante como se menciona en la descripción de la recompensa.
oscuro
0

Nota # 1 : No pude reproducir el problema propuesto porque la herramienta Diferencia funcionó bien para mí en varias pruebas que realicé (tal vez se debió a la geometría simple del problema o porque la herramienta se ha mejorado desde que la pregunta fue preguntado hace 1 año).

Sin embargo, propongo una solución alternativa en PyQGIS para evitar el uso de la herramienta Diferencia . Todo se basa en la intersección local entre dos capas de entrada (ver figura a continuación):

  1. una capa vectorial poligonal que representa los polígonos Voronoi;
  2. una capa vectorial poligonal que representa los agujeros / restricciones que deben excluirse del análisis.

ingrese la descripción de la imagen aquí

Nota # 2 : Como no quiero usar la herramienta Diferencia , no puedo evitar la creación de "astillas" (ver entonces), así que necesitaba ejecutar la v.cleanherramienta para eliminarlas. Además, como dijo @Chris W,

[...] pero el segundo resultado en astillas no es del todo inesperado: después de todo, esa área aún se asignaría al mismo punto, incluso si el agujero está presente. Puede usar ese método y luego incorporar las astillas en sus vecinos con un método de limpieza .

Después de estas premisas necesarias, publico mi código:

##Voronoi_Polygons=vector polygon
##Constraints=vector polygon
##Voronoi_Cleaned=output vector

from qgis.core import *

voronoi = processing.getObject(Voronoi_Polygons)
crs = voronoi.crs().toWkt()
ex = voronoi.extent()
extent = '%f,%f,%f,%f' % (ex.xMinimum(), ex.xMaximum(), ex.yMinimum(), ex.yMaximum())

constraints = processing.getObject(Constraints)

# Create the output layer
voronoi_mod = QgsVectorLayer('Polygon?crs='+ crs, 'voronoi' , 'memory')
prov = voronoi_mod.dataProvider()
fields = voronoi.pendingFields() # Fields from the input layer
prov.addAttributes(fields) # Add input layer fields to the outLayer
voronoi_mod.updateFields()

# Spatial index containing all the 'constraints'
index_builds = QgsSpatialIndex()
for feat in constraints.getFeatures():
    index_builds.insertFeature(feat)

final_geoms = {}
final_attrs = {}

for feat in voronoi.getFeatures():
    input_geom = feat.geometry()
    input_attrs = feat.attributes()
    final_geom = []
    multi_geom = input_geom.asPolygon()
    input_geoms = [] # edges of the input geometry
    for k in multi_geom:
        input_geoms.extend(k)
    final_geom.append(input_geoms)
    idsList = index_builds.intersects(input_geom.boundingBox())
    mid_geom = [] # edges of the holes/constraints
    if len(idsList) > 0:
        req = QgsFeatureRequest().setFilterFids(idsList)
        for ft in constraints.getFeatures(req):
            geom = ft.geometry()
            hole = []
            res = geom.intersection(input_geom)
            res_geom = res.asPolygon()
            for i in res_geom:
                hole.extend(i)
                mid_geom.append(hole)
        final_geom.extend(mid_geom)
    final_geoms[feat.id()] = final_geom
    final_attrs[feat.id()] = input_attrs

# Add the features to the output layer
outGeom = QgsFeature()
for key, value in final_geoms.iteritems():
    outGeom.setGeometry(QgsGeometry.fromPolygon(value))
    outGeom.setAttributes(final_attrs[key])
    prov.addFeatures([outGeom])

# Add 'voronoi_mod' to the Layers panel
QgsMapLayerRegistry.instance().addMapLayer(voronoi_mod)

# Run 'v.clean'
processing.runalg("grass7:v.clean",voronoi_mod, 2, 0.1, extent, -1, 0.0001, Voronoi_Cleaned, None)

# Remove 'voronoi_mod' to the Layers panel
QgsMapLayerRegistry.instance().removeMapLayer(voronoi_mod)

lo que lleva a este resultado:

ingrese la descripción de la imagen aquí

Solo por claridad, este sería el resultado sin el uso de la v.cleanherramienta:

ingrese la descripción de la imagen aquí

La diferencia con el resultado de @LeaningCactus es que, por ahora, las geometrías no están rotas y podrían "limpiarse" sin errores .

mgri
fuente
Haga agujeros más largos, por ejemplo, cortando todo el mapa, como un río, y verá el problema. Agregar astillas a los vecinos crea polígonos que se ven muy diferentes a un diagrama de Voronoi restringido adecuado. Lo intenté
oscuro
Lo siento, no entiendo: ¿encontró algún error en los resultados? Solo probé el código para los casos en que los polígonos eran similares a los propuestos en la pregunta.
mgri
Desafortunadamente, no puedo probar el código en este momento, pero ¿podría mostrar los resultados obtenidos con el cambio en los agujeros esbozados en i.stack.imgur.com/Jpfra.png ?
oscuro
Si extiendo la restricción hasta la característica de la derecha, obtengo esto . En cambio, si muevo directamente la restricción, obtengo esto .
mgri
El pequeño triángulo al que apunta la flecha roja en mi dibujo es el problema. No debería estar allí, pero también está en tus resultados. Parece que este enfoque resuelve el problema n. ° 1 de la pregunta, pero deja el n. ° 2 sin resolver.
oscuro