<jupyter_start><jupyter_text>Tracer une carte en PythonLe 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](https://matplotlib.org/basemap/) (*plus trop maintenu*) ou [cartopy](http://scitools.org.uk/cartopy/) qui sont des surcouches de [matplotlib](https://matplotlib.org/). Le second moyen est une carte animée ou non dans un notebook avec des modules comme [pygal](http://www.pygal.org/en/stable/documentation/types/maps/pygal_maps_world.html), [plotly](https://plot.ly/python/choropleth-maps/). La dernière consiste à insérer des éléments sur une carte en ligne telle que [OpenStreetMap](https://www.openstreetmap.org/) et le module [folium](https://github.com/python-visualization/folium) ou [ipyleaflet](https://github.com/ellisonbg/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](https://fr.wikipedia.org/wiki/Coordonn%C3%A9es_g%C3%A9ographiques). Le problème est chaque pays a son propre système adapté à sa position géographique. Il faut souvent convertir (voir [lambert93_to_WGPS](http://www.xavierdupre.fr/app/ensae_teaching_cs/helpsphinx/ensae_teaching_cs/helpers/geo_helper.html?highlight=basemapensae_teaching_cs.helpers.geo_helper.lambert93_to_WGPS), [pyproj](https://jswhit.github.io/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 : [départements français](http://www.xavierdupre.fr/app/ensae_teaching_cs/helpsphinx/notebooks/td1a_cenonce_session_12.html?highlight=basemapcartes-avec-les-departements). La troisième difficulté est qu'on veut tracer des cartes très chargées et cela prend un temps fou.<jupyter_code>from jyquickhelper import add_notebook_menu
add_notebook_menu()
%matplotlib inline<jupyter_output><empty_output><jupyter_text>données<jupyter_code>from papierstat.datasets import load_enedis_dataset
df = load_enedis_dataset()
df.head(n=2).T<jupyter_output><empty_output><jupyter_text>cartopy[basemap](https://matplotlib.org/basemap/users/intro.htmlcartopy-new-management-and-eol-announcement) est l'ancêtre des modules de tracé de cartes sous Python mais il n'est plus vraiment maintenu. Il faut utiliser [cartopy](http://scitools.org.uk/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](http://scitools.org.uk/cartopy/docs/latest/crs/projections.html) 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](https://fr.wikipedia.org/wiki/WGS_84) ou longitude, latitude. En France, l'[INSEE](https://www.insee.fr/fr/information/2526242) utilise aussi le système [Lambert 93](https://epsg.io/2154) ou *EPSG 2154*. Source : [Introduction à la manipulation de données cartographiques](https://medium.com/@_FrancoisM/introduction-%C3%A0-la-manipulation-de-donn%C3%A9es-cartographiques-23b4e38d8f0f). Tout n'est pas parfait dans Cartopy comme ce problème [Create Cartopy projection from pyproj.Proj](https://github.com/SciTools/cartopy/issues/813).<jupyter_code>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');<jupyter_output><empty_output><jupyter_text>On peut obtenir une carte un peu plus détaillée mais le module [cartopy](https://scitools.org.uk/cartopy/docs/latest/) télécharge des données pour ce faire. Cela se faire avec l'instruction [with_scale](https://scitools.org.uk/cartopy/docs/latest/matplotlib/feature_interface.html?highlight=with_scalecartopy.feature.NaturalEarthFeature.with_scale) qui propose trois résolutions :10m, 50m, 110m.<jupyter_code>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');<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>from pyproj import Proj, transform
p1 = Proj(init='epsg:4326')  # longitude / latidude
p2 = Proj(init='epsg:2154')  # Lambert 93
transform(p1, p2, -5, 42)
transform(p1, p2, 10, 52)<jupyter_output><empty_output><jupyter_text>On convertit toutes les coordonnées.<jupyter_code>lamb_x, lamb_y = transform(p1, p2, df.long.values, df.lat.values)<jupyter_output><empty_output><jupyter_text>Et on dessine deux cartes, la première en longitude, latitude, la seconde en Lambert 93.<jupyter_code>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 +towgs84=0,0,0,0,0,0,0 +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');<jupyter_output><empty_output><jupyter_text>plotly, gmaps, bingmapsIl faut s'authentifier. Voir [gmaps](https://github.com/pbugnion/gmaps), [bingmaps](http://bharadwajyarlagadda.github.io/bingmaps/index.html), [plotly](https://plot.ly/python/maps/).  geopandas[geopandas](http://geopandas.org/) est l'outil qui devient populaire. La partie cartes est accessible via l'[API de geopandas](http://geopandas.org/mapping.html). Il n'inclut moins de données que *basemap*.<jupyter_code>import geopandas as gpd
cities = gpd.read_file(gpd.datasets.get_path('naturalearth_cities'))
cities.head()
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.head()
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');<jupyter_output><empty_output><jupyter_text>On restreint à l'Europe et pas trop loin de la France métropole.<jupyter_code>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()
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()
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");<jupyter_output><empty_output><jupyter_text>folium<jupyter_code>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<jupyter_output><empty_output><jupyter_text>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](http://openstreetmap.fr/) avec un certain [niveau de détail](https://wiki.openstreetmap.org/wiki/Zoom_levels).<jupyter_code>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');<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup