Tech - carte

Links: notebook, html, PDF, python, slides, GitHub

Faire une carte, c’est toujours compliqué. C’est simple jusqu’à ce qu’on s’aperçoive qu’on doit récupérer la description des zones administratives d’un pays, fournies parfois dans des coordonnées autres que longitude et latitude. Quelques modules utiles :

  • cartopy : surcouche de matplotlib pour faire des dessins avec des coordonnées géographiques

  • bokeh : pour tracer des cartes interactives

  • pyproj : conversion entre systèmes de coordonnées

  • shapely : manipuler des polygones géographiques (union, intersection, …)

  • pyshp : lire ou écrire des polygones géographiques

  • geopandas : manipulation de dataframe avec des coordonnées géographiques

Quelques notebooks intéressants : * Tracer une carte en Python avec bokeh * Tracer une carte en Python * Données carroyées et OpenStreetMap * Carte de France avec les départements * Carte de France avec les départements (2)

from jyquickhelper import add_notebook_menu
add_notebook_menu()
%matplotlib inline

Exposé

On télécharge des données hospitalières par départements.

Données COVID

# https://www.data.gouv.fr/fr/datasets/donnees-hospitalieres-relatives-a-lepidemie-de-covid-19/
from pandas import read_csv
url = "https://www.data.gouv.fr/fr/datasets/r/63352e38-d353-4b54-bfd1-f1b3ee1cabd7"
covid = read_csv(url, sep=";")
covid.tail()
dep sexe jour hosp rea rad dc
76048 974 1 2020-11-23 47 8 295 20
76049 974 2 2020-11-23 40 7 298 24
76050 976 0 2020-11-23 11 4 505 35
76051 976 1 2020-11-23 7 2 242 21
76052 976 2 2020-11-23 4 2 257 14
last_day = covid.loc[covid.index[-1], "jour"]
last_day
'2020-11-23'
last_data = covid[covid.jour == last_day].groupby("dep").sum()
last_data.shape
(101, 5)
last_data.describe()
sexe hosp rea rad dc
count 101.0 101.000000 101.000000 101.000000 101.000000
mean 3.0 619.396040 87.584158 3013.990099 670.326733
std 0.0 642.394039 108.421333 3876.642059 847.575190
min 3.0 22.000000 2.000000 214.000000 30.000000
25% 3.0 196.000000 26.000000 759.000000 168.000000
50% 3.0 419.000000 48.000000 1480.000000 346.000000
75% 3.0 751.000000 116.000000 3565.000000 738.000000
max 3.0 3337.000000 544.000000 19416.000000 4803.000000
last_data.head()
sexe hosp rea rad dc
dep
01 3 904 52 2148 515
02 3 430 54 3173 749
03 3 461 52 1341 354
04 3 262 15 695 103
05 3 368 28 835 168
last_data.tail()
sexe hosp rea rad dc
dep
971 3 192 14 1393 326
972 3 112 34 616 80
973 3 26 2 3698 132
974 3 174 30 1186 88
976 3 22 8 1004 70

Données départements

On récupère ensuite la définition géographique des départements.

