# The topology of Paris as seen by my alarm clock

In this post, we will use a little geography to determine the topology of Paris as seen from the point of view of my alarm clock. Here's the context: I'd like to move flats in the Paris area and take advantage of the transportation means set up by my company. My company provides a number of bus routes travelling to our offices every morning, with each bus route consisting of a number of stops at given times.

Using some data that I collected on my company's internal website and coordinates from Google Maps, the goal of this post can be described as:

"which location within Paris would allow me to get out of the house the latest and still get to work on time?"

Let's see if how we can answer this question.

# The data¶

Firsts, we load the data for the bus departure times. We will encode the departure time as the minutes elapsed since midnight.

In [1]:
departure_times = list(map(lambda s: (int(s.split('\t')[0]) * 60 + int(s.split('\t')[1])),
"""7	30
7	35
7	36
7	40
7	45
7	15
7	27
7	25
7	29
7	27
7	45
7	45
7	40
7	30
7	35
7	33
7	30
7	27
7	20
7	13""".split('\n')))


We now define the location coordinates for each bus in (latitude, longitude) format.

In [2]:
location_coords = list(map(lambda s: (float(s.split(',')[0]), float(s.split(',')[1])),
"""48.878428, 2.282963
48.870967, 2.275274
48.863403, 2.268703
48.846252, 2.257876
48.837814, 2.256867
48.837877, 2.274259
48.836091, 2.277381
48.832199, 2.287949
48.826994, 2.300598
48.827143, 2.304493
48.819960, 2.324732
48.815806, 2.344237
48.819381, 2.359011
48.830428, 2.352359
48.831459, 2.398214
48.835386, 2.406239
48.846754, 2.411582
48.853348, 2.410273
48.863569, 2.408535
48.876956, 2.406550""".split('\n')))


# Visualizing the data¶

With the data points having been defined, we can now visualize tan approximate map of the bus departure times using the latitudes and longitudes of the bus stops as X,Y coordinates.

In [3]:
from pylab import *
%matplotlib inline

