Chapitre 3 : Manipuler et visualiser des données géographiques
Loin d’être exhaustif, ce tutoriel succint sur les données géographiques vise à présenter les bases pour manipuler et visualiser des données géographiques. Pour réaliser ce tutoriel, on s’est essentiellement appuyé sur la librarie geopandas dont la documentation officielle fournit d’ailleurs de nombreux exemples et explications.
Commençons par préalablement installer les packages utiles pour ce tutoriel :
!pip install geopandas mapclassify descartes pysal
Exercice 1 : Charger les données associées aux communes 2019 et celles associées à la population
import pandas as pd
communes=pd.read_csv("commune2019.csv")
pop=pd.read_excel("base-pop-historiques-1876-2015.xls", header=5, sheet_name="pop_1876_2015")
Exercice 2 : Créer un dataframe propre issu de la fusion des deux bases à partir du code commune puis agréger les données par département
data=pd.merge(communes, pop, left_on='com', right_on ="CODGEO")
data=data.groupby('dep', as_index=False)[['PMUN15', 'PMUN10']].sum()
Présentons les deux modes principaux de représentation numériques des données spatialisées :
- les images géoréférencées (photographies, images satellitaires …) aussi appelées rasters. Pour simplifier, un raster est une grille composée de cellules (ou pixels) organisées sous forme matricielle. Par exemple, une orthophotographie (photographie aérienne) est composée de 3 matrices empilées correspondant aux 3 couleurs primaires (rouges, vert, bleu).
- les couches vectorielles s’appuient sur le concept d’objets géométriques (points, lignes, polygones) pour représenter les entités géographiques. Un point est un doublet de coordonnées (longitude, latitude) auquel on rajoute parfois l’altitude. Une ligne est composée d’une suite de points. Un polygone est une ligne fermée. Afin de tenir compte des éventuelles enclaves et exclaves, les découpages territoriaux sont généralement décrits par des multipolygones qui sont des ensembles de polygones. On appelle données attributaires les variables non spatiales (par exemple le code géographique, le libellé, la région d’appartenance) qui décrivent les entités géographiques. En géopandas, la géométrie est intégrée sous la forme d’une colonne qui se rajoute aux données attributaires.
Vous pouvez télécharger les fonds de cartes (shapefile) sur le site creacartes.
Un shapefile est un format de stockage des données vectorielles. Il contient généralement les fichiers suivants :
- un fichier d’extension shp qui stocke les entités géographiques. Il s’agit du shapefile proprement-dit.
- un fichier d’extension dbf qui fichier qui contient les données attributaires relatives aux objets contenus dans le shapefile.
- un fichier d’extension shx qui stocke l’index de la géométrie
- un fichier d’extension * prj* qui stocke l’information sur le système de coordonnées
Pour comprendre les différents systèmes de projection, vous pouvez lire cet article.
Exercice 3 : Après avoir téléchargé le zip associé au shapefile dans l’environnement local, procéder au dézippage de ce dossier (à l’aide du package zipfile)
import zipfile
with zipfile.ZipFile('dep_francemetro_2018.zip', 'r') as zip_ref:
zip_ref.extractall()
Exercice 4 : Charger le fond de cartes (shapefile) des départements
import geopandas as gpd
dep = gpd.read_file('dep_francemetro_2018.shp')
dep.head()
Exercice 5 : Tracer la carte des départements à partir du shapefile précédent
dep.plot(figsize=(8, 8))
plt.axis('off')
La graduation automatique n’a pas de signification pour nos cartes. La commande plt.axis(‘off’) permet alors de la supprimer.
Exercice 6 : Tracer la carte des départements du la région Auvergne-Rhône-Alpes
dep[dep.reg=='84'].plot(figsize=(8, 8))
plt.axis('off')
Exercice 7 : Fusionner le shapefile avec la base de données
cartes = dep.merge(data, left_on = "code", right_on="dep")
print(cartes.shape)
cartes.head()
Exercice 8 : Tracer la carte de France indiquant la population par département
cartes.plot('PMUN10', legend=True, figsize=(8, 8))
plt.axis('off')
Exercice 9 : Pourquoi cette carte n’est-elle pas pertinente ?
Cette carte n’est pas pertinente pour deux raisons principales :
- Légende peu indicative : le changement de teinte en cartographie est généralement associé à un changement de signes (les couleurs chaudes sont associées aux valeurs positives)
- Cette carte ne réflète pas les différences de densité et est donc sémantiquement et sémiologiquement fausse. Une différence de population peut réfléter une différence d’aires.
Exercice 10 : Créer des variables relatives à la densité de la population en 2010 et 2015
cartes['densite_10'] = cartes.PMUN10 / (cartes.geometry.area / 10**6)
cartes['densite_15'] = cartes.PMUN15 / (cartes.geometry.area / 10**6)
Exercice 11 : Tracer la carte de densité de population par département. Il peut d’ailleurs être plus pertinent de tracer le log de la densité de population par département
cartes.plot('densite_10', legend=True, cmap='OrRd', scheme='Quantiles', figsize=(8, 8))
plt.axis('off')
Pour faciliter la visualisation, les variables quantitatives représentées sont discrétisées :
- soit par la méthode des quantiles (scheme='Quantiles’)
- soit par un découpage en intervalles de même taille (scheme='equal_interval’)
- soit à l’aide de l’algorithme de Jenks, c’est-à-dire un découpage en classes homogènes (scheme='Fisher_Jenks’).
En cartographique, on privilégie généralement la discrétisation selon l’algorithme de Jenks. Néanmoins, cette méthode peut s’avérer particluièrement lente sur des bases de données volumineuses. Dans ce cas, il est possible d’accélérer les calculs en discrétisant sur un échantillon avec l’option Fisher_Jenks_Sampled.
Exercice 12 : Calculer la différence de densité entre 2010 et 2015 puis calculer l’évolution en pourcentage de la densité de population
cartes['diff_15_10'] = ((cartes['densite_15']- cartes['densite_10']) )
cartes['evol_15_10'] = 100* cartes['diff_15_10'] / cartes['densite_10']
Exercice 13 : Tracer la carte de l’évolution de la densité de population par département entre 2010 et 2015.
cartes.plot('evol_15_10', legend=True, cmap='RdYlBu_r', scheme='Quantiles', figsize=(8, 8))
plt.axis('off')
Exercice 14 : Ajouter le numéro de département à la carte précédente
ax = cartes.plot('evol_15_10', legend=True, cmap='RdYlBu_r', scheme='Quantiles', figsize=(8, 8))
label=cartes.apply(lambda x: ax.annotate(s=x.code, xy=x.geometry.centroid.coords[0], ha='center'),axis=1)
plt.axis('off')
Exercice 15 : Agréger les données des départements pour obtenir celui des régions (incluant les données et la géométrie)
reg = cartes.dissolve(by='reg', aggfunc=sum).reset_index()
reg.plot(figsize=(8, 8))
plt.axis('off')
Exercice 16 : Représenter alors la densité de population par région en 2010
reg['densite_10'] = reg.PMUN10 / (reg.geometry.area / 10**6)
ax = reg.plot('densite_10', legend=True, cmap='OrRd', scheme='Quantiles', figsize=(8, 8), edgecolor='k')
label=reg.apply(lambda x: ax.annotate(s=x.reg, xy=x.geometry.centroid.coords[0], ha='center'),axis=1)
plt.axis('off')
Exercice 17 : Télécharger le zip associé au shapefile des fleuves français dans l’environnement local puis procéder au dézippage de ce dossier (à l’aide du package zipfile)
import zipfile
with zipfile.ZipFile('fleuve_francemetro_2018.zip', 'r') as zip_ref:
zip_ref.extractall()
Exercice 18 : Charger ce shapefile
fleuve = gpd.read_file('fleuve_francemetro_2018.shp')
fleuve.head(2)
Exercice 19 : Ne garder que les observations relatives à la Loire et tracer ce fleuve
loire = fleuve[fleuve.libgeo.str.lower().str.contains("loire")]
loire.plot(figsize=(8, 8))
plt.axis('off')
Exercice 20 : Tracer la carte des départements que la Loire traverse. Colorer les départements selon leur région d’appartenance.
Il faut commencer par installer le package rtree comme suit :
import os
os.environ['http_proxy'] = "http://proxy-rie.http.insee.fr:8080"
os.environ['https_proxy'] = "https://proxy-rie.http.insee.fr:8080"
!conda install -y -c conda-forge rtree
On peut alors évaluer l’intersection entre nos données et la Loire.
import rtree
dep_loire = gpd.sjoin(cartes, loire, op='intersects')
ax = dep_loire.plot(cmap='tab10', column='reg')
labels=dep_loire.apply(lambda x: ax.annotate(s=x.code, xy=x.geometry.centroid.coords[0], ha='center'),axis=1)
loire.plot(ax=ax, facecolor='none', edgecolor='k', figsize=(8, 8))
plt.axis('off')
Exercice 21 : Tracer la carte des départements qui se situent à moins de 10km autour de la Loire.
On sélectionne les données présentes dans le périmètre de 10 km autour de la Loire.
loire_buffer = gpd.geodataframe.GeoDataFrame(geometry=loire.geometry.buffer(10000))
dep_loire_buffer = gpd.sjoin(cartes, loire_buffer, op='intersects')
Puis, on représente la carte des départements concernés.
ax = dep_loire_buffer.plot(cmap='tab10', column='reg')
labels=dep_loire_buffer.apply(lambda x: ax.annotate(s=x.code, xy=x.geometry.centroid.coords[0], ha='center'),axis=1)
loire.plot(ax=ax, facecolor='none', edgecolor='k', figsize=(8, 8))
plt.axis('off')
Un petit aperçu de folium pour des cartes encore plus visuelles
L’objectif de folium est de réaliser des cartes dynamiques. Attention, contrairement aux manipulations précédentes, folium impose l’usage du système de projection WGS 84. Par défaut, les fonds de cartes de situation proposés sont issus d’OpenStreetMap.
Il faut commencer par installer et charger le package folium :
!pip install folium
import folium
Exercice 22 : Produire une carte dynamique centrée sur le bâtiment White (dont les coordonnées en WGS 84 sont [48.816207, 2.308189]).
m = folium.Map(location=[48.816207, 2.308189], zoom_start=19)
m
Exercice 23 : Creer une carte avec un marker sur notre lieu de travail.
m = folium.Map(
location=[48.816207, 2.308189],
zoom_start=18
)
folium.Marker([48.816285, 2.307647], popup='<i>SSP_Lab</i>', tooltip='Saurez-vous deviner notre lieu de travail?').add_to(m)
m
Exercice 24 : Représenter sur une carte les départements. (Attention : pour l’instant, on est encore obligé de transformer en geojson un fond de carte pour l’intégrer sous folium. Pour mémoire, il faut convertir en WGS 84 avant l’export en geoJson.)
cartes_json = cartes.to_crs({'init': 'epsg:4326'}).to_json()
m = folium.Map(
location=[48.816207, 2.308189],
zoom_start=6
)
folium.GeoJson(
cartes_json,
name='départements').add_to(m)
m
Voici une solution avec des otpions pour permettre l’affichage des noms des départements :
m = folium.Map(
location=[48.816207, 2.308189],
zoom_start=6
)
folium.GeoJson(
cartes_json,
name='départements',
show=True,
style_function=lambda feature: {
'fillColor': 'red',
'color': 'black',
'weight': 1,
'fillOpacity':0.1
},
highlight_function=lambda x: {'weight':3,
'color':'black',
'fillOpacity':1},
tooltip=folium.features.GeoJsonTooltip(
fields=['libelle'],
aliases=['Nom du département:'])
).add_to(m)
folium.LayerControl().add_to(m)
m
Exercice 25 : Représenter sur une carte l’évolution de la densité entre 2010 et 2015 par département.
bins = list(cartes['evol_15_10'].quantile([0, 0.25, 0.5, 0.75, 1]))
m = folium.Map(location=[48.816207, 2.308189], zoom_start=5)
folium.Choropleth(
geo_data=cartes_json,
data=cartes,
columns=['code', 'evol_15_10'],
key_on='feature.properties.code',
fill_color='RdYlBu_r',
fill_opacity=0.7,
line_opacity=0.5,
bins=bins,
reset=True
).add_to(m)
m