import geopandas
# dernier lien de la page (format shapefiles)
url = "https://www.data.gouv.fr/en/datasets/r/ed02b655-4307-4db4-b1ca-7939145dc20f"
geo = geopandas.read_file(url)
geo.tail()
code_depart departement code_region region code_ancien ancienne_re geometry
96 87 Haute-Vienne 75.0 Nouvelle-Aquitaine 74.0 Limousin POLYGON ((0.92278 45.94610, 0.92856 45.94808, ...
97 63 Puy-de-Dôme 84.0 Auvergne-Rhône-Alpes 83.0 Auvergne POLYGON ((2.45492 45.76131, 2.45169 45.76137, ...
98 14 Calvados 28.0 Normandie 25.0 Basse-Normandie MULTIPOLYGON (((-0.42990 48.86405, -0.42881 48...
99 07 Ardèche 84.0 Auvergne-Rhône-Alpes 82.0 Rhône-Alpes POLYGON ((4.30746 44.98596, 4.30412 44.98837, ...
100 32 Gers 76.0 Occitanie 73.0 Midi-Pyrénées POLYGON ((-0.24307 43.66404, -0.24328 43.66476...

Il faudrait aussi fusionner avec la population de chaque département. Ce sera pour une autre fois.

Carte

import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 1, figsize=(5, 4))
geo.plot(ax=ax, color='white', edgecolor='black');
../_images/2020_carte_16_0.png

On enlève tous les départements à trois chiffres.

codes = [_ for _ in set(geo.code_depart) if len(_) < 3]
metropole = geo[geo.code_depart.isin(codes)]
metropole.tail()
code_depart departement code_region region code_ancien ancienne_re geometry
96 87 Haute-Vienne 75.0 Nouvelle-Aquitaine 74.0 Limousin POLYGON ((0.92278 45.94610, 0.92856 45.94808, ...
97 63 Puy-de-Dôme 84.0 Auvergne-Rhône-Alpes 83.0 Auvergne POLYGON ((2.45492 45.76131, 2.45169 45.76137, ...
98 14 Calvados 28.0 Normandie 25.0 Basse-Normandie MULTIPOLYGON (((-0.42990 48.86405, -0.42881 48...
99 07 Ardèche 84.0 Auvergne-Rhône-Alpes 82.0 Rhône-Alpes POLYGON ((4.30746 44.98596, 4.30412 44.98837, ...
100 32 Gers 76.0 Occitanie 73.0 Midi-Pyrénées POLYGON ((-0.24307 43.66404, -0.24328 43.66476...
fig, ax = plt.subplots(1, 1, figsize=(5, 4))
metropole.plot(ax=ax, color='white', edgecolor='black')
ax.set_title("%s départements" % metropole.shape[0]);
../_images/2020_carte_19_0.png

Carte COVID

merged = last_data.reset_index(drop=False).merge(metropole, left_on="dep", right_on="code_depart")
merged.shape
(96, 13)
merged.tail()
dep sexe hosp rea rad dc code_depart departement code_region region code_ancien ancienne_re geometry
91 91 3 1311 178 8567 1741 91 Essonne 11.0 Île-de-France 11.0 Île-de-France POLYGON ((2.04261 48.62638, 2.04045 48.62754, ...
92 92 3 1683 308 14598 2830 92 Hauts-de-Seine 11.0 Île-de-France 11.0 Île-de-France POLYGON ((2.22040 48.92061, 2.23114 48.92773, ...
93 93 3 1691 200 13101 2681 93 Seine-Saint-Denis 11.0 Île-de-France 11.0 Île-de-France POLYGON ((2.45949 48.95505, 2.46072 48.95840, ...
94 94 3 1807 231 13667 3221 94 Val-de-Marne 11.0 Île-de-France 11.0 Île-de-France POLYGON ((2.32906 48.81378, 2.33190 48.81701, ...
95 95 3 1047 106 8701 2052 95 Val-d'Oise 11.0 Île-de-France 11.0 Île-de-France POLYGON ((2.43548 49.13411, 2.44084 49.13419, ...
fig, ax = plt.subplots(1, 1, figsize=(5, 4))
merged.hist('rea', bins=20, ax=ax)
ax.set_title("Distribution rea");
../_images/2020_carte_23_0.png

Les régions les plus peuplées ont sans doute la plus grande capacité hospitalière. Il faudrait diviser par cette capacité pour avoir une carte qui ait un peu plus de sens. Comme l’idée est ici de simplement tracer la carte, on ne calculera pas de ratio.

merged.sort_values('rea').tail()
dep sexe hosp rea rad dc code_depart departement code_region region code_ancien ancienne_re geometry
92 92 3 1683 308 14598 2830 92 Hauts-de-Seine 11.0 Île-de-France 11.0 Île-de-France POLYGON ((2.22040 48.92061, 2.23114 48.92773, ...
59 59 3 3158 507 12412 2892 59 Nord 32.0 Hauts-de-France 31.0 Nord-Pas-de-Calais MULTIPOLYGON (((3.04058 50.15986, 3.04806 50.1...
69 69 3 3337 508 13746 2771 69 Rhône 84.0 Auvergne-Rhône-Alpes 82.0 Rhône-Alpes POLYGON ((4.43893 46.16789, 4.43463 46.17034, ...
12 13 3 2627 519 17760 2790 13 Bouches-du-Rhône 93.0 Provence-Alpes-Côte d'Azur 93.0 Provence-Alpes-Côte d'Azur MULTIPOLYGON (((5.38534 43.18909, 5.38440 43.1...
75 75 3 2365 544 19416 4803 75 Paris 11.0 Île-de-France 11.0 Île-de-France POLYGON ((2.27995 48.87857, 2.28101 48.88297, ...
geomerged = geopandas.GeoDataFrame(merged)
from mpl_toolkits.axes_grid1 import make_axes_locatable

fig, ax = plt.subplots(1, 1)

# ligne à ajouter pour avoir une légende ajustée à la taille du graphe
cax = make_axes_locatable(ax).append_axes("right", size="5%", pad=0.1)

geomerged.plot(column="rea", ax=ax, edgecolor='black', legend=True, cax=cax)
ax.set_title("Réanimations pour les %d départements" % metropole.shape[0]);
../_images/2020_carte_27_0.png

La création de carte a toujours été plus ou moins compliqué. Les premiers notebooks que j’ai créés sur le sujet étaient beaucoup plus complexe. geopandas a simplifié les choses. Son développement a commencé entre 2013 et a bien évolué depuis. Et j’ai dû passer quelques heures à récupérer les contours des départements il y a cinq ans.

On peut également récupérer la capacité maximale de chaque département en regardant sur le passé.

capacite = covid.groupby(["jour", "dep"]).sum().groupby("dep").max()
capacite.head()
sexe hosp rea rad dc
dep
01 3 904 78 2148 515
02 3 607 94 3173 749
03 3 550 58 1341 354
04 3 280 26 695 103
05 3 380 50 835 168
capa_merged = merged.merge(capacite, left_on="dep", right_on="dep")
capa_merged["occupation"] = capa_merged["rea_x"] / capa_merged["rea_y"]
capa_merged.head(n=2).T
0 1
dep 01 02
sexe_x 3 3
hosp_x 904 430
rea_x 52 54
rad_x 2148 3173
dc_x 515 749
code_depart 01 02
departement Ain Aisne
code_region 84 32
region Auvergne-Rhône-Alpes Hauts-de-France
code_ancien 82 22
ancienne_re Rhône-Alpes Picardie
geometry POLYGON ((4.795393450203793 46.22006712607377,... POLYGON ((3.120799460421226 49.70572360290051,...
sexe_y 3 3
hosp_y 904 607
rea_y 78 94
rad_y 2148 3173
dc_y 515 749
occupation 0.666667 0.574468
geocapa = geopandas.GeoDataFrame(capa_merged)
fig, ax = plt.subplots(1, 1)

# ligne à ajouter pour avoir une légende ajustée à la taille du graphe
cax = make_axes_locatable(ax).append_axes("right", size="5%", pad=0.1)

geocapa.plot(column="occupation", ax=ax, edgecolor='black', legend=True, cax=cax)
ax.set_title("Occupations en réanimations pour les %d départements" % metropole.shape[0]);
../_images/2020_carte_33_0.png

On voit que les premières régions touchées, l’Est et Paris paraissent les moins occupés mais c’est aussi sans doute parce que leur capacité ont dû augmenter durant le premier confinement. Il faudrait sans doute ne pas prendre le maximum de lit occupé mais une moyenne sur l’année précédente.