Les cartes sont probablement les graphes les plus longs à mettre en place. Si les coordonnées longitude lattitude sont les plus utilisées, il est possible de les considérer comme des coordonnées carthésiennes sur une grande surface (sphérique). Il faut passer par une projection spécifique dans le plan. Il faut également ajouter des repères géographiques pour lire efficacement une carte comme les frontières des pays, les rivières... Tout ça prend du temps.
%matplotlib inline
from jyquickhelper import add_notebook_menu
add_notebook_menu()
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 :
Il en existe d'autres mais leur installation recèle quelques difficultés que je n'ai pas toujours eu la patience de contourner :
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");
# debug
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");
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");
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.
import os
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");
/tmp/ipykernel_31451/1518446146.py:25: UserWarning: L'image n'a pas été téléchargée. Lire le paragraphe qui précède. warnings.warn("L'image n'a pas été téléchargée. Lire le paragraphe qui précède.")
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='--');
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'avais choisi la première source de données, GADM. La licence n'est pas précisée explicitement (on pouvait 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.
L'accès aux données géographiques s'est démocratisé et nombreux pays disposent maintenant d'un site mettant à disposition gratuitement ce type de données. C'est sans doute la source de données à privilégié avec notemment la description des zones administratives comme les départements, régions dont les définitions peuvent changer.
from pyensae.datasource 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("GEOFLA_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}'\ncurrent folder: '{1}'".format(
departements, os.getcwd()))
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.shpLength, r.numRecords
(3134040, 96)
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], '_errors': {}, '_Shape__oid': 0, 'bbox': [688654.4, 6690595.300000001, 800332.3, 6811114.5]}
records[0]
Record #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.30350528758955), (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)
Record #0: ['DEPARTEM0000000000000004', '89', 'YONNE', '024', 'AUXERRE', 742447, 6744261, 748211, 6750855, '27', 'BOURGOGNE-FRANCHE-COMTE'] [0]
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.
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.