1A.data - Visualisation des données

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

Les tableaux et les graphes sont deux outils incontournables des statisticiens. Petite revue des graphes.

%matplotlib inline

Cette instruction fait apparaître les graphes dans le notebook. Si ce n’est pas le cas, il faut la réexécuer. Les deux lignes suivantes permettent de vérifier où matplotlib a prévu d’afficher ses résultats. Pour un notebook, cela doit être 'nbAgg' ou 'module://ipykernel.pylab.backend_inline'.

import matplotlib
matplotlib.get_backend()
'module://ipykernel.pylab.backend_inline'
import matplotlib.pyplot as plt
import matplotlib
matplotlib.get_backend()
'module://ipykernel.pylab.backend_inline'
from jyquickhelper import add_notebook_menu
add_notebook_menu()

Matplotlib, pandas

Récupération des données

On récupère les données disponibles sur le site de l’INSEE : Naissance, décès, mariages 2012. Il s’agit de récupérer la liste des mariages de l’année 2012. On souhaite représenter le graphe du nombre de mariages en fonction de l’écart entre les mariés.

from urllib.error import URLError
import pyensae
from pyensae.datasource import dBase2df, DownloadDataException
files = ["etatcivil2012_nais2012_dbase.zip",
         "etatcivil2012_dec2012_dbase.zip",
         "etatcivil2012_mar2012_dbase.zip" ]

try:
    pyensae.download_data(files[-1],
                          website='http://telechargement.insee.fr/fichiersdetail/etatcivil2012/dbase/')
except (DownloadDataException, URLError, TimeoutError):
    # backup plan
    pyensae.download_data(files[-1], website="xd")

df = dBase2df("mar2012.dbf")
df.shape, df.columns
((246123, 16),
 Index(['ANAISH', 'DEPNAISH', 'INDNATH', 'ETAMATH', 'ANAISF', 'DEPNAISF',
        'INDNATF', 'ETAMATF', 'AMAR', 'MMAR', 'JSEMAINE', 'DEPMAR', 'DEPDOM',
        'TUDOM', 'TUCOM', 'NBENFCOM'],
       dtype='object'))
df.head()
ANAISH DEPNAISH INDNATH ETAMATH ANAISF DEPNAISF INDNATF ETAMATF AMAR MMAR JSEMAINE DEPMAR DEPDOM TUDOM TUCOM NBENFCOM
0 1982 75 1 1 1984 99 2 1 2012 01 1 29 99 9 N
1 1956 69 2 4 1969 99 2 4 2012 01 3 75 99 9 N
2 1982 99 2 1 1992 99 1 1 2012 01 5 34 99 9 N
3 1985 99 2 1 1987 84 1 1 2012 01 4 13 99 9 N
4 1968 99 2 1 1963 99 2 1 2012 01 6 26 99 9 N

On récupère de la même manière la signification des variables :

from pyensae.datasource import dBase2df
vardf = dBase2df("varlist_mariages.dbf")
vardf.shape, vardf.columns
((16, 4), Index(['VARIABLE', 'LIBELLE', 'TYPE', 'LONGUEUR'], dtype='object'))
vardf
VARIABLE LIBELLE TYPE LONGUEUR
0 AMAR Année du mariage CHAR 4
1 ANAISF Année de naissance de l'épouse CHAR 4
2 ANAISH Année de naissance de l'époux CHAR 4
3 DEPDOM Département de domicile après le mariage CHAR 3
4 DEPMAR Département de mariage CHAR 3
5 DEPNAISF Département de naissance de l'épouse CHAR 3
6 DEPNAISH Département de naissance de l'époux CHAR 3
7 ETAMATF État matrimonial antérieur de l'épouse CHAR 1
8 ETAMATH État matrimonial antérieur de l'époux CHAR 1
9 INDNATF Indicateur de nationalité de l'épouse CHAR 1
10 INDNATH Indicateur de nationalité de l'époux CHAR 1
11 JSEMAINE Jour du mariage dans la semaine CHAR 1
12 MMAR Mois du mariage CHAR 2
13 NBENFCOM Enfants en commun avant le mariage CHAR 1
14 TUCOM Tranche de commune du lieu de domicile des époux CHAR 1
15 TUDOM Tranche d'unité urbaine du lieu de domicile de... CHAR 1