In [4]:
format_hour = lambda hour: "{0}h{1}".format(hour // 60, hour % 60)

for hour, (lat, lng) in zip(departure_times, location_coords):
text(lng, lat, format_hour(hour), bbox=dict(facecolor='red', alpha=0.5), horizontalalignment='center')
xlabel("longitude")
ylabel("latitude")
xlim(2.23, 2.45)
ylim(48.81, 48.89)

Out[4]:
(48.81, 48.89)

This being done, we can now start doing some computations!

# Computing the nearest distance to a bus stop¶

As we're going to do some computation on a map, let's define a mesh we can use to compute stuff:

In [5]:
lat_min = min([c[0] for c in location_coords])

In [6]:
lat_max = max([c[0] for c in location_coords])

In [7]:
long_min = min([c[1] for c in location_coords])

In [8]:
long_max = max([c[1] for c in location_coords])


Our bus stops are located within the following coordinates:

In [9]:
print(lat_min, lat_max, long_min, long_max)

48.815806 48.878428 2.256867 2.411582


However, to compute stuff, we'll use a grid that's a little bit larger than just the coordinates of the bus stops:

In [10]:
x = linspace(long_min - 0.05, long_max + 0.05, 250)
y = linspace(lat_min - .03, lat_max + 0.03, 200)

In [11]:
X, Y = meshgrid(x, y)


The next thing we will need is a function able to compute a distance between two points defined by their respective latitudes and longitudes. The math for this is a little complex but we're lucky: someone already wrote an implementation in Python, which we're going to use. Thanks John D. Cook: http://www.johndcook.com/blog/python_longitude_latitude/!

In [12]:
import math

def distance_on_unit_sphere(lat1, long1, lat2, long2):

# Convert latitude and longitude to

# phi = 90 - latitude

# theta = longitude

# Compute spherical distance from spherical coordinates.

# For two locations in spherical coordinates
# (1, theta, phi) and (1, theta, phi)
# cosine( arc length ) =
#    sin phi sin phi' cos(theta-theta') + cos phi cos phi'
# distance = rho * arc length

cos = (math.sin(phi1)*math.sin(phi2)*math.cos(theta1 - theta2) +
math.cos(phi1)*math.cos(phi2))
arc = math.acos( cos )

# Remember to multiply arc by the radius of the earth
# in your favorite set of units to get length.
return 6371 * arc # distance will be in km


Our first question is: how far is each point in our mesh from the nearest bus stop? It turns out this is quite easy to compute, using a little loop.

In [13]:
distance_to_bus = nan * zeros_like(X)

for i in range(x.size):
for j in range(y.size):
for k in range(len(location_coords)):
dist_to_loc = distance_on_unit_sphere(Y[j,i], X[j, i],
location_coords[k][0], location_coords[k][1])
if isnan(distance_to_bus[j, i]):
distance_to_bus[j, i] = dist_to_loc
else:
distance_to_bus[j, i] = min(distance_to_bus[j, i], dist_to_loc)


We can plot the computed distance to the nearest bus stops as an approximate map using latitude and longitude:

In [14]:
figure(figsize=(10, 6))
imshow(distance_to_bus, origin='lower', extent=(x.min(), x.max(), y.min(), y.max()), aspect='auto')
colorbar(shrink=0.5)
title('distance to bus stop (km)')
xlabel("longitude")
ylabel("latitude")

Out[14]:

This is a nice visualization, but we can do much better if we use the appropriate toolkit! The basemap package that works with matplotlib is exactly what we are looking for.

In [15]:
from mpl_toolkits.basemap import Basemap

In [16]:
from matplotlib.collections import LineCollection
import shapefile

figure(figsize=(10, 10))
ax = subplot(111)

# first, choose map coordinates and draw rivers
m = Basemap(projection='stere',
resolution='h',
llcrnrlon=x.min(),
llcrnrlat=y.min(),
urcrnrlon=x.max(),
urcrnrlat=y.max(),
lon_0=x.mean(),
lat_0=y.mean())
m.drawrivers() # the Seine

shapes = r.shapes()
records = r.records()
for record, shape in zip(records,shapes):
if record[6] == 75: # filtering Paris data
lons,lats = zip(*shape.points)
data = np.array(m(lons, lats)).T

if len(shape.parts) == 1:
segs = [data,]
else:
segs = []
for i in range(1,len(shape.parts)):
index = shape.parts[i-1]
index2 = shape.parts[i]
segs.append(data[index:index2])
segs.append(data[index2:])

lines = LineCollection(segs,antialiaseds=(1,))
lines.set_edgecolors('k')
lines.set_linewidth(0.1)

cb.set_label('km')
ax.set_title('distance to nearest bus stop')

Out[16]:

# Computing the time at which I need to leave to get to the bus on time¶

The distances to the nearest bus stop, which we have plotted in the previous image, are nice, but they don't take into account the time at which the bus leaves. To take this information into account in the plot, we can compute the time at which we need to leave a point on the mesh to reach the nearest bus on time.

We do this by using the following formula:

$$t_{\text{leave_appartment}}(x, y) = t_{\text{nearest_bus_leaves}} - \text{distance_to_bus} * \text{walking_speed}$$
In [17]:
walking_speed = 10. # minutes per kilometer

In [18]:
start_time = nan * zeros_like(X)

for i in range(x.size):
for j in range(y.size):
for k in range(len(location_coords)):
dist_to_loc = distance_on_unit_sphere(Y[j,i], X[j, i],
location_coords[k][0], location_coords[k][1])
time_to_leave = departure_times[k] - dist_to_loc * walking_speed
if isnan(start_time[j, i]):
start_time[j, i] = time_to_leave
else:
start_time[j, i] = max(start_time[j, i], time_to_leave)


We can plot the resulting departure times, just like for the distances:

In [19]:
labels = array([6.5, 6.75, 7, 7.25, 7.5, 7.75]) * 60
ticklabels = ['6h30', '6h45', '7h', '7h15', '7h30', '7h45']

figure(figsize=(10, 6))
imshow(start_time, origin='lower', extent=(x.min(), x.max(), y.min(), y.max()), aspect='auto',
vmin=labels.min(), vmax=labels.max())
cbar = colorbar(shrink=0.5)
cbar.set_ticks(labels)
cbar.set_ticklabels(ticklabels)
title('time to leave the house in the morning')
xlabel("longitude")
ylabel("latitude")

Out[19]: