7

Using Qgis 3.2 or python's geovoronoi, the Voronoi algorithm does not work on 2 points. I expected that the bounding box between the 2 points would have been cut in half but got errors instead. Qgis stating that at least 3 points are needed.

I'm interested in the reason why at least 3 points are needed for Voronoi?

1
  • 1
    This is all I could find. I am not sure whether the issue stems from the calculation method or the expected length of the tuples inside the function.
    – ChloeG
    Mar 7, 2022 at 13:02

1 Answer 1

4

I can't fully answer this, but I can confirm that it if you bypass the check for two points, it will almost always work, generating two Voronoi Polygons.

The source code for the Voronoi Polygons tool is in a file called VoronoiPolygons.py.

  • On Windows, this is in qgis\python\plugins\processing\algs\qgis.
  • On macOS it's in QGIS.app/Contents/Resources/python/plugins/processing/algs/qgis. (You can look inside a .app in the Finder by right-clicking on it and selecting "Show Package Contents".

As of QGIS 3.22.0, lines 151 to 160 of this file were:

if len(pts) < 3:
    raise QgsProcessingException(
        self.tr('Input file should contain at least 3 points. Choose '
                'another file and try again.'))
# Eliminate duplicate points
uniqueSet = set(item for item in pts)
ids = [pts.index(item) for item in uniqueSet]
sl = voronoi.SiteList([voronoi.Site(i[0], i[1], sitenum=j)
                       for (j, i) in enumerate(uniqueSet)])
voronoi.voronoi(sl, c)

Note that it checks the number of points, then removes duplicates. This means there are 3 ways you can see what happens running the algorithm on 2 points:

  1. Create a point layer with 3 points, 2 of them with duplicated coordinates. The uniqueSet line will reduce this to 2 points.
  2. Close QGIS, back up your existing VoronoiPolygons.py, edit the if len(pts) < 3: to if len(pts) < 2, and restart QGIS. When you finish, restore the original VoronoiPolygons.py. Seriously. Don't leave edited files within your QGIS installation. You could cause a lot of damage that way.
  3. Create a script or a whole plugin that duplicates VoronoiPolygons.py, including the voronoi.py file it imports. Then experiment on that as in step 2, but without messing with any of your existing files.

What I've seen from any of these methods is results like the ones below. The top row shows 2 points in various orientations. The bottom row shows 3 points for comparison.

  • Most of the time, it generates the Voronoi Polygons we expect.
  • If the points are all in line horizontally or vertically, it outputs the error message The extent of the input points is not a polygon (all the points are on a vertical or horizontal line) - cannot make a Voronoi diagram! and exits without generating anything.
  • In the specific case I've highlighted in red, where the two points are in line southwest to northeast, it outputs an error message of the type Feature could not be written to Voronoi_polygons_db6ff839_d2b5_45fe_a6ab_5ea2c5613d41: Could not add feature with geometry type LineString to layer of type Polygon and generates the southwest triangle but not the northeast one.

Points in pairs and triples in different orientations, and the Voronoi Polygons generated around them. One is highlighted in red.

This happened at a variety of origin positions and scales, but if someone wants to duplicate it and has trouble, below is a file where it happened to me.

Longitude,Latitude,Name
10,10,Northeast
-10,-10,Southwest

I've investigated some of this, and it seems that the output from one of the algorithm steps only defines 2 points when it should be at least 4 (a triangle with the start and end point repeated). This is why it gets treated as a LineString, but can't by added to a Polygon layer.

Were the developers of the Voronoi polygons tool aware of this error and prevented it running on 2 points to avoid it happening? Possibly.

I hope this answers some of your questions. If you experiment with running the algorithm on 2 points, please let us know if you get the answers you expect.

Edited 11th March 2022

I have run many test cases through VoronoiPolygons.py. The VoronoiPolygons.clip_voronoi() method doesn't allow for a few unusual cases. These are only likely to appear with "toy problems" created to test it, and are highly improbable with "real world" inputs. I will be writing up an issue report for these, but it's relevant to the question, so I will document it here first. As also noted in the comments, none of these problems happen with QgsGeometry.voronoiDiagram().

As of QGIS 3.22.0, here are the problems:

Voronoi diagrams from 3 colinear points with corners missing in some cases.

  1. If a single polygon should have points in diagonally opposite corners, it will be created missing both of them, as in the diagrams above.

  2. If a polygon should be a triangle going to the very northeast, southeast and southwest corners of the bounding box, it will miss the northeast corner and therefore be treated as a LineString rather than a Polygon. This can only happen with a square bounding box and 2 points, one in the northeast and the other in the southwest corner. The existing check in VoronoiPolygons.processAlgorithm() requiring 3 or more points therefore prevents this from happening in most cases except in issue 3 below.

  3. VoronoiPolygons.processAlgorithm() counts the number of points before removing duplicates, which is the only way to cause issue 2.

3
  • Thanks. I will see if I need to change the .py files in the future. In my case I needed to loop through different regions and create constraint Voronoi areas per region. I decided to go the python route and used geovoronoi. Where there are only 2 points I created a long line perpendicular to the line between a b, where a is the first point and b the center point of the original 2 points. This line was then used to cut the region. see: stackoverflow.com/questions/57065080/…
    – user19349
    Mar 10, 2022 at 7:44
  • Just how do you access geovoronoi? I found the PyQGIS QgsGeometry.voronoiDiagram(), which is definitely a different implementation from the one in the toolkit, and doesn't throw up any of the weirdness described above.
    – Cowirrie
    Mar 11, 2022 at 11:13
  • 1
    Sorry if there were any confusion but I ended up using a standalone python script, not Qgis or PyQGIS. The geovoronoi library can be found at pypi.org/project/geovoronoi. So I used that in conjunction with geopandas and shapely. I won't mind sharing the script but if the QgsGeometry.voronoiDiagram() works, I think that is the way to go then. I just thought that if it does not work in Qgis it won't work in PyQGIS.
    – user19349
    Mar 11, 2022 at 13:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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