Skip to article frontmatterSkip to article content

ce TP vise à vous faire découvrir quelques possibilités de manipulation et d’affichage de données géographiques

le résultat final sera de produire une carte interactive; pour des raisons techniques on ne peut pas l’afficher directement ici, mais pour avoir un aperçu, cherchez dans le zip téléchargé le fichier media/addresses-final.html et ouvrez-le dans un navigateur web

l’API de géolocalisation

pour obtenir des coordonnées latitude/longitude à partir d’une adresse, en France on peut utiliser un service public et gratuit:

# for starters, for now we only need regular pandas
# we'll go to geopandas a bit later

import pandas as pd

vous pouvez charger le fichier data/addresses.csv; toutes ces adresses sont situées à PARIS

# load the data in data/addresses.csv
# and display a few first lines

addresses = pd.read_csv('data/addresses.csv')
addresses.head(4)
Loading...

et la première chose qu’on va faire, c’est naturellement d’utiliser cette API pour géolocaliser ces adresses

un petit extrait

mais avant cela, je vous recommande de produire un fichier addresses-small.csv qui contient un petit extrait, disons les 10 ou 20 premières lignes; ce sera très utile pour débugger

# produce a small extract into data/addresses-small.csv

# your code here

une par une

c’est très pratique de pouvoir faire une recherche des adresses ‘une par une’; voici comment ça on s’y prendrait:

fort de toutes ces observations, allez chercher sur Internet la doc d’une librairie Python qui s’appelle requests
et notamment vous devriez tomber sur une page qui commence avec cet exemple

>>> r = requests.get('https://api.github.com/user', auth=('user', 'pass'))
>>> r.status_code
200
>>> r.headers['content-type']
'application/json; charset=utf8'
>>> r.encoding
'utf-8'
>>> r.text
'{"type":"User"...'
>>> r.json()
{'private_gists': 419, 'total_private_repos': 77, ...}

qui devrait vous éclairer sur l’usage de base de la librairie
(sachant que dans notre cas c’est encore plus simple puisque nous n’avons pas besoin de nous authentifier)

# requests is the swiss knife for doing http

import requests
# to localize ONE address
# we just need to forge the right URL

def localize_one(num, typ, nom):

    # we build the URL which directly contains the address pieces
    url = f"https://data.geopf.fr/geocodage/search/?q={num}+{typ}+{nom},Paris&limit=1"
    
    # print(f"localize_one is fetching page\n{url}")

    # sending request to the web server
    response = requests.get(url)

    # if all is OK, http returns a code in the [200 .. 300[ range
    if not (200 <= response.status_code < 300):
        print(f"WHOOPS.... {response.status_code=}")
        return

    # we can then read the answer 
    # remember it's a JSON string
    # so we can decode it on the fly
    return response.json()
# and here is how we could use it

localize_one(18, 'rue', 'BERNARDINS')
{'type': 'FeatureCollection', 'features': [{'type': 'Feature', 'geometry': {'type': 'Point', 'coordinates': [2.350728, 48.850345]}, 'properties': {'label': '18 Rue des Bernardins 75005 Paris', 'score': 0.7545049650349649, 'housenumber': '18', 'id': '75105_0896_00018', 'name': '18 Rue des Bernardins', 'postcode': '75005', 'citycode': '75105', 'x': 652355.3, 'y': 6861340.65, 'city': 'Paris', 'district': 'Paris 5e Arrondissement', 'context': '75, Paris, Île-de-France', 'type': 'housenumber', 'importance': 0.68417, 'street': 'Rue des Bernardins', '_type': 'address'}}], 'query': '18 rue BERNARDINS,Paris'}
# try to estimate how long it would take 
# to resolve 20_000 addresses this way

# your code here

en un seul coup

si vous avez bien lu la page qui décrit l’API, vous devez avoir remarqué qu’il y a une autre façon de lui soumettre une recherche
c’est ce qui est indiqué ici (cherchez search/csv dans la page de l’API)

curl -X POST -F data=@path/to/file.csv -F columns=voie -F columns=ville https://data.geopf.fr/geocodage/search/csv/

en Python

