Tracer une carte en Python#

Le notebook propose plusieurs façons de tracer une carte en Python.

Il y a principalement trois façons de tracer une carte. La première est statique avec des modules comme basemap ou cartopy qui sont des surcouches de matplotlib. Le second moyen est une carte animée ou non dans un notebook avec des modules comme pygal, plotly. La dernière consiste à insérer des éléments sur une carte en ligne telle que OpenStreetMap et le module folium ou ipyleaflet.

Il y a souvent trois problèmes avec les cartes. Le premier sont avec les coordonnées. Les plus utilisées sont les coordonnées longitude / latitude. Le problème est que chaque pays a son propre système adapté à sa position géographique. Il faut souvent convertir (voir lambert93_to_WGPS, pyproj). Le second problème est l’ajout de repères géographiques (frontières, fleuves, …). Certains modules contiennent certaines informations, souvent pour les Etats-Unis. Mais souvent, il faut récupérer ces informations sur les sites open data de chaque pays. La troisième difficulté est qu’on veut tracer des cartes très chargées et cela prend un temps fou pour construire une carte lisible.

[2]:

%matplotlib inline

données#

[2]:

from teachpyx.datasets import load_enedis_dataset

df = load_enedis_dataset()
df.head(n=2).T

[2]:
0 1
Année 2016 2016
Nom commune Ponteilla Varreddes
Code commune 66145 77483
Nom EPCI CU Perpignan Méditerranée (Pmcu) CA Pays de Meaux
Code EPCI 200027183 247700628
Type EPCI CU CA
Nom département Pyrénées-Orientales Seine-et-Marne
Code département 66 77
Nom région Occitanie Île-de-France
Code région 76 11
Domaine de tension BT > 36 kVA BT <= 36 kVA
Nb sites Photovoltaïque Enedis 73.0 10.0
Energie produite annuelle Photovoltaïque Enedis (MWh) 10728.620366 21.41684
Nb sites Eolien Enedis 0.0 0.0
Energie produite annuelle Eolien Enedis (MWh) 0.0 0.0
Nb sites Hydraulique Enedis 0.0 0.0
Energie produite annuelle Hydraulique Enedis (MWh) 0.0 0.0
Nb sites Bio Energie Enedis 0.0 0.0
Energie produite annuelle Bio Energie Enedis (MWh) 0.0 0.0
Nb sites Cogénération Enedis 0.0 0.0
Energie produite annuelle Cogénération Enedis (MWh) 0.0 0.0
Nb sites Autres filières Enedis 0.0 0.0
Energie produite annuelle Autres filières Enedis (MWh) 0.0 0.0
Geo Point 2D 42.6323626522, 2.82631103755 49.0059497861, 2.92725176893
long 2.826311 2.927252
lat 42.632363 49.00595

cartopy#

basemap est l’ancêtre des modules de tracé de cartes sous Python. Ce notebook utilise cartopy. Contrairement à basemap, cartopy n’installe pas toutes les données dont il a besoin mais télécharge celle dont il a besoin pour tracer une carte. La projection indique comment la surface de la terre, sphérique, sera projetée dans le plan. Ensuite, il faut un système de coordonnées pour localiser un point sur la surface. Le plus utilisée est WGS_84 ou longitude, latitude. En France, l”INSEE utilise aussi le système Lambert 93 ou EPSG 2154. Source : Introduction à la manipulation de données cartographiques. Tout n’est pas parfait dans Cartopy comme ce problème Create Cartopy projection from pyproj.Proj.

[3]:

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_extent([-5, 10, 42, 52])

ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=":")
ax.plot(df.long, df.lat, ".")
ax.set_title("France");

/home/xadupre/.local/lib/python3.10/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_physical/ne_10m_ocean.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)
/home/xadupre/.local/lib/python3.10/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_physical/ne_10m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)
/home/xadupre/.local/lib/python3.10/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_cultural/ne_10m_admin_0_boundary_lines_land.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)
../_images/c_data_enedis_cartes_6_1.png

On peut obtenir une carte un peu plus détaillée mais le module cartopy télécharge des données pour ce faire. Cela se faire avec l’instruction with_scale qui propose trois résolutions :10m, 50m, 110m.

[4]:

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_extent([-5, 10, 42, 52])

