Trajectories in the French department landscape

First, let's load a list of the French departments.

In [1]:
src = u"""01	Ain	Bourg-en-Bresse	Rhône-Alpes
02	Aisne	Laon	Picardie
03	Allier	Moulins	Auvergne
04	Alpes de Hautes-Provence	Digne	Provence-Alpes-Côte d'Azur
05	Hautes-Alpes	Gap	Provence-Alpes-Côte d'Azur
06	Alpes-Maritimes	Nice	Provence-Alpes-Côte d'Azur
07	Ardèche	Privas	Rhône-Alpes
08	Ardennes	Charleville-Mézières	Champagne-Ardenne
09	Ariège	Foix	Midi-Pyrénées
10	Aube	Troyes	Champagne-Ardenne
11	Aude	Carcassonne	Languedoc-Roussillon
12	Aveyron	Rodez	Midi-Pyrénées
13	Bouches-du-Rhône	Marseille	Provence-Alpes-Côte d'Azur
14	Calvados	Caen	Basse-Normandie
15	Cantal	Aurillac	Auvergne
16	Charente	Angoulême	Poitou-Charentes
17	Charente-Maritime	La Rochelle	Poitou-Charentes
18	Cher	Bourges	Centre
19	Corrèze	Tulle	Limousin
2A	Corse-du-Sud	Ajaccio	Corse
2B	Haute-Corse	Bastia	Corse
21	Côte-d'Or	Dijon	Bourgogne
22	Côtes d'Armor	Saint-Brieuc	Bretagne
23	Creuse	Guéret	Limousin
24	Dordogne	Périgueux	Aquitaine
25	Doubs	Besançon	Franche-Comté
26	Drôme	Valence	Rhône-Alpes
27	Eure	Évreux	Haute-Normandie
28	Eure-et-Loir	Chartres	Centre
29	Finistère	Quimper	Bretagne
30	Gard	Nîmes	Languedoc-Roussillon
31	Haute-Garonne	Toulouse	Midi-Pyrénées
32	Gers	Auch	Midi-Pyrénées
33	Gironde	Bordeaux	Aquitaine
34	Hérault	Montpellier	Languedoc-Roussillon
35	Ille-et-Vilaine	Rennes	Bretagne
36	Indre	Châteauroux	Centre
37	Indre-et-Loire	Tours	Centre
38	Isère	Grenoble	Rhône-Alpes
39	Jura	Lons-le-Saunier	Franche-Comté
40	Landes	Mont-de-Marsan	Aquitaine
41	Loir-et-Cher	Blois	Centre
42	Loire	Saint-Étienne	Rhône-Alpes
43	Haute-Loire	Le Puy-en-Velay	Auvergne
44	Loire-Atlantique	Nantes	Pays de la Loire
45	Loiret	Orléans	Centre
46	Lot	Cahors	Midi-Pyrénées
47	Lot-et-Garonne	Agen	Aquitaine
48	Lozère	Mende	Languedoc-Roussillon
49	Maine-et-Loire	Angers	Pays de la Loire
50	Manche	Saint-Lô	Basse-Normandie
51	Marne	Châlons-en-Champagne	Champagne-Ardenne
52	Haute-Marne	Chaumont	Champagne-Ardenne
53	Mayenne	Laval	Pays de la Loire
54	Meurthe-et-Moselle	Nancy	Lorraine
55	Meuse	Bar-le-Duc	Lorraine
56	Morbihan	Vannes	Bretagne
57	Moselle	Metz	Lorraine
58	Nièvre	Nevers	Bourgogne
59	Nord	Lille	Nord-Pas-de-Calais
60	Oise	Beauvais	Picardie
61	Orne	Alençon	Basse-Normandie
62	Pas-de-Calais	Arras	Nord-Pas-de-Calais
63	Puy-de-Dôme	Clermont-Ferrand	Auvergne
64	Pyrénées-Atlantiques	Pau	Aquitaine
65	Hautes-Pyrénées	Tarbes	Midi-Pyrénées
66	Pyrénées-Orientales	Perpignan	Languedoc-Roussillon
67	Bas-Rhin	Strasbourg	Alsace
68	Haut-Rhin	Colmar	Alsace
69	Rhône	Lyon	Rhône-Alpes
70	Haute-Saône	Vesoul	Franche-Comté
71	Saône-et-Loire	Mâcon	Bourgogne
72	Sarthe	Le Mans	Pays de la Loire
73	Savoie	Chambéry	Rhône-Alpes
74	Haute-Savoie	Annecy	Rhône-Alpes
75	Paris	Paris	Ile-de-France
76	Seine-Maritime	Rouen	Haute-Normandie
77	Seine-et-Marne	Melun	Ile-de-France
78	Yvelines	Versailles	Ile-de-France
79	Deux-Sèvres	Niort	Poitou-Charentes
80	Somme	Amiens	Picardie
81	Tarn	Albi	Midi-Pyrénées
82	Tarn-et-Garonne	Montauban	Midi-Pyrénées
83	Var	Toulon	Provence-Alpes-Côte d'Azur
84	Vaucluse	Avignon	Provence-Alpes-Côte d'Azur
85	Vendée	La Roche-sur-Yon	Pays de la Loire
86	Vienne	Poitiers	Poitou-Charentes
87	Haute-Vienne	Limoges	Limousin
88	Vosges	Épinal	Lorraine
89	Yonne	Auxerre	Bourgogne
90	Territoire-de-Belfort	Belfort	Franche-Comté
91	Essonne	Évry	Ile-de-France
92	Hauts-de-Seine	Nanterre	Ile-de-France
93	Seine-Saint-Denis	Bobigny	Ile-de-France
94	Val-de-Marne	Créteil	Ile-de-France
95	Val-d'Oise	Pontoise	Ile-de-France"""
In [2]:
depts = dict([line.split('\t')[:2] for line in src.split('\n')])
In [3]:
for d in ['01', '06', '23', '56', '95']:
    print depts[d]
