11

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

enter image description here

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.

  • This would be similar to but the opposite of (I think) using an environment mask in ArcGIS. That would let you constrain the created polygons to within a particular boundary. However I'm not aware of any tool that would make use of complex boundaries/holes (though maybe in ArcGIS the mask could be that complex - I haven't tested it and I might give it a try later if I have time). – Chris W Apr 27 '15 at 21:33
  • I've tested ArcGIS theory, and it won't work. Per the linked question, you can constrain results to an outer shape. However, a hole cut in the shape is ignored by the resulting polys. Furthermore, if that hole also has some of the points in it, the tool errors out and fails to run. I can't explain your first issue with difference, 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. – Chris W May 1 '15 at 23:43
  • 2
    You could potentially solve this by going to raster. With a raster mask, Euclidean distance would go out from your points until it hit either cells coming out from another point or your mask raster (your boundary slam description). Then you do some zonal cleanup and vectorize the result to get polygons. – Chris W May 1 '15 at 23:49
  • 1
    I would make sure the voronoi Geometry is valid by running v.clean then check geometry. Finally, execute Difference to create the holes. – klewis Dec 7 '16 at 21:33
  • What is Voronoi about these holes? Aren't you wanting to punch holes cleanly? Why wouldn't any polygon layer do? – mdsumner Dec 7 '16 at 23:00
3
+50

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.

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.

enter image description here

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

enter image description here

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.

enter image description here

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

  • 1
    Thanks Steven, I have a working GRASS raster work-around but I was hoping for a more elegant solution as mentioned in the bounty description. – underdark Dec 11 '16 at 20:54
0

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.

enter image description here

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:

enter image description here

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

enter image description here

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

  • Make to holes longer, e.g. cutting through the whole map, like a river, and you will see the issue. Adding slivers to neighbors creates polygons that look very different than a proper constrained Voronoi diagram. I tried that. – underdark Dec 11 '16 at 14:48
  • Sorry, I don't understand: did you find any error in the results? I only tested the code for that cases in which the polygons were similar to that ones proposed in the question. – mgri Dec 11 '16 at 15:07
  • Can't test the code right now unfortunately, but could you show the results obtained with the change in holes sketched in i.stack.imgur.com/Jpfra.png? – underdark Dec 11 '16 at 15:22
  • If I extend the constraint up to the feature on the right, I obtain this. Instead, if I directly move the constraint, I obtain this. – mgri Dec 11 '16 at 15:33
  • The small triangle that the red arrow in my drawing points to is the issue. It shouldn't be there but it's also in your results. Seems like this approach solves problem #1 of the question but leaves #2 unsolved. – underdark Dec 11 '16 at 18:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.