From Voronoi tessellation to Shapely polygons From Voronoi tessellation to Shapely polygons python python

From Voronoi tessellation to Shapely polygons


If you're just after a collection of polygons you don't need to pre-order the point to build them.

The scipy.spatial.Voronoi object has a ridge_vertices attribute containing indices of vertices forming the lines of the Voronoi ridge. If the index is -1 then the ridge goes to infinity.

First start with some random points to build the Voronoi object.

import numpy as npfrom scipy.spatial import Voronoi, voronoi_plot_2dimport shapely.geometryimport shapely.opspoints = np.random.random((10, 2))vor = Voronoi(points)voronoi_plot_2d(vor)

Voronoi plot from scipy.spatial

You can use this to build a collection of Shapely LineString objects.

lines = [    shapely.geometry.LineString(vor.vertices[line])    for line in vor.ridge_vertices    if -1 not in line]

The shapely.ops module has a polygonize that returns a generator for Shapely Polygon objects.

for poly in shapely.ops.polygonize(lines):    #do something with each polygon

Polygons from Voronoi with some sample points

Or if you wanted a single polygon formed from the region enclosed by the Voronoi tesselation you can use the Shapely unary_union method:

shapely.ops.unary_union(list(shapely.ops.polygonize(lines)))

Merged Voronoi tesselation polygon


As others have said, it is because you have to rebuild the polygons from the resulting points correctly based on indexes. Although you have the solution, I thought I should mention there is also another pypi supported tesselation package called Pytess (Disclaimer: I am the package maintainer) where the voronoi function returns the voronoi polygons fully built for you.


The library is able to generate ordered list of coordinates, you just need to make use of the index lists provided:

import numpy as npfrom scipy.spatial import Voronoi...ids = np.array(my_points_list)vor = Voronoi(points)polygons = {}for id, region_index in enumerate(vor.point_region):    points = []    for vertex_index in vor.regions[region_index]:        if vertex_index != -1:  # the library uses this for infinity            points.append(list(vor.vertices[vertex_index]))    points.append(points[0])    polygons[id]=points

each polygon in the polygons dictionary can be exported to geojson or brought into shapely and I was able to render them properly in QGIS