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

dev , python

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.

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

Figure 1 : 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")
Figure 2 : Frontières de la Belgique

Figure 2 : 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.

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

Figure 3 : 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

  • Folium pour générer des cartes leaflet en python.