Considérant les trous / contraintes dans la création de polygones Voronoi dans QGIS?

J'essaie de créer des polygones voronoi dans QGIS qui considéreraient des "trous" dans le domaine général. Un exemple serait:

J'ai en fait créé le Voronois dans cette image en utilisant QGIS via la commande GRASS, puis en utilisant l'outil "Différence" pour créer les trous. Un fichier de formes polygonal distinct, qui contient l'étendue des trous, a été utilisé comme couche "Différence". Un exemple d'application serait de créer des polygones autour des points d'échantillonnage qui ont été collectés entre les structures qui devraient être exclues de l'analyse.

Deux problèmes se posent ici:

1. La fonction "différence" ne semble pas fonctionner à 100% correctement, avec certaines limites de polygone s'étendant dans les "trous". Cela peut être résolu en recherchant une ligne dans la table d'attributs qui n'a pas de numéro d'identification de polygone (ou d'ID "0").

2. Ce type de "perforation" après coup peut entraîner des polygones discontinus, comme le montre la flèche rouge sur l'image.

Ma question est la suivante: existe-t-il un outil ou un plugin Voronoi qui peut considérer la présence de "trous" au centre du domaine, comme un processus en une étape, et également éliminer la génération de polygones discontinus? J'imagine qu'un tel outil étendrait une frontière de polygone à l'intersection la plus proche avec une autre frontière, à moins que cette autre frontière ne claque d'abord contre une frontière de "trou".

Cela peut être possible à l'aide de rasters. Convertissez d'abord vos points et polygones de limite en un raster haute résolution. Définissez un masque pour vos limites à l'aide de r.mask. Ensuite, exécutez r.grow.distanceGRASS et utilisez le Value= output. Cela vous donnera pour chaque pixel, qui est le point le plus proche. Convertissez-le en polygones vectoriels. Il pourrait y avoir des étapes supplémentaires nécessaires pour se débarrasser des polygones de ruban.

C'est certainement possible avec les rasters.

Nous espérons que cette capture d'écran montre le problème plus clairement. La partie B du voronoi est plus proche «à vol d'oiseau» du centre d'origine du voronoi, mais cela ne tient pas compte du fait qu'il faudrait plus de temps pour se promener dans le bâtiment. Ma compréhension de la question du PO est que les voronoi doivent prendre en compte cette distance supplémentaire pour se promener dans le bâtiment.

J'aime la suggestion de @Guillaume. Cependant, quand je l'ai essayé, j'ai eu du mal r.grow.distanceà honorer le masque (voir ci-dessous. Les ondulations ne devraient pas traverser les bâtiments).

Ma connaissance de Grass n'est pas aussi solide qu'elle pourrait l'être, alors peut-être que je fais quelque chose de stupide. Certainement, consultez d'abord cette suggestion car ce sera beaucoup moins de travail que le mien ;-)

Étape 1 - Créer une surface de coût

La première étape consiste à créer une surface de coût. Cela a seulement besoin d'être fait une fois.

• créer un calque modifiable, des trous et tout.
• ajoutez un champ appelé 'unité', réglez-le sur 1.
• en utilisant polygone-à-raster sur votre couche vectorielle "perforée" (celle qui a les trous), en utilisant le champ "unité". Vous avez maintenant un "masque" de calque, où 1 est un espace libre et 0 est un bâtiment.

• utilisez la calculatrice raster pour en faire une surface de coût. Je définirai «extérieur» à 1 et «intérieur» à 9999. Cela rendra les déplacements dans les bâtiments extrêmement difficiles.

  (("masque @ 1" = 1) * 1) + (("masque @ 1" = 0) * 9999)

Vous pouvez obtenir des résultats plus «organiques» en ajoutant un peu de bruit à la surface des coûts (par exemple, utilisez un nombre aléatoire de 1 à 3, plutôt que seulement 1 pour les pxiels extérieurs.)

Étape 2. Créez des rasters de coûts cumulés pour chaque centre de voronoi

Maintenant, nous pouvons exécuter (pour une cellule de voronoï à la fois) l'algorithme GRASS r.cost.coordinatescontre notre couche de surface de coût.

Pour la coordonnée de départ, utilisez le centre de vornoi. Pour la coordonnée de fin, choisissez l'un des coins de votre zone. Je suggère d'utiliser "Knights Tour" car cela donne des résultats plus fluides.

Le résultat montre des lignes de temps de trajet égal à partir d'un centre voronoi. Notez comment les bandes s'enroulent autour des bâtiments.

Je ne sais pas comment automatiser cela au mieux. Peut-être le traitement en mode batch, ou fait dans pyqgis.

Étape 3. Fusionnez les rasters

Cela aura probablement besoin de code. L'algorithme serait

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

Cette approche devrait produire un raster où chaque cellule est classée par le centre de voronoï dont elle est la plus proche, en tenant compte des obstacles.

Vous pouvez ensuite utiliser le tramage sur polygone. Vous pouvez ensuite utiliser le plug-in Généraliser pour supprimer les artefacts d'effet "pas" du raster.

Toutes mes excuses pour l'imprécision des étapes 2 et 3 ... j'espère que quelqu'un interviendra avec une solution plus élégante :)

Note # 1 : Je n'ai pas pu reproduire le problème proposé car l' outil Différence a bien fonctionné pour moi dans plusieurs tests que j'ai effectués (peut-être à cause de la géométrie simple du problème ou parce que l'outil a été amélioré depuis que la question était demandé il y a 1 an).

Cependant, je propose une solution de contournement dans PyQGIS pour éviter l'utilisation de l' outil Différence . Tout est basé sur l'intersection locale entre deux couches d'entrée (voir figure ci-dessous):

1. une couche vectorielle polygonale représentant les polygones de Voronoi;
2. une couche vectorielle polygonale représentant les trous / contraintes qui doivent être exclus de l'analyse.

Note # 2 : Comme je ne veux pas utiliser l' outil Différence , je ne peux pas éviter la création de "slivoïdes" (voir alors), j'ai donc dû exécuter l' v.cleanoutil pour les éliminer. En outre, comme l'a dit @Chris W,

[...] mais le second résultant en éclats n'est pas entièrement inattendu - après tout, cette zone serait toujours allouée au même point même si le trou est présent. Vous pouvez utiliser cette méthode, puis incorporer les éclats à leurs voisins avec une méthode de nettoyage .

Après ces locaux nécessaires, je poste mon code:

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

ce qui conduit à ce résultat:

Pour plus de clarté, ce serait le résultat sans l'utilisation de l' v.cleanoutil:

La différence avec le résultat de @LeaningCactus est que, maintenant, les géométries ne sont pas cassées et peuvent être "nettoyées" sans erreur .