#!/usr/bin/env python # coding: utf-8 # > This is one of the 100 recipes of the [IPython Cookbook](http://ipython-books.github.io/), 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 get_ipython().run_line_magic('matplotlib', 'inline') # 2. 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) # 3. 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) # 4. 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 # 5. We also extract the coordinates of all metro stations. # In[ ]: lon = metro[1] lat = metro[2] # 6. 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) # 7. 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]) # 8. 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) # 9. 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] # 10. 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]] # 11. 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](http://ipython-books.github.io/), by [Cyrille Rossant](http://cyrille.rossant.net), Packt Publishing, 2014 (500 pages).