This is one of the 100 recipes of the IPython Cookbook, the definitive guide to high-performance scientific computing and data science in Python.

14.5. Computing the Voronoi diagram of a set of points

You will need the Smopy module to display the OpenStreetMap map of Paris. You can install this package with pip install smopy.

Download the Metro dataset on the book's website and extract it in the current directory. (https://ipython-books.github.io)

  1. Let's import NumPy, Pandas, matplotlib, and SciPy, which implements a Voronoi diagram algorithm.
In [ ]:
import numpy as np
import pandas as pd
import scipy.spatial as spatial
import matplotlib.pyplot as plt
import matplotlib.path as path
import matplotlib as mpl
import smopy
%matplotlib inline
  1. Let's load the dataset with Pandas.
In [ ]:
df = pd.read_csv('data/ratp.csv', sep='#', header=None)
In [ ]:
df[df.columns[1:]].tail(3)
  1. The DataFrame contains the coordinates, name, city, district and type of stations from RATP (the society that manages public transportations in Paris). Let's select the metro stations only.
In [ ]:
metro = df[(df[5] == 'metro')]
In [ ]:
metro[metro.columns[1:]].tail(3)
  1. We are going to extract the district number of the stations that are located in Paris. With Pandas, we can use vectorized string operations using the str attribute of the corresponding column.
In [ ]:
# We only extract the district from stations in Paris.
paris = metro[4].str.startswith('PARIS').values
In [ ]:
# We create a vector of integers with the district number of
# the corresponding station, or 0 if the station is not
# in Paris.
districts = np.zeros(len(paris), dtype=np.int32)
districts[paris] = metro[4][paris].str.slice(6, 8).astype(np.int32)
districts[~paris] = 0
ndistricts = districts.max() + 1
  1. We also extract the coordinates of all metro stations.
In [ ]:
lon = metro[1]
lat = metro[2]
  1. Now, let's retrieve the map of Paris with OpenStreetMap. We specify the map's boundaries with the extreme latitude and longitude coordinates of all our metro stations. We use the lightweight smopy module for generating the map.
In [ ]:
box = (lat[paris].min(), lon[paris].min(), 
       lat[paris].max(), lon[paris].max())
m = smopy.Map(box, z=12)
  1. We now compute the Voronoi diagram of the stations using SciPy. A Voronoi object is created with the points coordinates, and contains several attributes we will use for display.
In [ ]:
vor = spatial.Voronoi(np.c_[lat, lon])
  1. We need to define a custom function to display the Voronoi diagram. Although SciPy comes with its own display function, it does not take infinite points into account. This function can be used as is every time you need to display a Voronoi diagram.
In [ ]:
def voronoi_finite_polygons_2d(vor, radius=None):
    """Reconstruct infinite Voronoi regions in a 
    2D diagram to finite regions.
    Source: http://stackoverflow.com/questions/20515554/colorize-voronoi-diagram
    """
    if vor.points.shape[1] != 2:
        raise ValueError("Requires 2D input")
    new_regions = []
    new_vertices = vor.vertices.tolist()
    center = vor.points.mean(axis=0)
    if radius is None:
        radius = vor.points.ptp().max()
    # Construct a map containing all ridges for a given point
    all_ridges = {}
    for (p1, p2), (v1, v2) in zip(vor.ridge_points, 
                                  vor.ridge_vertices):
        all_ridges.setdefault(p1, []).append((p2, v1, v2))
        all_ridges.setdefault(p2, []).append((p1, v1, v2))
    # Reconstruct infinite regions
    for p1, region in enumerate(vor.point_region):
        vertices = vor.regions[region]
        if all(v >= 0 for v in vertices):
            # finite region
            new_regions.append(vertices)
            continue
        # reconstruct a non-finite region
        ridges = all_ridges[p1]
        new_region = [v for v in vertices if v >= 0]
        for p2, v1, v2 in ridges:
            if v2 < 0:
                v1, v2 = v2, v1
            if v1 >= 0:
                # finite ridge: already in the region
                continue
            # Compute the missing endpoint of an infinite ridge
            t = vor.points[p2] - vor.points[p1] # tangent
            t /= np.linalg.norm(t)
            n = np.array([-t[1], t[0]])  # normal
            midpoint = vor.points[[p1, p2]].mean(axis=0)
            direction = np.sign(np.dot(midpoint - center, n)) * n
            far_point = vor.vertices[v2] + direction * radius
            new_region.append(len(new_vertices))
            new_vertices.append(far_point.tolist())
        # Sort region counterclockwise.
        vs = np.asarray([new_vertices[v] for v in new_region])
        c = vs.mean(axis=0)
        angles = np.arctan2(vs[:,1] - c[1], vs[:,0] - c[0])
        new_region = np.array(new_region)[np.argsort(angles)]
        new_regions.append(new_region.tolist())
    return new_regions, np.asarray(new_vertices)
  1. The voronoi_finite_polygons_2d function returns a list of regions and a list of vertices. Every region is a list of vertex indices. The coordinates of all vertices are stored in vertices. From these structures, we can create a list of cells. Every cell represents a polygon as an array of vertex coordinates. We also use the to_pixels method of the smopy.Map instance which converts latitude and longitude geographical coordinates to pixels in the image.
In [ ]:
regions, vertices = voronoi_finite_polygons_2d(vor)
cells = [m.to_pixels(vertices[region]) for region in regions]
  1. Now we compute the color of every polygon.
In [ ]:
cmap = plt.cm.Set3
# We generate colors for districts using a color map.
colors_districts = cmap(np.linspace(0., 1., ndistricts))[:,:3]
# The color of every polygon, grey by default.
colors = .25 * np.ones((len(districts), 3))
# We give each polygon in Paris the color of its district.
colors[paris] = colors_districts[districts[paris]]
  1. Finally, we can display the map with the Voronoi diagram, using the show_mpl method of the Map instance.
In [ ]:
ax = m.show_mpl();
ax.add_collection(mpl.collections.PolyCollection(cells,
                  facecolors=colors, edgecolors='k',
                  alpha=.35,));

You'll find all the explanations, figures, references, and much more in the book (to be released later this summer).

IPython Cookbook, by Cyrille Rossant, Packt Publishing, 2014 (500 pages).