Ain
Alpes-Maritimes
Creuse
Morbihan
Val-d'Oise

Now, on to some parsing of phone numbers.

In [4]:
phone_number = '0676237181'
In [7]:
def number_to_depts(number):
    output = []
    for i in range(5):
        try:
            output.append(depts[number[2*i:2*(i+1)]])
        except KeyError:
            output.append('**')
    return output
In [8]:
number_to_depts(phone_number)
Out[8]:
[u'Alpes-Maritimes',
 u'Seine-Maritime',
 u'Creuse',
 u'Sa\xf4ne-et-Loire',
 u'Tarn']

Finally, we can represent this on a map, using coordinates of the prefectures of each department. The first step for this is to generate all coordinates for the prefecture towns of France.

In [9]:
from geopy import geocoders
In [10]:
g = geocoders.GoogleV3()
adr = 'Bourg-en-Bresse, France'
place, (lat, lng) = g.geocode(adr)
In [11]:
place, lat, lng
Out[11]:
(u'Bourg-en-Bresse, France', 46.20279, 5.219246)

As this seems to work, I'm gonna generate a table of coordinates for all prefectures:

In [12]:
prefectures = dict([(line.split('\t')[0], ', '.join((line.split('\t')[2], line.split('\t')[1], 'France'))) for line in src.split('\n')])
In [13]:
for d in ['01', '06', '23', '56', '95']:
    print prefectures[d]
Bourg-en-Bresse, Ain, France
Nice, Alpes-Maritimes, France
Guéret, Creuse, France
Vannes, Morbihan, France
Pontoise, Val-d'Oise, France
In [14]:
coords = {}
In [16]:
import time
for item in prefectures.keys():
    if not item in coords:
        adr = prefectures[item]
        place, (lat, lng) = g.geocode(adr, exactly_one=False)[0]
        coords[item] = [place, (lat, lng)]
        time.sleep(0.2)
In [17]:
for d in ['01', '06', '23', '56', '95']:
    print coords[d]
[u'Bourg-en-Bresse, France', (46.20279, 5.219246)]
[u'Nice, France', (43.696036, 7.265592)]
[u'23000 Gu\xe9ret, France', (46.16959900000001, 1.871452)]
[u'Vannes, France', (47.658236, -2.760847)]
[u'Pontoise, France', (49.050966, 2.100645)]

Actually, we're going to save the data about geolocation to disk, so that we don't need to query Google Maps every time we want to execute this notebook. We'll do this with in a plain text file.

In [35]:
import os

filename = "files/pref_geolocalization.txt" 

if not os.path.exists(filename):
    with open(filename, 'w') as f:
        for item in [', '.join((item[0], str(item[1][1][0]), str(item[1][1][1]))) for item in coords.items()]:
            f.write(item +'\n')

We can now plot the trajectories on the map. First, let's setup the basic map we're going to use to plot the trajectories.

In [36]:
from mpl_toolkits.basemap import Basemap
In [37]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib
WARNING: pylab import has clobbered these variables: ['place', 'f']
`%pylab --no-import-all` prevents importing * from pylab and numpy
In [38]:
def plot_French_map():
    # setup Lambert Conformal basemap.
    figure(figsize=(10, 10))
    m = Basemap(width=1000000, height=1000000, projection='lcc',
                resolution='l',
                lat_0=47.,lon_0=2.)
    # draw coastlines.
    m.drawcoastlines()
    # draw a boundary around the map, fill the background.
    # this background will end up being the ocean color, since
    # the continents will be drawn on top.
    m.drawmapboundary(fill_color='aqua')
    # fill continents, set lake color same as ocean color.
    m.fillcontinents(color='coral',lake_color='aqua')
    m.drawcountries()
    return m
plot_French_map()
Out[38]:
<mpl_toolkits.basemap.Basemap at 0x4a09170>

Looks good enough for our purpose, now on to the real plotting with arrows.

In [41]:
def plot_number_on_map(number):
    double_digits = []
    for i in range(5):
        double_digits.append(number[2*i:2*(i+1)])
    trajectory = []
    for i in range(5):
        trajectory.append(coords[double_digits[i]])
    # plot background map of France
    m = plot_French_map()
    # plot arrows between two successive cities
    for i in range(4):
        lat0, lon0 = trajectory[i][1]
        x0, y0 = m(lon0, lat0)
        lat1, lon1 = trajectory[i + 1][1]
        x1, y1 = m(lon1, lat1)
        #m.plot([x0, x1], [y0, y1], color='g', linewidth=5)
        arrow(x0, y0, x1-x0, y1-y0) #, linewidth=3, head_width=20000, head_length=20000)
    # plot dots for the prefectures and text labels
    for i in range(5):
        latitudes = trajectory[i][1][0]
        longitudes = trajectory[i][1][1]
        x, y = m(longitudes, latitudes)
        m.plot(x, y, 'bo', markersize=10)
        text(x - 100000, y + 20000, prefectures[double_digits[i]])
        title("phone number: %s" % number)
plot_number_on_map(phone_number)
 

We can now test some other phone numbers...

In [42]:
plot_number_on_map("0641677342")

Comments