Exercice 1 : écart entre les mariés

  1. En ajoutant une colonne et en utilisant l’opération group by, on veut obtenir la distribution du nombre de mariages en fonction de l’écart entre les mariés. Au besoin, on changera le type d’une colone ou deux.
  2. On veut tracer un nuage de points avec en abscisse l’âge du mari, en ordonnée, l’âge de la femme. Il faudra peut-être jeter un coup d’oeil sur la documentation de la méthode plot.
# df["colonne"] = df.apply (lambda r:  int(r["colonne"]), axis=1)  # pour changer de type
# df["difference"] = ...

Exercice 2 : graphe de la distribution avec pandas

Le module pandas propose un panel de graphiques standard faciles à obtenir. On souhaite représenter la distribution sous forme d’histogramme. A vous de choisir le meilleure graphique depuis la page Visualization.

# df.plot(...)

matplotlib

matplotlib est le module qu’utilise pandas. Ainsi, la méthode plot retourne un objet de type Axes qu’on peut modifier par la suite via les méthodes suivantes. On peut ajouter un titre avec set_title ou ajouter une grille avec grid. On peut également superposer deux courbes sur le même graphique, ou changer de taille de caractères. Le code suivant trace le nombre de mariages par département.

df["nb"] = 1
dep = df[["DEPMAR","nb"]].groupby("DEPMAR", as_index=False).sum().sort_values("nb",ascending=False)
ax = dep.plot(kind = "bar", figsize=(14,6))
ax.set_xlabel("départements", fontsize=16)
ax.set_title("nombre de mariages par départements", fontsize=16)
ax.legend().set_visible(False)  # on supprime la légende

# on change la taille de police de certains labels
for i,tick in enumerate(ax.xaxis.get_major_ticks()):
    if i > 10 :
        tick.label.set_fontsize(8)
../_images/td1a_cenonce_session_12_18_0.png

Quand on ne sait pas, le plus simple est d’utiliser un moteur de recherche avec un requête du type : matplotlib + requête. Pour créer un graphique, le plus courant est de choisir le graphique le plus ressemblant d’une gallerie de graphes puis de l’adapter à vos données.

Exercice 3 : distribution des mariages par jour

On veut obtenir un graphe qui contient l’histogramme de la distribution du nombre de mariages par jour de la semaine et d’ajouter une seconde courbe correspond avec un second axe à la répartition cumulée.

Réseaux, graphes

networkx

Le module networkx permet de représenter un réseau ou un graphe de petite taille (< 500 noeuds). Un graphe est défini par un ensemble de noeuds (ou vertex en anglais) reliés par des arcs (ou edge en anglais). La gallerie vous donnera une idée de ce que le module est capable de faire.

import random
import networkx as nx
G=nx.Graph()
for i in range(15) :
    G.add_edge ( random.randint(0,5),  random.randint(0,5) )

import matplotlib.pyplot as plt
f, ax = plt.subplots(figsize=(8,4))
nx.draw(G, ax = ax)
../_images/td1a_cenonce_session_12_22_0.png

Graphviz

Graphviz est un outil développé depuis plusieurs années déjà qui permet de réprésenter des graphes plus conséquents (> 500 noeuds). Il propose un choix plus riche de graphes : gallerie. Il est utilisable via le module graphviz. Son installation requiert l’installation de l’outil Graphviz qui n’est pas inclus. La différence entre les deux modules tient dans l’algorithme utilisé pour assigner des coordonnées à chaque noeud du graphe de façon à ce que ses arcs se croisent le moins possibles. Au delà d’une certaine taille, le dessin de graphe n’est plus lisible et nécessite quelques tatônnements. Cela peut passer par une clusterisation du graphe (voir la méthode Louvain) de façon à colorer certains noeuds proches voire à les regrouper.

import random, os
from graphviz import Digraph
from IPython.display import Image
from pyquickhelper.helpgen import find_graphviz_dot
bin = os.path.dirname(find_graphviz_dot())
if bin not in os.environ["PATH"]:
    os.environ["PATH"] = os.environ["PATH"] + ";" + bin

dot = Digraph(comment='random graph', format="png")
for i in range(15) :
    dot.edge ( str(random.randint(0,5)),  str(random.randint(0,5)) )

img = dot.render('t_random_graph.gv')
Image(img)
../_images/td1a_cenonce_session_12_24_0.png

Exercice 4 : dessin d’un graphe avec networkx

