Sèkun Blog

Création d'un scatter plot de coordonnées avec matplotlib sur une carte openstreetmap

Sekun   29-03-2020 11:18 informatique

L'objectif va être de réprésenter des coordonnées gps sur une carte quelconque, ici on va utiliser une carte de la Belgique, mais le code s'adapte pour tout autre pays.

Carte de la belgique avec des plusieurs coordonnées en bleu situées dessus

Installation

Les librairies utilisées sont principalement:

  • Pandas l'outil de manipulation de données en python.
  • Geopandas afin de représenter des données géo-spatiales, basé sur pandas.
  • Matplotlib comme outil de création de graphe dont la syntaxe est inspirée par Matlab (utilisé par pandas).
  • Contextily pour rajouter un fond en provenance de Openstreetmap.
  • Shapely pour convertir des données d'une projection à une autre.

Le plus simple pour mettre en place l'environnement de travail est d'utiliser Anaconda.

On peut alors installer les dépendances nécessaires, en créant un nouvel environnement virtuel via conda:

conda create -n drawit matplotlib geopandas pandas shapely
conda activate drawit
pip install contextily==1.0rc2

Les imports utilisés sont les suivants:

import pandas as pd
import matplotlib.pyplot as plt
import geopandas as gpdd
import json
from random import choices
from shapely.geometry import Point

Données externes

Pour la carte de la Belgique, on va utiliser un fichier GeoJSON ou shp (ESRI Shapefile), que Geopandas peut lire. Le fichier utilisé est le suivant: BEL_adm0.shp, et provient de ce repo github.

Pour les coordonnées GPS à représenter, au lieu de les générer aléatoirement dans les bornes délimitant la Belgique, on peut plutôt utiliser la liste des coordonnées des villes du pays, obtenue via le ce repo github.

belgium = gp.read_file("BEL_adm_shp/BEL_adm0.shp")
with open("zipcode-belgium.json") as f:
    gps = json.load(f)

On peut déjà afficher le rendu du fichier shp:

belgium.plot(facecolor="none", edgecolor="gray")
plt.axis("off")
Frontières de la Belgique

Les données gps obtenues contiennent l'ensemble des villes, leur code postal ainsi que leurs coordonnées gps.

[{'zip': '1000', 'city': 'Bruxelles', 'lng': 4.351697, 'lat': 50.8465573},
 {'zip': '1020', 'city': 'Laeken', 'lng': 4.3487134, 'lat': 50.883392},
 {'zip': '1030', 'city': 'Schaerbeek', 'lng': 4.3737121, 'lat': 50.8676041},
 ...

Représentation des coordonnées sur la carte

Nous allons d'abord prendre au hasard 50 données parmi ces dernières, que nous instancions via la classe DataFrame de pandas.

data = pd.DataFrame(choices(gps, k=50))

Il nous reste à représenter ces coordonnées sur la carte.

# Création d'une seule figure de taille 10x10
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(1,1,1)
plt.axis("off") # masquer les axes
belgium.plot(ax=ax, facecolor="none", edgecolor="gray", linewidth=0.4)
data.plot(ax=ax, x="lng", y="lat", kind='scatter', alpha=0.7)
plt.show()

Le facecolor="none" est la façon avec Matplotlib de rendre le fond de la carte transparent.

Frontières de la Belgique avec des coordonnées dessinées dessus

Openstreetmap pour le fond de la carte

Nous allons rajouter comme fond une carte en provenance de Openstreetmap. Pour ce faire, nous allons utiliser la librairie contextily, conseillée dans la doc de geopandas.

Le problème est que les coordonnées GPS fournies, et celles utilisées dans les cartes Openstreetmap ne sont pas obtenues via la même projection.

Pour Openstreetmap on a une projection «Spherical Mercator» (EPSG:3857), pour les coordonnées GPS: «World Geodetic System» (EPSG:4326). Nous allons devoir convertir les données de la variable belgium, ainsi que celle de gps de EPSG:4326 vers EPSG:3857. (Nous allons utiliser la librairie Shapely).

belgium = belgium.to_crs(epsg=3857)
geometry = data2["coords"].apply(Point)
data = gpd.GeoDataFrame(data, crs={"init": 'epsg:4326'}, geometry=geometry).to_crs(epsg=3857)

Reste à rajouter une carte comme fond ; nous avons plusieurs choix fournis par la librairie contextily. Celle-ci utilise en réalité une librairie basée sur leaflet. On choisit un modèle sur cette page, et le nom peut être récupéré via le fichier _providers.py de contextily généré automatiquement.

On utilise simplement la méthode add_basemap de contextily.

fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(1,1,1)
plt.axis("off")
belgium.plot(ax=ax, facecolor="none", edgecolor="black", linewidth=1, alpha=0.8)
data2.geometry.plot(ax=ax, alpha=0.5)
ctx.add_basemap(ax, url=ctx.providers.Hydda.Base, alpha=0.8)

Autres outils

Des commentaires ? J'y répondrai sur mastodon .

Ce texte est publié sous la licence cc-by-sa.