From Voronoi tessellation to Shapely polygons

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.

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