This is one of the 100 recipes of the IPython Cookbook, the definitive guide to high-performance scientific computing and data science in Python.
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)
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
df = pd.read_csv('data/ratp.csv', sep='#', header=None)
df[df.columns[1:]].tail(3)
metro = df[(df[5] == 'metro')]
metro[metro.columns[1:]].tail(3)
str
attribute of the corresponding column.# We only extract the district from stations in Paris.
paris = metro[4].str.startswith('PARIS').values
# 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
lon = metro[1]
lat = metro[2]
box = (lat[paris].min(), lon[paris].min(),
lat[paris].max(), lon[paris].max())
m = smopy.Map(box, z=12)
Voronoi
object is created with the points coordinates, and contains several attributes we will use for display.vor = spatial.Voronoi(np.c_[lat, lon])
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)
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.regions, vertices = voronoi_finite_polygons_2d(vor)
cells = [m.to_pixels(vertices[region]) for region in regions]
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]]
show_mpl
method of the Map
instance.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).