sauf que nous, on ne veut pas utiliser curl, on veut faire cette requête en Python; voici comment on va s’y prendre

  1. en une seule requête à l’API, on va envoyer tout le fichier csv, en lui indiquant juste quelles sont les colonnes qui contiennent les morceaux de l’adresse

  2. le résultat - toujours au format csv - pourra être également transformé en dataframe

  3. qu’il ne restera plus qu’à merge (ou join si vous préférez) avec la dataframe de départ, pour ajouter les résultats de la géolocalisation dans les données de départ pour cette étape on peut envisager de ne garder que certaines colonnes de la géolocalisation (assez bavarde par ailleurs), je vous recommande de conserver uniquement:

    • latitude, longitude - of course

    • result_city pour pouvoir vérifier la validité des résultats - ici on devrait toujours avoir Paris

    • result_type qui devrait toujours renvoyer housenumber, ça permet à nouveau de vérifier qu’on a bien une adresse connue

je vous donne ma solution; si vous essayez de la chercher par vous-même, je vous recommande d’y aller pas à pas, commencez par juste l’étape 1, puis 1 et 2, et enfin de 1 à 3
c’est utile aussi de commencer par une toute petite dataframe pour ne pas attendre des heures pendant la mise au point...

import requests
# this one is to fake a file from a string
import io

def mass_post(filename, col_number, col_type, col_name, col_city):

    # we don't care about the index for now
    # but note that if we want to merge the result
    # with the initial data, there is a need 
    # for a unique index...

    with open(filename) as feed:
        response = requests.post(
            "https://data.geopf.fr/geocodage/search/csv/", 
            files={'data': feed},
            data={'columns': [col_number, col_type, col_name, col_city]})
    if not (200 <= response.status_code < 300):
        print(f"OOPS, got {response.status_code}")
        return None

    # this time the output is csv, not json
    # note that read_csv requires a file-like object
    # and we have a str, so io.StringIO comes to the rescue
    result = pd.read_csv(io.StringIO(response.text))

    return result


def localize_many(filename, col_number, col_type, col_name, col_city):
    df = pd.read_csv(filename)
    geoloc = mass_post(filename, col_number, col_type, col_name, col_city)
    # because we don't set an index, it is safe to merge
    geoloc = geoloc[['latitude', 'longitude', 'result_city', 'result_type']]
    return df.merge(geoloc, left_index=True, right_index=True)
# try this code on the small sample for starters

addresses_small = localize_many("data/addresses-small.csv", "number", "type", "name", "city")
addresses_small
Loading...
# sanity check : make sure that all the entries have 
# result_city == 'Paris'
# and
# result_type == 'housenumber'

# your code
# when you think you're ready, go full scale
# be ready to wait for at least 40-60s 

# optional: try to record the time it takes !
# tip: see for example time.time()

# addresses = localize_many("data/addresses.csv", "number", "type", "name", "city")
# sanity check

# len(addresses)
# at this point you should store the data
# it's just good practice, as you've done one important step of the process
# plus it takes quite some time to redo

# so: store geolocalized addresses in data/addresses-geoloc.csv

# your code

afficher sur une carte

à présent qu’on a une position, on va pouvoir afficher ces adresses sur une carte

et pour ça il y a plein de libs disponibles, on va choisir folium

si nécessaire, il faut l’installer (comment on fait déjà ?)

import folium

pour commencer, allez chercher la documentation; le plus simple c’est de demander à google folium python

Question subsidiaire: comment j’ai fait d’après vous pour embarquer cette recherche dans un lien hypertexte ?

et surtout (regardez les deux premiers pour l’instant):

# ~ chatelet

CENTER = 48.856542, 2.347614

le fond de carte

# pour commencer on va recharger la dataframe précédente
import pandas as pd

addresses = pd.read_csv("data/addresses-geoloc.csv")

# et en faire un petit échantillon

addresses_small = addresses.iloc[:20]
# créez une map centrée sur ce point et de zoom 13
# n'oubliez pas de l'afficher,
# et vérifiez que vous voyez bien Paris, 
# que vous pouvez zoomer et vous déplacer, ...

# votre code
def paris_map():
    ...

# pour l'afficher
paris_map()

on ajoute les adresses

pareil mais vous ajoutez les adresses qui se trouvent dans la dataframe
éventuellement, vous pouvez comme sur l’exemple du Getting Started ajouter un tooltip avec l’adresse complète

# your code here

def map_addresses(geoloc):
    """
    creates folium map centered on Paris, with the various addresses
    contained in the input dataframe shown as a marker
    """
    pass
# and try it out

# make sure you use A SMALL DATAFRAME because with this method
# trying to display tens of thousands addresses
# is going to be SUPER SLOOOOOW !

map_addresses(addresses_small)

on sauve la carte

une fonction très sympatique de folium, c’est qu’on peut sauver cette carte sous la forme d’un fichier html, on dit standalone, c’est-à-dire autosuffisant, pas besoin de Python ni de Jupyter pour la regarder

map_small = map_addresses(addresses_small)

# je commente car à ce stade dans le notebook étudiant
# on a map_small is None et donc si on essaie de faire cela
# le notebook plante
# map_small.save("addresses-small.html")

et maintenant pour voir le résultat, ouvrez le fichier dans un autre tab du navigateur

les quartiers de Paris

maintenant on va ranger les adresses par quartier; pour cela nous avons dans le dossier data/ le fichier quartier_paris.zip qui contiennent la définition des 80 quartiers qui recouvrent Paris intra-muros

ça se présente comme ceci:

# à ce stade on a besoin de geopandas
# installez-le si besoin

import geopandas as gpd
# ça se lit facilement
quartiers = gpd.read_file("data/quartier_paris.zip", encoding="utf8")

# et le résultat est .. une dataframe
type(quartiers)
geopandas.geodataframe.GeoDataFrame
# qu'on peut afficher normalement

quartiers
Loading...
# et on peut même l'afficher sommairement avec matplotlib

# bon c'est loin d'être parfait mais ça nous suffit à ce stade
# on va voir plus loin comment incorporer ça dans la carte

quartiers.plot();
<Figure size 640x480 with 1 Axes>

pour afficher cela sur la carte, il faut se livrer à une petite gymnastique

  1. on construit la map comme on l’a vu jusqu’ici

  2. on passe à folium.GeoJson() la geo-dataframe, pour créer un objet

  3. qu’on ajoute dans la map

ce qui nous donne ceci:

def paris_map():
    """
    create a map of Paris with its 80 quartiers
    """
    map = folium.Map(
        location=CENTER,
        zoom_start=13,
        # width='80%',
    )

    # create a folium JSON object from the geo df
    folium.GeoJson(
        data=quartiers,
        # optionnally we could also tweak things on the way
        # try to uncomment these
        # style_function=lambda x: {"fillColor": "#45a012", "color": "#881212"}
        # the name is in the l_qu column
        # tooltip=folium.GeoJsonTooltip(fields=['l_qu'], aliases=['quartier']),
    # and add it to the map
    ).add_to(map)

    # and that's it
    return map
paris_map()
Loading...

un peu de couleurs

maintenant c’est à vous: il s’agit d’améliorer cela pour ajouter des couleurs aux quartiers:

  1. écrivez une fonction random_color qui renvoie une couleur au hasard, i.e. une chaine comme #12f285

  2. ajoutez dans la dataframe quartiers une colonne contenant une couleur aléatoire

  3. écrivez map_addresses_in_quartiers_2, une variante qui affiche les quartiers avec chacun une couleur

couleurs: étape 1

# in order to display a number in hexadecimal, you can use this 

x = 10
f"{x:02x}"
'0a'
# your code here

def random_color():
    """
    returns a random color as a string like e.g. #12f285
    that is to say containing 3 bytes in hexadecimal
    """
    # of course this is not the right answer
    return "#12f285"
random_color(), random_color()
('#12f285', '#12f285')

couleurs: étape 2

# add a random color column in quartiers

# your code here
quartiers.head(3)
Loading...

couleurs: étape 3

# display the quartiers with their individual color

# your code

def paris_map():
    """
    create a map of Paris with its 80 quartiers
    each quartier is shown in its individual color
    """
    pass
# display it
paris_map()

spatial join

ce qu’on aimerait bien faire à présent, c’est de trouver le quartier de chaque adresse
c’est-à-dire en pratique de rajouter dans la dataframe des adresses une ou des colonnes indiquant le quartier

et pour faire cela il existe un outil très intéressant dans geopandas qui s’appelle le spatial join et qui est décrit ici

en deux mots, l’idée c’est de faire comme un JOIN usuel (ou un pd.merge() si vous préférez)
mais pour décider si on doit aligner deux lignes (une dans la df de gauche et l’autre dans la df de droite):

ce qui signifie, en pratique, qu’on va faire ceci

  1. transformer la dataframe d’adresses en une GeoDataFrame et remplacer les colonnes latitude et longitude par une colonne position, cette fois dans un format connu de geopandas

  2. et ainsi on pourra appliquer un spatial join entre la (géo)dataframe d’adresses et la (géo)dataframe des quartiers, en choisissant ce prédicat contains

  3. et du coup modifier la carte pour afficher les adresses dans la bonne couleur - celle du quartier - pour vérifier qu’on a bien fait correctement le classement en quartiers

on doit donc obtenir autant de ligne que d’adresses, mais avec une ou des colonnes en plus (couleur, nom du quartier, ...) qui caractérisent le quartier dans lequel se trouve l’adresse

spatial join: étape 1

je vous donne le code; ce qu’il faut savoir notamment c’est qu’en geopandas il y a la notion de ‘colonne active’, celle qui contient les informations géographiques; ça n’est pas forcément hyper-intuitif la première fois...

https://geopandas.org/en/stable/docs/reference/api/geopandas.GeoDataFrame.set_geometry.html

def convert_lat_lon(df):
    """
    the input df is expected to have 2 columns named 'latitude' and 'longitude'

    this function will create a GeoDataFrame based on that, where 'latitude' and 'longitude'
    are replaced by a new column named 'position' 
    also this new column becomes the active geometry for future operations 
    (in our case the spatial join)
    """

    # we need a geopandas-friendly df
    geo_df = gpd.GeoDataFrame(df)

    # beware, here LONGITUDE comes first !
    geo_df['position'] = gpd.points_from_xy(df.longitude, df.latitude)

    # it would make sense to clean up
    # BUT
    # our map-production code relies on these columns, so...
    # del geo_df['latitude']
    # del geo_df['longitude']

    # declare the new column as the active one
    geo_df.set_geometry('position', inplace=True)

    # a technicality: turns out the shapefile was created
    # withthis coordinate reference system
    geo_df.set_crs(epsg=4326, inplace=True)

    return geo_df
# let's apply that to our small input

geoaddresses_small = convert_lat_lon(addresses_small)
# we will also need to set the active column in the quartiers (geo)dataframe

quartiers.set_geometry('geometry', inplace=True)

spatial join: étape 2

Il ne nous reste plus qu’à faire ce fameux spatial join, je vous laisse trouver le code pour faire ça

# spatial join allows to extend the addresses dataframe
# with quartier / arrondissement information

# your code

def add_quartiers(gdf):
    """
    given an addresses geo-dataframe, 
    (i.e. typically produced as an output of convert_lat_lon)
    this function will use spatial join with the quartiers information
    and return a copy of gdf extended with columns such as
    l_qu : name of the quartier
    c_ar : arrondissement number
    color: the (random) color of that quartier
    ...

    """
    # this is not the right answer...
    return gdf
# try your code

# xxx you can safely ignore this warning...
# UserWarning: CRS mismatch between the CRS of left geometries and the CRS of right geometries

geoaddresses_small_extended = add_quartiers(geoaddresses_small)
geoaddresses_small_extended.head(2)
Loading...
# verify your code

# make sure you have the right number of lines in the result

geoaddresses_small_extended.shape
(20, 9)

spatial join: étape 3

on va donc maintenant récrire map_addresses; la logique reste la même mais on veut afficher chaque adresse avec une couleur qui provient de son quartier

comme vous allez le voir, l’objet folium.Marker ne peut pas s’afficher avec une couleur arbitraire - il semble qu’il y ait seulement une petite liste de couleurs supportées

pour contourner ça, utilisez à la place un objet de type folium.CircleMarker

