Considering holes/constraints in Voronoi polygon creation in QGIS?

Question

I am trying to create voronoi polygons in QGIS that would consider "holes" in the general domain. An example would be:

I actually created the Voronois in this image using QGIS through the GRASS command, then using the "Difference" tool to create the holes. A separate polygon shapefile, which contains the extents of the holes, was used as the "Difference" layer. An example application would be creating polygons around sampling points that were collected between structures that should be excluded from the analysis.

Two problems arise here:

1. The "difference" function does not appear to work 100% properly, with some polygon boundaries extending into the "holes". This can be fixed by finding row in the Attribute Table which do not have a polygon ID number (or ID of "0").

2. This type of after-the-fact "hole-punching" can result in discontinuous polygons, as shown by the red arrow in the image.

My question is: is there a Voronoi tool or plugin that can consider the presence of "holes" in the center of the domain, as a one-step process, and also eliminate the generation of discontinuous polygons? I envision that such a tool would extend a polygon boundary to the nearest intersection with another boundary, unless that other boundary slams against a "hole" boundary first.

Answer 1

This might be possible using rasters. First convert your points and boundary polygons into a high resolution raster. Set a mask for your boundaries using r.mask. Then, run r.grow.distance in GRASS and use the Value= output. This will give you for each pixel, which is the closest point. Convert this back into vector polygons. There might be extra steps needed to get rid of sliver polygons.

Answer 2

This is certainly possible with rasters.

This screenshot hopefully shows the issue more clearly. The portion B of the voronoi is closer 'as the crow flies' to the original voronoi centre, but this doesn't take into account the fact that it would take longer to walk around the building. My understanding of the OP's question is that the voronoi needs to take into account this extra distance to walk around the building.

I like the suggestion from @Guillaume. However, when I tried it I had problems getting r.grow.distance to honour the mask (see below. The ripples shouldn't pass through the buildings).

My Grass knowledge isn't as strong as it could be, so maybe I'm doing something stupid. Definitely, check out that suggestion first as it'll be a lot less work than mine ;-)

Step 1 - Create a cost surface

First step is to create a cost surface. This only needs to be done once.

• create an editable layer, holes and all.
• add a field called 'unit', set it to 1.
• using polygon-to-raster on your "punched out" vector layer (the one which has the holes), using field 'unit'. You now have a layer "mask", where 1 is free space and 0 is building.

• use raster calculator to turn this into a cost surface. I'll set 'outdoors' to 1 and 'indoors' to 9999. This will make moving through buildings prohibitively difficult.

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

You can get more 'organic' results by adding a bit of noise to the cost surface (e.g. use random number from 1 to 3, rather than just 1 for outdoor pxiels.)

Step 2. Create cumulative cost rasters for each voronoi center

Now we can run (for one voronoi cell at a time) the GRASS algorithm r.cost.coordinates against our cost surface layer.

For the start coordinate, use the vornoi center. For the end coordinate, choose one of the corners of your area. I suggest using 'Knights Tour' as this gives smoother results.

The result shows lines of equal travel time from one voronoi center. Note how the bands wrap around the buildings.

Not sure how best to automate this. Maybe processing batch mode, or done in pyqgis.

Step 3. Merge down the rasters

This will probably need code. The algorithm would be

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

That approach should yield a raster where each cell is categorised by the voronoi center it's closest to, taking into account obstacles.

You could then use raster-to-polygon. You could then use the Generalise plugin to remove the "step" effect artifacts from the raster.

Apologies for vagueness on steps 2 and 3... I'm hoping someone chimes in with a more elegant solution :)

Answer 3

Note #1: I was not able to reproduce the proposed issue because the Difference tool worked well for me in several tests that I performed (maybe it was due to the simple geometry of the issue or because the tool has been improved since the question was asked 1 year ago).

However, I propose a workaround in PyQGIS to avoid the using of the Difference tool. Everything is based on the local intersection between two input layers (see figure below):

1. a polygon vector layer representing the Voronoi polygons;
2. a polygon vector layer representing the holes/constraints which need to be excluded from the analysis.

Note #2: Since I don't want to use the Difference tool, I'm not able to avoid the creation of "slivers" (see then), so I needed to run the v.clean tool to eliminate them. Furthermore, as @Chris W said,

[...] but the second resulting in slivers is not entirely unexpected - after all, that area would still be allocated to the same point even if the hole is present. You could use that method and then incorporate the slivers into their neighbors with a cleanup method.

After these necessary premises, I post my 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)

which leads to this result:

Just for clearness, this would be the result without the using of the v.clean tool:

The difference with the result by @LeaningCactus is that, by now, the geometries are not broken and they could be "cleaned" without errors.