ax.add_feature(cfeature.OCEAN.with_scale("50m"))
ax.add_feature(cfeature.COASTLINE.with_scale("50m"))
ax.add_feature(cfeature.RIVERS.with_scale("50m"))
ax.add_feature(cfeature.BORDERS.with_scale("50m"), linestyle=":")
ax.plot(df.long, df.lat, ".")
ax.set_title("France");

/home/xadupre/.local/lib/python3.10/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_physical/ne_50m_ocean.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)
/home/xadupre/.local/lib/python3.10/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_physical/ne_50m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)
/home/xadupre/.local/lib/python3.10/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_physical/ne_50m_rivers_lake_centerlines.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)
/home/xadupre/.local/lib/python3.10/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_cultural/ne_50m_admin_0_boundary_lines_land.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)
../_images/c_data_enedis_cartes_8_1.png

On ajoute un système de coordonnées français particulièrement intéressant pour la France. On convertit d’abord longitude, latitude en Lambert 93.

[11]:

from pyproj import Proj, Transformer

transform = Transformer.from_crs(
    "EPSG:4326", "EPSG:2154"  # longitude / latidude  # Lambert 93
)
transform.transform(-5, 42)

[11]:

(6634541.621546051, 1633374.3871915536)

[12]:

transform.transform(10, 52)

[12]:

(6723037.367738617, 4228679.19768552)

On convertit toutes les coordonnées.

[13]:

lamb_x, lamb_y = transform.transform(df.long.values, df.lat.values)

Et on dessine deux cartes, la première en longitude, latitude, la seconde en Lambert 93.

[24]:

import cartopy.crs as ccrs
from cartopy.crs import CRS, Globe
import cartopy.feature as cfeature
import matplotlib.pyplot as plt


def parse_option_pyproj(s):
    r = s.strip("+").split("=")
    if len(r) == 2:
        if "," in r[1]:
            return r[0], tuple(int(_) for _ in r[1].split(","))
        try:
            return r[0], float(r[1])
        except ValueError:
            return r[0], r[1]
    else:
        return r[0], True


class MyCRS(CRS):
    def __init__(self, proj4_params, globe=None):
        super(MyCRS, self).__init__(proj4_params, globe or Globe())


# voir https://epsg.io/2154, cliquer sur proj.4
proj4_params = (
    "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 "
    + "+y_0=6600000 +ellps=GRS80 +units=m +no_defs"
)

# Ne sert à rien si ce n'est à vérifier que le format est correct.
import pyproj

lambert93 = pyproj.Proj(proj4_params)

# Système de coordonnées de cartopy.
proj4_list = [(k, v) for k, v in map(parse_option_pyproj, proj4_params.split())]
crs_lambert93 = MyCRS(proj4_list, globe=None)

fig = plt.figure(figsize=(14, 7))
ax = fig.add_subplot(1, 2, 1, projection=ccrs.PlateCarree())
ax.set_extent([-5, 10, 42, 52])

ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=":")
ax.plot(df.long, df.lat, ".")  #
ax.set_title("France - Long/Lat")

df = df.copy()


ax = fig.add_subplot(1, 2, 2, projection=ccrs.PlateCarree())
ax.set_extent([36954, 1181938, 6133555, 7233428], crs_lambert93)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=":")
ax.plot(
    lamb_x, lamb_y, ".", transform=crs_lambert93
)  # ne pas oublier transform=crs_lambert93
ax.set_title("France - Lambert93");
../_images/c_data_enedis_cartes_15_0.png

plotly, gmaps, bingmaps#

Il faut s’authentifier. Voir gmaps, bingmaps, plotly.

geopandas#

geopandas est l’outil qui devient populaire. La partie cartes est accessible via l”API de geopandas. Il n’inclut moins de données que basemap.

[16]:

import geopandas as gpd

cities = gpd.read_file(gpd.datasets.get_path("naturalearth_cities"))
cities.head()

/tmp/ipykernel_25365/1282256058.py:2: FutureWarning: The geopandas.dataset module is deprecated and will be removed in GeoPandas 1.0. You can get the original 'naturalearth_cities' data from https://www.naturalearthdata.com/downloads/110m-cultural-vectors/.
  cities = gpd.read_file(gpd.datasets.get_path('naturalearth_cities'))

[16]:
name geometry
0 Vatican City POINT (12.45339 41.90328)
1 San Marino POINT (12.44177 43.93610)
2 Vaduz POINT (9.51667 47.13372)
3 Lobamba POINT (31.20000 -26.46667)
4 Luxembourg POINT (6.13000 49.61166) 
[17]:

world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
world.head()

/tmp/ipykernel_25365/913829029.py:1: FutureWarning: The geopandas.dataset module is deprecated and will be removed in GeoPandas 1.0. You can get the original 'naturalearth_lowres' data from https://www.naturalearthdata.com/downloads/110m-cultural-vectors/.
  world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

[17]:
pop_est continent name iso_a3 gdp_md_est geometry
0 889953.0 Oceania Fiji FJI 5496 MULTIPOLYGON (((180.00000 -16.06713, 180.00000...
1 58005463.0 Africa Tanzania TZA 63177 POLYGON ((33.90371 -0.95000, 34.07262 -1.05982...
2 603253.0 Africa W. Sahara ESH 907 POLYGON ((-8.66559 27.65643, -8.66512 27.58948...
3 37589262.0 North America Canada CAN 1736425 MULTIPOLYGON (((-122.84000 49.00000, -122.9742...
4 328239523.0 North America United States of America USA 21433226 MULTIPOLYGON (((-122.84000 49.00000, -120.0000... 
[18]:

import matplotlib.pyplot as plt

fig, ax = plt.subplots(1, 1, figsize=(6, 6))
ax.set_aspect("equal")
world.plot(ax=ax, color="white", edgecolor="black")
cities.plot(ax=ax, marker="o", color="red", markersize=5)
ax.set_title("Carte avec geopandas");
../_images/c_data_enedis_cartes_20_0.png

On restreint à l’Europe et pas trop loin de la France métropole.

[19]:

from shapely.geometry import Polygon

europe = world[world.continent == "Europe"].copy()
europe["geometry"] = europe.geometry.intersection(
    Polygon([(-10, 35), (50, 35), (50, 70), (-10, 70)])
)
europe.head()

[19]:
pop_est continent name iso_a3 gdp_md_est geometry
18 144373535.0 Europe Russia RUS 1699876 MULTIPOLYGON (((47.67591 45.64149, 46.68201 44...
21 5347896.0 Europe Norway NOR 403336 MULTIPOLYGON (((29.39955 69.15692, 28.59193 69...
43 67059887.0 Europe France FRA 2715518 MULTIPOLYGON (((6.65823 49.20196, 8.09928 49.0...
110 10285453.0 Europe Sweden SWE 530883 POLYGON ((11.46827 59.43239, 12.30037 60.11793...
111 9466856.0 Europe Belarus BLR 63080 POLYGON ((29.22951 55.91834, 29.37157 55.67009... 
[20]:

from shapely.geometry import Point

points = [Point(lon, lat) for ind, lat, lon in df[["lat", "long"]][:1000].itertuples()]
enedis = gpd.GeoDataFrame(data=dict(geometry=points))
enedis.head()

[20]:
geometry
0 POINT (2.82631 42.63236)
1 POINT (2.92725 49.00595)
2 POINT (4.21389 44.46046)
3 POINT (0.97421 47.12047)
4 POINT (5.08532 48.61706) 
[21]:

fig, ax = plt.subplots(1, 1, figsize=(5, 4))
europe.plot(ax=ax, color="white", edgecolor="black")
enedis.plot(ax=ax, marker="o", color="red", markersize=2)
ax.set_title("Carte de l'Europe avec geopandas");
../_images/c_data_enedis_cartes_24_0.png

folium#

[16]:

import folium

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

for ind, lat, lon, com in df[["lat", "long", "Nom commune"]][:50].itertuples():
    map_osm.add_child(
        folium.RegularPolygonMarker(
            location=[lat, lon], popup=com, fill_color="#132b5e", radius=5
        )
    )

map_osm

[16]:

cartopy avec les données d’OpenStreetMap#

On peut choisir également d’inclure ces détails dans une image fixe si l’image va dans un rapport écrit. On utilise les données d”OpenStreetMap avec un certain niveau de détail.

[22]:

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
from cartopy.io.img_tiles import OSM

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_extent([-5, 10, 42, 52])

imagery = OSM()
ax.add_image(imagery, 5)
# plus c'est grand, plus c'est précis, plus ça prend du temps

ax.plot(df.long[:5], df.lat[:5], ".")
ax.set_title("France");
../_images/c_data_enedis_cartes_29_0.png

Notebook on github