# rewrite map_addresses so that each address is shown 
# in the color of its quartier

def map_addresses(gdf):
    """
    (slightly) rewrite the first version of this function

    your input is now a geopandas dataframe, with the information
    about the 'quartier'

    and your job is to display all the addresses on the map, now with
    the color of the 'quartier'

    return a folium map of paris with the adresses displayed
    """
    return
# test the new function

map = map_addresses(geoaddresses_small_extended)

map

on sauve la carte

c’est sans doute un bon moment pour sauver tout ce qu’on a fait:

# votre code

on clusterise

pour pouvoir passer à l’échelle, il est indispensable de clusteriser; c’est-à-dire de grouper les points en fonction du niveau de zoom; voyons ce que ça donne

ne perdez pas de temps à chercher le code vous-même, c’est un peu hacky comme on dit, voici comment il faut faire

from folium import plugins

def map_addresses(gdf):

    # start like before
    if 'human' not in gdf.columns:
        # create a column with a human-readable address
        gdf['human'] = gdf['number'].astype(str) + ', ' + gdf['type'] + ' ' + gdf['name']

    map = paris_map()


    # we need a JavaScript function...
    # this is required to get  good performances
    callback = """
function (row) {
  // if you need to debug it
  // console.log(row)
  let circle = L.circleMarker(
      // the position
      new L.LatLng(row[0], row[1]),
      // styling
     {color: row[3], radius: 8},
    )
  // add the tooltip
  circle.bindTooltip(row[2]).openTooltip()
  return circle
}
"""

    # compute an extract with fewer columns
    # not strictly necessary but easier to deal with column names
    extract = gdf[['latitude', 'longitude', 'human', 'color']]
    cluster = plugins.FastMarkerCluster(
        extract, callback=callback,
        # nicer behaviour
        disableClusteringAtZoom=18,
    ).add_to(map)

    return map

sur l’échantillon

# let's first test it on the small extract

map_addresses(geoaddresses_small_extended)
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[44], line 3
      1 # let's first test it on the small extract
----> 3 map_addresses(geoaddresses_small_extended)

Cell In[43], line 33, in map_addresses(gdf)
     15     callback = """
     16 function (row) {
     17   // if you need to debug it
   (...)     28 }
     29 """
     31     # compute an extract with fewer columns
     32     # not strictly necessary but easier to deal with column names
---> 33     extract = gdf[['latitude', 'longitude', 'human', 'color']]
     34     cluster = plugins.FastMarkerCluster(
     35         extract, callback=callback,
     36         # nicer behaviour
     37         disableClusteringAtZoom=18,
     38     ).add_to(map)
     40     return map

File ~/.local/lib/python3.12/site-packages/geopandas/geodataframe.py:1896, in GeoDataFrame.__getitem__(self, key)
   1890 def __getitem__(self, key):
   1891     """
   1892     If the result is a column containing only 'geometry', return a
   1893     GeoSeries. If it's a DataFrame with any columns of GeometryDtype,
   1894     return a GeoDataFrame.
   1895     """
-> 1896     result = super().__getitem__(key)
   1897     # Custom logic to avoid waiting for pandas GH51895
   1898     # result is not geometry dtype for multi-indexes
   1899     if (
   1900         pd.api.types.is_scalar(key)
   1901         and key == ""
   (...)   1904         and not is_geometry_type(result)
   1905     ):

File ~/.local/lib/python3.12/site-packages/pandas/core/frame.py:4119, in DataFrame.__getitem__(self, key)
   4117     if is_iterator(key):
   4118         key = list(key)
-> 4119     indexer = self.columns._get_indexer_strict(key, "columns")[1]
   4121 # take() does not accept boolean indexers
   4122 if getattr(indexer, "dtype", None) == bool:

File ~/.local/lib/python3.12/site-packages/pandas/core/indexes/base.py:6212, in Index._get_indexer_strict(self, key, axis_name)
   6209 else:
   6210     keyarr, indexer, new_indexer = self._reindex_non_unique(keyarr)
-> 6212 self._raise_if_missing(keyarr, indexer, axis_name)
   6214 keyarr = self.take(indexer)
   6215 if isinstance(key, Index):
   6216     # GH 42790 - Preserve name from an Index

