.. sectionauthor:: Daniel Courivaud .. _geo: *************** Géolocalisation *************** Parmi le volume exponentiel de données générées par l'activité humaine, un très grand nombre d'entre elles sont géolocalisées. Ceci signifie, qu'au delà de la donnée elle même, on dispose de l'information (éventuellement partielle) du lieu géographique auquel elle est attachée. On entend par information partielle, une adresse, une région administrative, un pays, etc... A cette information partielle, pour les besoins de l'affichage, on doit faire correspondre des coordonnées géographiques, qui sont le plus souvent exprimées par un couple ``(latitude, longitude)``. Pour rappel, la `latitude `_ et la `longitude `_ sont deux mesures angulaires définissant les `coordonnées sphériques `_ d'un point à la surface du globe terrestre. Géocoding ========= Le `geocoding `_ est l'action de faire correspondre une chaîne de caractère caractéristique d'un lieu géographique avec le couple ``(latitude, longitude)`` évoqué ci dessus. Nous allons travailler avec la `Base Adresse Nationale `_. Elle est co-produite par la Direction Interministérielle du Numérique (la DINSIC), par l’IGN, par La Poste, par la Direction Générale des Finances Publiques (DGFiP) et par l’association OpenStreetMap France, dans le cadre d’une convention. Elle fait partie du Service Public des Données de référence. L'API correspondante est disponible à cette adresse : `https://geo.api.gouv.fr/adresse `_. Par exemple, le code suivant permet de récupérer une chaîne de caractères contenant une structure `GeoJSON `_:: import requests, json import urllib.parse api_url = "https://api-adresse.data.gouv.fr/search/?q=" adr = "2, boulevard Blaise Pascal, 93160 Noisy le Grand" r = requests.get(api_url + urllib.parse.quote(adr)) print(r.content.decode('unicode_escape')) Plusieurs packages ou fonctions sont utilisés: - `requests `_ permet de faire des requêtes HTTP - la fonction :func:`~urllib.parse.quote` permet d'encoder l'adresse en une chaîne de caractères valide pour le protocole HTTP en remplaçant les caractères réservés par un code. En particulier, le protocole HTTP interdisant les espaces et la ponctuation dans les requêtes, l'adresse doit être transformée (url encoding) avant transmission effective sur le réseau. Ainsi, la chaîne ``"2, boulevard Blaise Pascal, 93160 Noisy le Grand"`` est remplacée par ``"2%2C%20boulevard%20Blaise%20Pascal%2C%2093160%20Noisy%20le%20Grand"`` - la fonction :func:`decode` utilisée avec le paramètre ``unicode_escape`` permet de récupérer les caractères accentués. Ainsi ``"\u00cele-de-France"`` est transformé en ``"Île-de-France"``. La requête :: https://api-adresse.data.gouv.fr/search/?q=2%2C+boulevard+Blaise+Pascal%2C+93160+Noisy+le+Grand renvoie .. code-block:: javascript :emphasize-lines: 10,11 { "type": "FeatureCollection", "version": "draft", "features": [ { "type": "Feature", "geometry": { "type": "Point", "coordinates": [ 2.583275, 48.838151 ] }, "properties": { "label": "2 Boulevard Blaise Pascal 93160 Noisy-le-Grand", "score": 0.9577360238468353, "housenumber": "2", "id": "93051_0258_00002", "type": "housenumber", "name": "2 Boulevard Blaise Pascal", "postcode": "93160", "citycode": "93051", "x": 669412.64, "y": 6859869.58, "city": "Noisy-le-Grand", "context": "93, Seine-Saint-Denis, Île-de-France", "importance": 0.5350962623151893, "street": "Boulevard Blaise Pascal" } } ], "attribution": "BAN", "licence": "ODbL 1.0", "query": "2, boulevard Blaise Pascal, 93160 Noisy le Grand", "limit": 5 } On trouvera plus d'informations sur le format GeoJSON sur la page `Wikipedia `_ et se familiariser avec le format en générant ses propres données sur `geojson.io `_. On peut vérifier la validité de l'opération de geocoding avec Google Maps. ``http://maps.google.com/?q=48.838151,2.583275``. .. image:: ../images/15-esiee-map.png :scale: 75 % :align: center Reverse geocoding ================= L'opération de `reverse geocoding `_ permet de convertir une coordonnée (latitude, longitude) en une chaîne de caractères contenant une adresse ou un lieu. Le processus est tout à fait similaire au geocoding et nécessite de forger une URL à partir d'un couple (latitude, longitude). Ainsi la requête https://api-adresse.data.gouv.fr/reverse/?lon=2.583&lat=48.838 fournit une structure JSON contenant le résultat du reverse geocoding: .. code-block:: javascript "label": "2 Boulevard Blaise Pascal 93160 Noisy-le-Grand" Automatiser le processus ------------------------ Pour automatiser le processus de (reverse) geocoding, un script Python devra: - forger l'URL à partir de l'adresse à géocoder, ou des coordonnées à "reverse géocoder" - encoder cette URL avec :func:`urllib.parse.urlencode` - utiliser le module `requests `_ pour récupérer le résultat - convertir l'objet JSON retourné en un dictionnaire Python avec la fonction :func:`json.loads` - exploiter la structure du dictionnaire obtenu pour récupérer les données intéressantes Tracer une carte ================ Python dispose de plusieurs modules permettant la représentation de données géolocalisées. Celui que nous allons utiliser dans ce cours est `Folium `_ qui permet d'utiliser l'écosystème Python pour créer des cartes interactives basées sur la bibliothèque JavaScript `Leaflet.js `_. Ces cartes sont visualisables dans un navigateur web. Cette bibliothèque est très largement utilisée: - https://whatismyipaddress.com/ - https://fr.foursquare.com - https://www.accuweather.com - https://www.openstreetmap.org - http://www.ign.com/maps/the-elder-scrolls-5-skyrim/skyrim - https://earthquake.usgs.gov/earthquakes/eventpage/usp000jkn8#map - https://www.velib-metropole.fr/map#/ Visitez la `gallerie `_ pour avoir une idée plus globale de ce que l'on peut faire. Et de comment on peut le faire... Afficher une carte ------------------ La figure ci dessous est une capture d'écran d'une carte générée avec Folium dont on trouvera la documentation complète `ici `_. .. image:: ../images/15-map1.png :scale: 50 % :align: center Le code utilisé pour produire cette figure est donné ci dessous: .. code-block:: python :linenos: import folium coords = (48.8398094,2.5840685) map = folium.Map(location=coords, tiles='OpenStreetMap', zoom_start=15) map.save(outfile='map.html') Comme on le voit, afficher une carte nécessite simplement deux actions : - instancier la classe :class:`Map` en passant au constructeur le paramètre ``location`` qui définit le centre de la carte. ``tiles`` définit la cartographie utilisée et ``zoom_start`` le niveau de zoom appliqué. - appeler la méthode :meth:`save` sur l'objet créé à l'étape précédente. La carte ainsi générée est visible dans un navigateur web. Pour plus de renseignements, `consulter la documentation de la classe Map `_. Positionner un marqueur ----------------------- Folium permet de positionner un marqueur sur un fond de carte, pour identifier un POI (Point Of Interest) particulier. .. image:: ../images/15-map2.png :scale: 50 % :align: center Le code utilisé pour produire cette figure est donné ci dessous: .. code-block:: python :linenos: :emphasize-lines: 6 import folium coords = (48.8398094,2.5840685) map = folium.Map(location=coords, tiles='OpenStreetMap', zoom_start=15) coords = [48.8490591,2.577023] folium.Marker(location=coords, popup = "ESIEE Paris").add_to(map) map.save(outfile='map.html') Ajouter un marqueur, c'est créer un objet de type :class:`Marker` dont au moins l'attribut ``location`` est renseigné. Pour plus de renseignements, `consulter la documentation de la classe Marker `_. Afficher des données sur la carte --------------------------------- Les données sont représentées sur une carte grâce à des symboles de forme, de taille et de couleur variables. Le fichier :download:`meteo2014.json <../data/meteo2014.json>` contient les températures relevées en 2014 pour l'ensemble les stations météo du territoire français. Une fois importées dans Python avec la fonction :func:`json.load`, les données sont organisées dans des dictionnaires imbriqués au format ``METEO[year][month][day][hour][code]``. Pour ce qui va suivre, on va considérer uniquement un sous ensemble de ces données: - les stations météo de France métropolitaine (ordonnées alphabétiquement) placées dans la liste STATIONS - les latitudes de ces stations placées dans la liste ordonnée LATS - les longitudes de ces stations placées dans la liste ordonnée LONGS - les températures moyennes du mois de janvier 2014, relevées à 12:00, placées dans la liste ordonnée TEMPS Les données extraites du fichier JSON sont les suivantes: .. code-block:: python STATIONS = ['ABBEVILLE', 'AJACCIO', 'ALENCON', 'BALE-MULHOUSE', 'BELLE ILE-LE TALUT', 'BORDEAUX-MERIGNAC', 'BOURGES', 'BREST-GUIPAVAS', 'CAEN-CARPIQUET', 'CAP CEPET', 'CLERMONT-FD', 'DIJON-LONGVIC', 'EMBRUN', 'GOURDON', 'LE PUY-LOUDES', 'LILLE-LESQUIN', 'LIMOGES-BELLEGARDE', 'LYON-ST EXUPERY', 'MARIGNANE', 'MILLAU', 'MONT-DE-MARSAN', 'MONTELIMAR', 'MONTPELLIER', 'NANCY-OCHEY', 'NANTES-BOUGUENAIS', 'NICE', 'ORLY', 'PERPIGNAN', "PLOUMANAC'H", 'POITIERS-BIARD', 'PTE DE CHASSIRON', 'PTE DE LA HAGUE', 'REIMS-PRUNAY', 'RENNES-ST JACQUES', 'ROUEN-BOOS', 'ST GIRONS', 'STRASBOURG-ENTZHEIM', 'TARBES-OSSUN', 'TOULOUSE-BLAGNAC', 'TOURS', 'TROYES-BARBEREY'] LATS = [50.136, 41.918, 48.4455, 47.614333, 47.294333, 44.830667, 47.059167, 48.444167, 49.18, 43.079333, 45.786833, 47.267833, 44.565667, 44.745, 45.0745, 50.57, 45.861167, 45.7265, 43.437667, 44.1185, 43.909833, 44.581167, 43.577, 48.581, 47.15, 43.648833, 48.716833, 42.737167, 48.825833, 46.593833, 46.046833, 49.725167, 49.209667, 48.068833, 49.383, 43.005333, 48.5495, 43.188, 43.621, 47.4445, 48.324667] LONGS = [1.834, 8.792667, 0.110167, 7.51, -3.218333, -0.691333, 2.359833, -4.412, -0.456167, 5.940833, 3.149333, 5.088333, 6.502333, 1.396667, 3.764, 3.0975, 1.175, 5.077833, 5.216, 3.0195, -0.500167, 4.733, 3.963167, 5.959833, -1.608833, 7.209, 2.384333, 2.872833, -3.473167, 0.314333, -1.4115, -1.939833, 4.155333, -1.734, 1.181667, 1.106833, 7.640333, 0.0, 1.378833, 0.727333, 4.02] TEMPS = [7.6, 13.5, 7.6, 6.8, 10.5, 11.5, 8.5, 9.7, 8.6, 11.8, 9.1, 7.2, 5.7, 9.2, 6.0, 7.2, 7.6, 8.4, 12.0, 6.1, 11.6, 9.6, 11.7, 6.5, 10.0, 11.7, 8.1, 12.6, 9.9, 9.1, 10.8, 9.5, 7.4, 9.0, 7.1, 10.3, 6.7, 10.8, 10.6, 8.4, 8.1] Pour tracer une carte, il faut importer le module :mod:`folium` et créer un objet :class:`Map`: .. code-block:: python import folium coords = (46.539758, 2.430331) map = folium.Map(location=coords, tiles='OpenStreetMap', zoom_start=6) .. note:: Les coordonnées utilisées ici sont celles du `centre de la France `_. Pour chaque station on crée un :class:`CircleMarker` dont on définit: - la position ; - le rayon (proportionnel à la température) ; - la couleur. Ce :class:`CircleMarker` est ajouté à la carte avec la méthode :meth:`add_to`. .. code-block:: python for i in range(len(STATIONS)): folium.CircleMarker( location = (LATS[i], LONGS[i]), radius = TEMPS[i]*2, color = 'crimson', fill = True, fill_color = 'crimson' ).add_to(map) Le résultat obtenu est conforme aux attentes: .. image:: ../images/15-map3.png :scale: 75 % :align: center Customiser l'affichage ---------------------- Pour produire quelque chose de plus visuel, on peut utiliser également la couleur pour coder la température. L'utilisation d'une ``colormap`` est appropriée. Une ``colormap`` (ou palette de couleurs) est un ensemble prédéfini de couleurs que l'on peut mapper (linéairement par exemple) avec des valeurs numériques (les données). On utilisera ici la convention classique, d'utiliser le bleu pour les températures basses, et le rouge pour les températures élevées. .. code-block:: python :emphasize-lines: 1, 6, 7, 13, 15 import folium, branca coords = (46.539758, 2.430331) map = folium.Map(location=coords, tiles='OpenStreetMap', zoom_start=6) cm = branca.colormap.LinearColormap(['blue', 'red'], vmin=min(TEMPS), vmax=max(TEMPS)) map.add_child(cm) # add this colormap on the display for lat, lng, size, color in zip(LATS, LONGS, TEMPS, TEMPS): folium.CircleMarker( location=[lat, lng], radius=size, color=cm(color), fill=True, fill_color=cm(color), fill_opacity=0.6 ).add_to(map) map.save(outfile='map.html') Le résultat obtenu montre sans surprise que les températures sont plus élevées dans le sud de la France et sur la côte Atlantique: .. image:: ../images/15-map4.png :scale: 75 % :align: center Ajout d'objets géométriques --------------------------- Il peut être parfois nécessaire de superposer à la carte des informations géographiques complémentaires. Par exemple, les contours administratifs des entités administratives (communes, cantons, départements). Le format le plus répandu pour le stockage des informations géographiques est le format `shapefile `_ (SHP) développé par l'entreprise ESRI (`La spécification complète `_). Cependant, ce format présente `un certain nombre d'inconvénients `_. Pour le développement web, le format `GeoJSON `_ est utilisé. .. tip:: Comme un grand nombre de ressources géographiques sont encore au format SHP, il pourra être utile d'utiliser le convertisseur en ligne `MapShaper `_ pour les convertir au format GeoJSON. Les données ........... `data.gouv.fr `_ fournit le découpage administratif des communes françaises au format GeoJSON. Les données sont disponibles dans le fichier :download:`datagouv-communes.geojson <../data/datagouv-communes.geojson>` (310 MB). Pour illustrer l'ajout d'objets géométriques, on va travailler avec un sous ensemble des données globales et ne conserver que les données relatives aux contours administratifs des communes d'Ile de France. L'extraction des données qui concernent l'Ile de France est effectuée à partir du fichier global en utilisant le package :mod:`geopandas`. .. tip:: Les deux méthodes ci dessous donnent de bons résultats pour l'installation de `geopandas `_ qui peut s'avérer complexe à cause des dépendances. Dans un environnement Anaconda (Windows, Mac, Linux) :: $ conda install geopandas Avec une installation Python (Windows) standard :: $ python -m pip install pipwin $ python -m pipwin install gdal $ python -m pipwin install fiona $ python -m pip install geopandas :mod:`geopandas` permet de lire, traiter et écrire les données au format ``GeoJSON`` de façon simple et efficace en reprenant le concept de data frames de :mod:`pandas`. En particulier, la fonction :func:`pandas.concat` permet d'agréger plusieurs data frames. Le code correspondant:: import geojson, geopandas, pandas # lecture du fichier global france = geopandas.read_file("datagouv-communes.geojson") l = [] # sélection des données d'Ile de France for dpt in ["75", "77", "78", "91", "92", "93", "94", "95"]: dptidf = france[france["code_commune"].str.startswith(dpt)] l.append(dptidf) # construction de la GeoDataFrame correspondante idf = pandas.concat(l) # écriture dans un fichier with open("idf.geojson", "w") as f: geojson.dump(idf, f) Après exécution du code ci dessus, le fichier :file:`idf.geojson` doit être présent dans le répertoire de travail. On peut également le trouver :download:`ici <../data/idf.geojson>` (6.5 MB). L'affichage ........... A partir des données du fichier :file:`idf.geojson`, on peut afficher le contour des communes d'Ile de France sur la carte. .. code-block:: python import folium coords = (48.7453229,2.5073644) map = folium.Map(location=coords, tiles='OpenStreetMap', zoom_start=9) #style function sf = lambda x :{'fillColor':'#E88300', 'fillOpacity':0.5, 'color':'#E84000', 'weight':1, 'opacity':1} folium.GeoJson( data="idf.geojson", name="idf", style_function= sf ).add_to(map) map.save(outfile='map.html') Le résultat obtenu: .. image:: ../images/15-map5.png :scale: 75 % :align: center .. note:: https://htmlcolorcodes.com/, https://color.adobe.com, http://www.paletton.com sont d'excellentes ressources pour le choix des couleurs définies dans la *lambda* function :func:`sf`. Cartes choroplèthes ------------------- `Une carte choroplèthe `_ est une carte pour laquelle des surfaces élémentaires (définies par un ou plusieurs fichiers GeoJSON) sont colorées conformément à une donnée (une statistique par exemple). On peut par exemple produire une carte de la population par commune en Ile de France en croisant le contour administratif des communes avec les données de recensement de l'`INSEE `_. Ces données de recensement sont disponibles dans le fichier :download:`insee-pop-communes.csv <../data/insee-pop-communes.csv>` (1 MB). .. image:: ../images/15-map6.png :scale: 75 % :align: center Tracer une telle carte requiert plusieurs étapes: - préparation des données géographiques - préparation des données numériques - création d'une instance de :class:`Folium.Map` - utilisation de la classe :class:`Choropleth` à cette instance Préparation des données géographiques ..................................... L'observation de la structure du fichier GeoJSON montre que: - les données sont contenues dans une liste de dictionnaires, un par commune - pour chaque chaque dictionnaire - la clé ``properties`` contient le code de la commune. Ce code servira à faire le lien avec le dataset de l'INSEE - la clé ``geometry`` contient les données géographiques de la commune sous la forme d'un polygone dont les sommets sont définis par les coordonnées GPS. .. code-block:: javascript { "type":"FeatureCollection", "features": [ { "type": "Feature", "geometry": { "type": "Polygon", "coordinates": [[ [2.304785272419358, 48.67849158613953], ... [2.304785272419358, 48.67849158613953] ]] }, "properties": { "code_epci": "200056232", "commune": "Longjumeau", ... "code_commune": "91345", ... "intitule_commune": "91345 Longjumeau"} }, }, ... ] } Préparation des données numériques .................................. Les données doivent être lues et stockées dans une data frame `Pandas `_:: >>> pop_data = pd.read_csv('insee-pop-communes.csv', sep=';') >>> pop_data DEPCOM COM PMUN PCAP PTOT 0 01001 L' Abergement-Clémenciat 776 18 794 1 01002 L' Abergement-de-Varey 248 1 249 2 01004 Ambérieu-en-Bugey 14035 393 14428 3 01005 Ambérieux-en-Dombes 1689 34 1723 4 01006 Ambléon 111 6 117 ... ... ... ... ... ... 34990 97419 Sainte-Rose 6418 79 6497 34991 97420 Sainte-Suzanne 23505 199 23704 34992 97421 Salazie 7312 75 7387 34993 97422 Le Tampon 78629 1076 79705 34994 97423 Les Trois-Bassins 7139 95 7234 [34995 rows x 5 columns] La colonne ``DEPCOM`` contient le code INSEE des communes:: >>> dpt_code = pop_data['DEPCOM'] >>> dpt_code 0 01001 1 01002 2 01004 3 01005 4 01006 ... 34990 97419 34991 97420 34992 97421 34993 97422 34994 97423 Name: DEPCOM, Length: 34995, dtype: object La création d'un masque booléen permet de construire le sous ensemble recherché:: >>> mask = ( ( dpt_code.str.startswith('75') ) ... | ( dpt_code.str.startswith('77') ) ... | ( dpt_code.str.startswith('78') ) ... | ( dpt_code.str.startswith('91') ) ... | ( dpt_code.str.startswith('92') ) ... | ( dpt_code.str.startswith('93') ) ... | ( dpt_code.str.startswith('94') ) ... | ( dpt_code.str.startswith('95') ) ) >>> mask 0 False 1 False 2 False 3 False 4 False ... 34990 False 34991 False 34992 False 34993 False 34994 False Name: DEPCOM, Length: 34995, dtype: bool Il est utilisé pour effectuer le filtrage par ligne:: >>> pop_data = pop_data[mask] >>> pop_data DEPCOM COM PMUN PCAP PTOT 29297 75101 Paris 1er Arrondissement 16266 129 16395 29298 75102 Paris 2e Arrondissement 20900 142 21042 29299 75103 Paris 3e Arrondissement 34115 274 34389 29300 75104 Paris 4e Arrondissement 28088 282 28370 29301 75105 Paris 5e Arrondissement 58850 781 59631 ... ... ... ... ... ... 34878 95676 Villers-en-Arthies 503 10 513 34879 95678 Villiers-Adam 859 11 870 34880 95680 Villiers-le-Bel 27676 132 27808 34881 95682 Villiers-le-Sec 185 3 188 34882 95690 Wy-dit-Joli-Village 333 7 340 [1287 rows x 5 columns] La carte proprement dite ........................ Après création de la carte, il faut utiliser la classe :class:`Choropleth` dont les paramètres importants sont: - :data:`geo_data` : les données géographiques au format GeoJSON ; - :data:`pop_data` : les données numériques au format `Pandas `_ ; - :data:`columns` : la paire clé/valeur des données numériques ; - :data:`key_on` : nom de la variable géographique permettant de faire le lien entre données géographiques et numériques. .. code-block:: python coords = (48.7190835,2.4609723) map = folium.Map(location=coords, tiles='OpenStreetMap', zoom_start=9) folium.Choropleth( geo_data=geo_data, # geographical data name='choropleth', data=pop_data, # numerical data columns=['DEPCOM', 'PTOT'], # numerical data key/value pair key_on='feature.properties.code_commune', # geographical property used to establish correspondance with numerical data fill_color='YlGn', fill_opacity=0.7, line_opacity=0.2, legend_name='Population' ).add_to(map) map.save(outfile='map.html') .. note:: La ressource `Python Folium: Create Web Maps From Your Data `_ est un bon complément à ce qui précède. .. warning:: Le fichier GeoJSON décrit l'agglomération parisienne comme une seule et même entité alors que les données numériques de population découpent la capitale en arrondissements. Lorsque la correspondance ne peut être établie entre les données géographiques et les données numériques, aucune donnée n'est affichée pour le polygone en question. Deux solutions sont possibles : - on insère les données géographiques des arrondissements dans le fichier GeoJSON ; - ou on aggrège les données de population des arrondissements pour constituer une population unique pour Paris dans les données numériques. Exercice ======== Il n'y a pas de fichier squelette associé à cet exercice. La population totale de la commune est une donnée importante mais elle ne reflète pas complètement l’urbanisation du territoire. La densité de population serait une grandeur plus intéressante dans ce cas. La surface de chacune des communes est une donnée qui n'apparait pas dans les datasets utilisés dans le cours. On pourrait imaginer la calculer à partir des coordonnées du polygone correspondant. Mais cette opération est complexe (voir `l'algorithme du lacet `_). On peut également la rechercher dans une source de données différente. Utiliser ce `dataset `_ pour produire une carte choroplèthe de la densité de population des communes d'Ile de France.