On construit un graphe aléatoire, ses 20 arcs sont obtenus en tirant 20 fois deux nombres entiers entre 1 et 10. Chaque arc doit avoir une épaisseur aléatoire. On regardera les fonctions spring_layout, draw_networkx_nodes, draw_networkx_edges. La gallerie peut aider aussi.

Cartographie

Dessiner une carte n’est pas difficile en soit. Toute la difficulté vient du fait qu’on a besoin pour lire cette carte de point de référence : les rues pour une ville, les frontières pour un département, une région, un pays, les fleuves et montagnes pour une carte représentation des données démographiques. Ces informations sont importantes pour situer l’information représentée par le graphique. Si vous n’avez pas besoin de tout ça, les formules suivantes vous suffiront :

Ces fonctionnalités sont disponibles via le module geopy. Dans le cas contraire, voici quelques directions :

  • cartopy : le module est une surcouche du module matplotlib, il convertit les coordonnées géographiques et ajoutent quelques éléments de paysages (frontières, rivières…). Il ne contient pas tout comme la définition des départements français.
  • shapely : ce module est utile pour dessiner des aires sur des cartes. Sous Windows, il faut l’installer depuis Unofficial Windows Binaries for Python Extension Packages car il inclut la DLL geos_c.dll qui vient de GEOS. Dans le cas contraire, il faut installer GEOS, ce qui prend pas mal de temps. Il est utilisé par cartopy.

Il en existe d’autres mais leur installation recèle quelques difficultés que je n’ai pas toujours eu la patience de contourner :

  • mapnik : l’installation sur Windows est réservée aux connaisseurs, trop compliquée sous Windows.
  • geopandas manipule les coordonnées géographiques dans des dataframe et implémente des graphiques simples.

L’exemple qui suit utilise cartopy. Comme beaucoup de modules, il contient principalement des informations que le territoire américains.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
projection = ccrs.PlateCarree()
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(1, 2, 1, projection=projection)
ax.add_feature(cfeature.BORDERS)
ax.add_feature(cfeature.LAKES)
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.RIVERS)
ax.add_feature(cfeature.COASTLINE)
ax.set_extent([-5, 9, 42, 52])
ax.set_title("France")

ax = fig.add_subplot(1, 2, 2, projection=projection)
ax.add_feature(cfeature.BORDERS)
ax.add_feature(cfeature.LAKES)
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.RIVERS)
ax.add_feature(cfeature.COASTLINE)
ax.set_extent([-125, -66.5, 20, 50])
ax.set_title("USA");
../_images/td1a_cenonce_session_12_27_0.png

On peut améliorer la précision en précisant l’échelle avec la méthode with_scale qui propose trois résolution 10m, 50m, 110m. La module cartopy n’inclut que la résolution 110m, le reste doit être téléchargé. La première exécution prend un peu de temps.