File ~/.local/lib/python3.12/site-packages/pandas/core/indexes/base.py:6264, in Index._raise_if_missing(self, key, indexer, axis_name)
   6261     raise KeyError(f"None of [{key}] are in the [{axis_name}]")
   6263 not_found = list(ensure_index(key)[missing_mask.nonzero()[0]].unique())
-> 6264 raise KeyError(f"{not_found} not in index")

KeyError: "['color'] not in index"

sur le dataset entier

# and if all goes well we can try and display the full monty

# first prepare the full dataset
geoaddresses = add_quartiers(convert_lat_lon(addresses))

# and create the map from that
final_map = map_addresses(geoaddresses)

# display the map
final_map
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[45], line 7
      4 geoaddresses = add_quartiers(convert_lat_lon(addresses))
      6 # and create the map from that
----> 7 final_map = map_addresses(geoaddresses)
      9 # display the map
     10 final_map

Cell In[43], line 33, in map_addresses(gdf)
     15     callback = """
     16 function (row) {
     17   // if you need to debug it
   (...)     28 }
     29 """
     31     # compute an extract with fewer columns
     32     # not strictly necessary but easier to deal with column names
---> 33     extract = gdf[['latitude', 'longitude', 'human', 'color']]
     34     cluster = plugins.FastMarkerCluster(
     35         extract, callback=callback,
     36         # nicer behaviour
     37         disableClusteringAtZoom=18,
     38     ).add_to(map)
     40     return map

File ~/.local/lib/python3.12/site-packages/geopandas/geodataframe.py:1896, in GeoDataFrame.__getitem__(self, key)
   1890 def __getitem__(self, key):
   1891     """
   1892     If the result is a column containing only 'geometry', return a
   1893     GeoSeries. If it's a DataFrame with any columns of GeometryDtype,
   1894     return a GeoDataFrame.
   1895     """
-> 1896     result = super().__getitem__(key)
   1897     # Custom logic to avoid waiting for pandas GH51895
   1898     # result is not geometry dtype for multi-indexes
   1899     if (
   1900         pd.api.types.is_scalar(key)
   1901         and key == ""
   (...)   1904         and not is_geometry_type(result)
   1905     ):

File ~/.local/lib/python3.12/site-packages/pandas/core/frame.py:4119, in DataFrame.__getitem__(self, key)
   4117     if is_iterator(key):
   4118         key = list(key)
-> 4119     indexer = self.columns._get_indexer_strict(key, "columns")[1]
   4121 # take() does not accept boolean indexers
   4122 if getattr(indexer, "dtype", None) == bool:

File ~/.local/lib/python3.12/site-packages/pandas/core/indexes/base.py:6212, in Index._get_indexer_strict(self, key, axis_name)
   6209 else:
   6210     keyarr, indexer, new_indexer = self._reindex_non_unique(keyarr)
-> 6212 self._raise_if_missing(keyarr, indexer, axis_name)
   6214 keyarr = self.take(indexer)
   6215 if isinstance(key, Index):
   6216     # GH 42790 - Preserve name from an Index

File ~/.local/lib/python3.12/site-packages/pandas/core/indexes/base.py:6264, in Index._raise_if_missing(self, key, indexer, axis_name)
   6261     raise KeyError(f"None of [{key}] are in the [{axis_name}]")
   6263 not_found = list(ensure_index(key)[missing_mask.nonzero()[0]].unique())
-> 6264 raise KeyError(f"{not_found} not in index")

KeyError: "['color'] not in index"
# makes sense to save the hard work

# geoaddresses.to_csv("addresses-final.csv")
# final_map.save("addresses-final.html")

références

le jeu de données utilisé ici provient à l’origine de Brando et al. (2024), légèrement retravaillé pour les besoins du TP


References
  1. Brando, C., Elgarrista, G., Lapatniova, A., Mélanie-Becquet, F., & Parmentelat, A. (2024). Vol 1898 - index adresses - Annuaire des propriétaires et des propriétés de Paris et du département de la Seine. NAKALA - https://nakala.fr (Huma-Num - CNRS). 10.34847/NKL.3038F62V