From Voronoi tessellation to Shapely polygons

Question

from a set of points I built the Voronoi tessellation using scipy:

from scipy.spatial import Voronoi
vor = Voronoi(points)

Now I would like to build a Polygon in Shapely from the regions the Voronoi algorithm created. The problem is that the Polygon class requires a list of counter-clockwise vertices. Although I know how to order these vertices, I can't solve the problem because often this is my result:

(overlapping polygon). This is the code (ONE RANDOM EXAMPLE):

def order_vertices(l):
    mlat = sum(x[0] for x in l) / len(l)
    mlng = sum(x[1] for x in l) / len(l)

    # https://mcmap.net/q/659209/-how-can-i-sort-a-coordinate-list-for-a-rectangle-counterclockwise
    def algo(x):
        return (math.atan2(x[0] - mlat, x[1] - mlng) + 2 * math.pi) % 2*math.pi

    l.sort(key=algo)
    return l

a = np.asarray(order_vertices([(9.258054711746084, 45.486245994138976),
 (9.239284166975443, 45.46805963143515),
 (9.271640747003861, 45.48987234571072),
 (9.25828782103321, 45.44377372506324),
 (9.253993275176263, 45.44484395950612),
 (9.250114174032936, 45.48417979682819)]))
plt.plot(a[:,0], a[:,1])

How can I solve this problem?

Whitsun

That doesn't look like a counterclockwise ordering of vertices. Are you sure you implemented the ordering algorithm properly? – Andino

@Andino I think I implemented it correctly... I updated the question with an example – Whitsun

Is that the code that generated your graph? Or is it a separate example? – Andino

@Andino that code generated the graph :) – Whitsun

Answer 1

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 np
from scipy.spatial import Voronoi, voronoi_plot_2d
import shapely.geometry
import shapely.ops

points = np.random.random((10, 2))
vor = Voronoi(points)
voronoi_plot_2d(vor)

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

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

Hoashis

How can you save the shapefile with the polygons? – Oconnor

@ValerioD.Ciotti Use shapely.geometry.mapping on the polygon(s) to cast them as geojson objects, and then use a library like fiona to write the geojson to a new shapefile. From there if you're stuck ask questions on http://gis.stackexchange.com/. – Hoashis

@ValerioD.Ciotti Or you can make geopandas.GeoSeries from the polygons and save it directly into shapefile. – Hirohito

You can also use the fiona library directly, which GeoPandas already depends on. – Jellaba

See this answer for a way to handle the infinite Voronoi regions. – Exodus

Answer 2

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.

Maimaia

Answer 3

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

import numpy as np
from 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

Baileybailie

Answer 4

The function you implemented (order_vertices() ), cannot work in your case, because it simply takes a already ordered sequence of coordinates, which builds a rectangle, and inverts the direction of the polygon (and maybe works only with rectangles...). But you have a not ordered sequence of coordinates

Generally speaking, you can not build a polygon from a arbitrary sequence of not ordered vertices because there is no a unique solution for concave polygons, as showed in this example: https://mcmap.net/q/659211/-given-a-vector-of-points-possibly-out-of-order-find-polygon-not-convex-hull

However, if you are sure that your polygons are always convex, you can build a convex hull with this code: https://mcmap.net/q/659212/-convex-hull-and-scipy (tested right now, it worked for me)

Probably you can build the convex hull with scipy as well, but I dint test it: scipy.spatial.ConvexHull

Revere