projection = ccrs.PlateCarree()
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(1, 2, 1, projection=projection)
ax.add_feature(cfeature.BORDERS.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.set_extent([-5, 9, 42, 52])
ax.set_title("France")

ax = fig.add_subplot(1, 2, 2, projection=projection)
ax.add_feature(cfeature.BORDERS.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.set_extent([-125, -66.5, 20, 50])
ax.set_title("USA");
../_images/td1a_cenonce_session_12_29_0.png

On peut ajouter un peu plus d’information avec la méthode stock_img.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
projection = ccrs.PlateCarree()
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(1, 2, 1, projection=projection)
ax.add_feature(cfeature.BORDERS.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.set_extent([-5, 9, 42, 52])
ax.set_title("France")
ax.stock_img()

ax = fig.add_subplot(1, 2, 2, projection=projection)
ax.add_feature(cfeature.BORDERS.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.set_extent([-125, -66.5, 20, 50])
ax.stock_img()
ax.set_title("USA");
../_images/td1a_cenonce_session_12_31_0.png

Cette fonction plaque une portion d’une grande image représentant la terre dans cette projection. On peut en obtenir d’autre à des résolutions différentes sur 1:10m Natural Earth I ou encore 1:10m Gray Earth mais ces images sont conséquentes. On s’inspire du code de la fonction stock_img et on écrit le code qui suit. D’autres informations sont disponibles sur github/natural-earth-vector.

from matplotlib.image import imread

projection = ccrs.PlateCarree()
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(1, 2, 1, projection=projection)
ax.add_feature(cfeature.BORDERS.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.set_extent([-5, 9, 42, 52])
ax.set_title("France 50")
ax.stock_img()

ax = fig.add_subplot(1, 2, 2, projection=projection)

fname = "NE1_LR_LC.tif"
if os.path.exists(fname):
    ax.imshow(imread(fname), origin='upper', transform=projection,
              extent=[-180, 180, -90, 90]);
else:
    import warnings
    warnings.warn("L'image n'a pas été téléchargée. Lire le paragraphe qui précède.")

ax.add_feature(cfeature.BORDERS.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.set_extent([-5, 9, 42, 52])
ax.set_title("France 10");
c:python370_x64libsite-packagesPILImage.py:2546: DecompressionBombWarning: Image size (131220000 pixels) exceeds limit of 89478485 pixels, could be decompression bomb DOS attack.
  DecompressionBombWarning)
../_images/td1a_cenonce_session_12_33_1.png

Presque parfait. Le suivant montre l’Europe et ses pays puis ajoute les méridiens et parallèles.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
projection = ccrs.PlateCarree()
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(1, 1, 1, projection=projection)
ax.add_feature(cfeature.BORDERS)
ax.add_feature(cfeature.LAKES)
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.RIVERS)
ax.add_feature(cfeature.COASTLINE)
ax.set_extent([-10, 31, 35, 71])
ax.stock_img()
ax.set_title("Europe")

x, y = 2.3488000, 48.853410
ax.plot(x, y, 'bo', markersize=6)
ax.text(x, y, "Paris")

ax.gridlines(crs=projection, draw_labels=True,
             linewidth=2, color='gray', alpha=0.5, linestyle='--');
../_images/td1a_cenonce_session_12_35_0.png

Cartes avec les départements

Pour dessiner des formes sur une carte, il faut connaître les coordonnées de ces formes. La première chose à faire est de récupérer leurs données géographiques. Une fois simple de les trouver est d’utiliser un moteur de recherche avec le mot clé shapefile inclus dedans : c’est le format du fichier. shapefile france permet d’obtenir quelques sources. En voici d’autres :

La première chose à vérifier est la licence associées aux données : on ne peut pas en faire ce qu’on veut. Pour cet exemple, j’ai choisi la première source de données, GADM. La licence n’est pas précisée explicitement (on peut trouver happy to share sur le site, la page wikipedia GADM précise : GADM is not freely available for commercial use. The GADM project created the spatial data for many countries from spatial databases provided by national governments, NGO, and/or from maps and lists of names available on the Internet (e.g. from Wikipedia). C’est le choix que j’avais fait en 2015 mais l’accès à ces bases a probablement changé car l’accès est restreint. J’ai donc opté pour les bases accessibles depuis data.gouv.fr. Leur seul inconvénient est que les coordonnées sont exprimées dans une projection de type Lambert 93. Cela nécessite une conversion.

from pyensae import download_data
try:
    download_data("GEOFLA_2-1_DEPARTEMENT_SHP_LAMB93_FXX_2015-12-01.7z",
                  website="https://wxs-telechargement.ign.fr/oikr5jryiph0iwhw36053ptm/telechargement/inspire/" + \
                          "GEOFLA_THEME-DEPARTEMENTS_2015_2$GEOFLA_2-1_DEPARTEMENT_SHP_LAMB93_FXX_2015-12-01/file/")
except Exception as e:
    # au cas le site n'est pas accessible
    download_data("GGEOFLA_2-1_DEPARTEMENT_SHP_LAMB93_FXX_2015-12-01.7z", website="xd")
from pyquickhelper.filehelper import un7zip_files
try:
    un7zip_files("GEOFLA_2-1_DEPARTEMENT_SHP_LAMB93_FXX_2015-12-01.7z", where_to="shapefiles")
    departements = 'shapefiles\\GEOFLA_2-1_DEPARTEMENT_SHP_LAMB93_FXX_2015-12-01\\GEOFLA\\1_DONNEES_LIVRAISON_2015\\' + \
                   'GEOFLA_2-1_SHP_LAMB93_FR-ED152\\DEPARTEMENT\\DEPARTEMENT.shp'
except FileNotFoundError as e:
    # Il est possible que cette instruction ne fonctionne pas.
    # Dans ce cas, on prendra une copie de ce fichier.
    import warnings
    warnings.warn("Plan B parce que " + str(e))
    download_data("DEPARTEMENT.zip")
    departements = "DEPARTEMENT.shp"

if not os.path.exists(departements):
    raise FileNotFoundError("Impossible de trouver '{0}'".format(departements))

La license accompagne les données : ce produit est téléchargeable et utilisable gratuitement sous licence Etalab. Pour un usage commercial, il faut faire attentation à la license associée aux données. Le seul inconvénient des données GEOFLA est que certaines sont données dans le système de coordonnées Lambert 93 (voir aussi Cartographie avec R).

shp = departements
import shapefile
r = shapefile.Reader(shp)
shapes = r.shapes()
records = r.records()
len(shapes), len(records)
(96, 96)
r.measure
[0.0, 0.0]
r.bbox
[99217.1, 6049646.300000001, 1242417.2, 7110480.100000001]
r.shapeType
5

On regarde une zone en particulier mais on réduit la quantité de données affichées :

d = shapes[0].__dict__.copy()
d["points"] = d["points"][:10]
d
{'shapeType': 5,
 'points': [(701742.0, 6751181.100000001),
  (701651.9, 6751166.9),
  (701552.0, 6751162.7),
  (700833.7000000001, 6751313.7),
  (700669.4, 6751380.0),
  (700475.4, 6751476.600000001),
  (700400.7000000001, 6751517.2),
  (700098.3, 6751789.600000001),
  (699993.8, 6751845.4),
  (699874.1000000001, 6751876.4)],
 'parts': [0],
 'bbox': [688654.4, 6690595.300000001, 800332.3, 6811114.5]}
records[0]
['DEPARTEM0000000000000004',
 '89',
 'YONNE',
 '024',
 'AUXERRE',
 742447,
 6744261,
 748211,
 6750855,
 '27',
 'BOURGOGNE-FRANCHE-COMTE']

J’utilise une fonction de conversion des coordonnées récupérée sur Internet.

import math

def lambert932WGPS(lambertE, lambertN):

    class constantes:
        GRS80E = 0.081819191042816
        LONG_0 = 3
        XS = 700000
        YS = 12655612.0499
        n = 0.7256077650532670
        C = 11754255.4261

    delX = lambertE - constantes.XS
    delY = lambertN - constantes.YS
    gamma = math.atan(-delX / delY)
    R = math.sqrt(delX * delX + delY * delY)
    latiso = math.log(constantes.C / R) / constantes.n
    sinPhiit0 = math.tanh(latiso + constantes.GRS80E * math.atanh(constantes.GRS80E * math.sin(1)))
    sinPhiit1 = math.tanh(latiso + constantes.GRS80E * math.atanh(constantes.GRS80E * sinPhiit0))
    sinPhiit2 = math.tanh(latiso + constantes.GRS80E * math.atanh(constantes.GRS80E * sinPhiit1))
    sinPhiit3 = math.tanh(latiso + constantes.GRS80E * math.atanh(constantes.GRS80E * sinPhiit2))
    sinPhiit4 = math.tanh(latiso + constantes.GRS80E * math.atanh(constantes.GRS80E * sinPhiit3))
    sinPhiit5 = math.tanh(latiso + constantes.GRS80E * math.atanh(constantes.GRS80E * sinPhiit4))
    sinPhiit6 = math.tanh(latiso + constantes.GRS80E * math.atanh(constantes.GRS80E * sinPhiit5))

    longRad = math.asin(sinPhiit6)
    latRad = gamma / constantes.n + constantes.LONG_0 / 180 * math.pi

    longitude = latRad / math.pi * 180
    latitude = longRad / math.pi * 180

    return longitude, latitude

lambert932WGPS(99217.1, 6049646.300000001), lambert932WGPS(1242417.2, 7110480.100000001)
((-4.1615802638173065, 41.303505287589545),
 (10.699505053975292, 50.85243395553585))

Après plusieurs recherches sur internet, un grand nombre, on aboutit à la carte suivante non sans avoir redéfini des contours et frontières plus précis inspiré du code feature.py.

from cartopy.feature import NaturalEarthFeature, COLORS
resolution = "50m"
BORDERS = NaturalEarthFeature('cultural', 'admin_0_boundary_lines_land',
                              resolution, edgecolor='black', facecolor='none')
STATES = NaturalEarthFeature('cultural', 'admin_1_states_provinces_lakes',
                             resolution, edgecolor='black', facecolor='none')
COASTLINE = NaturalEarthFeature('physical', 'coastline', resolution,
                                edgecolor='black', facecolor='none')
LAKES = NaturalEarthFeature('physical', 'lakes', resolution,
                            edgecolor='face',
                            facecolor=COLORS['water'])
LAND = NaturalEarthFeature('physical', 'land', resolution,
                           edgecolor='face',
                           facecolor=COLORS['land'], zorder=-1)
OCEAN = NaturalEarthFeature('physical', 'ocean', resolution,
                            edgecolor='face',
                            facecolor=COLORS['water'], zorder=-1)
RIVERS = NaturalEarthFeature('physical', 'rivers_lake_centerlines', resolution,
                             edgecolor=COLORS['water'],
                             facecolor='none')
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
projection = ccrs.PlateCarree()
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(1, 1, 1, projection=projection)
ax.add_feature(BORDERS)
ax.add_feature(LAKES)
ax.add_feature(LAND)
ax.add_feature(OCEAN)
ax.add_feature(RIVERS)
ax.add_feature(COASTLINE)
ax.set_extent([-5, 12, 40, 54])
ax.set_title("Europe")
ax.gridlines(crs=projection, draw_labels=True,
             linewidth=2, color='gray', alpha=0.5, linestyle='--')

from matplotlib.collections import LineCollection
import shapefile
import geopandas
from shapely.geometry import Polygon
from shapely.ops import cascaded_union, unary_union

shp = departements
r = shapefile.Reader(shp)
shapes = r.shapes()
records = r.records()

polys = []
for i, (record, shape) in enumerate(zip(records, shapes)):
    # les coordonnées sont en Lambert 93
    if i == 0:
        print(record, shape.parts)
    geo_points = [lambert932WGPS(x,y) for x, y in shape.points]
    if len(shape.parts) == 1:
        # Un seul polygone
        poly = Polygon(geo_points)
    else:
        # Il faut les fusionner.
        ind = list(shape.parts) + [len(shape.points)]
        pols = [Polygon(geo_points[ind[i]:ind[i+1]]) for i in range(0, len(shape.parts))]
        try:
            poly = unary_union(pols)
        except Exception as e:
            print("Cannot merge: ", record)
            print([_.length for _ in pols], ind)
            poly = Polygon(geo_points)
    polys.append(poly)

data = geopandas.GeoDataFrame(geometry=polys)
# cmap -> voir https://matplotlib.org/users/colormaps.html
data.plot(ax=ax, cmap='tab20', edgecolor='black');
# Ou pour définir des couleurs spécifiques.
# geopandas.plotting.plot_polygon_collection(ax, data['geometry'], facecolor=data['colors'], values=None)
['DEPARTEM0000000000000004', '89', 'YONNE', '024', 'AUXERRE', 742447, 6744261, 748211, 6750855, '27', 'BOURGOGNE-FRANCHE-COMTE'] [0]
../_images/td1a_cenonce_session_12_51_1.png

Cartes interactives

La vidéo Spatial data and web mapping with Python vous en dira un peu plus sur la cartographie. Lorsqu’on dessine une carte avec les rues d’une ville, on veut pouvoir zoomer et dézoomer facilement pour avoir une vue d’ensemble ou une vue détaillé. Dans ce cas là, il faut utiliser un service externe telle que Gmap API, Bing Map API, Yahoo Map API ou OpenStreetMap qui est une version open source. Dans tous les cas, il faut faire attention si les données que vous manipulez dans la mesure où elles transitent par un service externe. L’article Busy areas in Paris est un exemple d’utilisation d”OpenStreetMap. Ce qu’on cherche avant tout, c’est un graphique interactif. Il existe des modules qui permettent d’utiliser ces services directement depuis un notebook python. smopy crée une carte non interactive :

import smopy
map = smopy.Map((42., -1., 55., 3.), z=5)
map.show_ipython()
../_images/td1a_cenonce_session_12_54_0.png

Le module folium insert du javascript dans le notebook lui-même. Voici un exemple construit à partir de ce module : Creating Interactive Election Maps Using folium and IPython Notebooks. Il est prévu pour fonctionner comme suit. D’abord, une étape d’initialisation :

import folium

Et si elle fonctionne (un jour peut-être), la suite devrait être :

map_osm = folium.Map(location=[48.85, 2.34])
map_osm

Donc, on prend un raccourci et on en profite pour ajouter un triangle à l’emplacement de l’ENSAE :

import folium
map_osm = folium.Map(location=[48.711478,2.207708])
from pyensae.notebookhelper import folium_html_map
map_osm.add_child(folium.RegularPolygonMarker(location=[48.711478,2.207708], popup='ENSAE',
                       fill_color='#132b5e', radius=10))
from IPython.display import HTML
#HTML(folium_html_map(map_osm))
map_osm

Un moteur de recherche vous donnera rapidement quelques exemples de cartes utilisant ce module : Folium Polygon Markers, Easy interactive maps with folium, Plotting shapefiles with cartopy and folium.

Graphiques interactifs

Ce sujet sort du cadre de cette séance. Un graphique interactif réagit à des événéments comme le passage du curseur de la souris, un clic, un zoom. Il est difficile de lire un graphique trop chargé, c’est pourquoi en rendant le graphique interactifs, on cherche à donner plus d’information sans nuire à sa lisibilité. Voici quelques scénarios :

  • On veut représenter une des dimensions du problème en animant le graphique. C’est fréquent en 3D où un des axes est celui du temps. On préfèrera un graphique en 2D évoluant dans le temps.
  • Lorsqu’il y a trop de courbes à dessiner, le lecteur peut activer ou désactiver certaines courbes pour pouvoir les comparer.
  • On peut permettre de changer d’échelle (logarithmique ou changer la base 100 à différents endroits).
  • On veut donner une vue d’ensemble et en même temps un niveau de détails plus fin si le lecteur le demande.

Ces animations pris leur essor avec internet et le langage javascript. Concevoir un graphique animé nécessite plus de temps car il faut prévoir ce qu’il doit se passer en cas d’action du lecteur (souris, touche, …). Les modules python permettant de les créer construisent en fait un code javascript qu’il faut ensuite exécuter dans un navigateur (ou directement dans un notebook comme celui-ci). La librairie javascript qui a changé la façon de les concevoir est d3.js. Beaucoup d’autres librairies sont des surcouches de celle-ci nvd3.

Le module plotly est intéressant mais toutes les fonctionnalités ne sont pas gratuites. De plus, il faut créer un compte pour pouvoir s’en servir. Le module mpld3 est intéressant dans le sens où il convertit un graphique créé avec un matplotlib en un graphique javascript réalisé avec d3.js. Il faut d’abord activer la sortie notebook (voir D3 Plugins: Truly Interactive Matplotlib In Your Browser). Le module n’est plus maintenu et il vaut mieux installer la version disponible sur GitHub/mpld3.

import mpld3
mpld3.enable_notebook()

Ensuite, les graphiques réalisés avec matplotlib seront affichés en javascript (si cela ne fonctionne pas, il ne faut pas hésiter à redémarrer le Kernel). Il faut cliquer sur le graphique pour que celui-ci deviennent zoomable.

import matplotlib.pyplot as plt
import numpy as np
fig, ax = plt.subplots()
ax.grid(color='white', linestyle='solid')

N = 50
scatter = ax.scatter(np.random.normal(size=N),
                     np.random.normal(size=N),
                     c=np.random.random(size=N),
                     s = 1000 * np.random.random(size=N),
                     alpha=0.3,
                     cmap=plt.cm.jet)

ax.set_title("D3 Scatter Plot", size=18);

Le module mpld3 n’est plus maintenu par son auteur qui a basculé sur altair. Il n’existe pas encore un module incontournable. Un de ceux qui pourrait devenir une référence est le module bokeh. Il n’utilise pas d3.js mais le principe est le même. Il existe aussi pyecharts.

import bokeh, bokeh.io as bio
bio.output_notebook()
Loading BokehJS ...

Une fois que le module est initialisé, on peut afficher son graphique.

import bokeh.plotting as bplt
p = bplt.figure(title = "Exemple")
p.xaxis.axis_label = 'X'
p.yaxis.axis_label = 'Y'

p.circle([1,2,3,4],[4,5,6,5.5], fill_color="red", color="red", size=12)
bplt.show(p)
from pyquickhelper.helpgen import NbImage  # seconde image
NbImage("pngbokeh.png")            # pour la conversion des notebooks au format HTML
../_images/td1a_cenonce_session_12_70_0.png

Pour sauver le graph sous forme de fichier HTML :

import os
bplt.output_file("example_bokeh.html")
bplt.save(p)
print([ _ for _ in os.listdir(".") if "html" in _ ] )
['example_bokeh.html']