#!/usr/bin/env python # coding: utf-8 # # Données de l'IREP et Devoir¶ # # I. Création du jeu de données # In[1]: get_ipython().run_line_magic('matplotlib', 'inline') # In[2]: import pandas as pd import urllib.request import zipfile from tqdm import tqdm from pyproj import Proj, transform import cartopy.crs as ccrs import cartopy.feature as cfeature import matplotlib.pyplot as plt # ### Téléchargement des données # # On commence par télécharger toutes les données. # # En TP vous avez téléchargé manuellement le .zip correspondant à chaque année, puis vous l'avez renommé et dézippé. C'est fastidieux ! # # Dans la cellule suivante je vous montre comment le faire de manière automatisée. # In[ ]: import requests headers = { 'Connection': 'keep-alive', 'Upgrade-Insecure-Requests': '1', 'DNT': '1', 'User-Agent': 'Mozilla/5.0 (Macintosh; Intel Mac OS X 10_14_6) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/77.0.3865.90 Safari/537.36', 'Accept': 'text/html,application/xhtml+xml,application/xml;q=0.9,image/webp,image/apng,*/*;q=0.8,application/signed-exchange;v=b3', 'Referer': 'http://www.georisques.gouv.fr/dossiers/irep/telechargement', 'Accept-Encoding': 'gzip, deflate', 'Accept-Language': 'en-US,en;q=0.9,fr-FR;q=0.8,fr;q=0.7,ru;q=0.6,de;q=0.5,pt;q=0.4', } url = 'http://www.georisques.gouv.fr/irep/data/' for i in tqdm(range(2003,2018)): response = requests.get(url+str(i), headers=headers, verify=False) with open('./'+str(i)+'.zip', mode='wb') as localfile: localfile.write(response.content) with zipfile.ZipFile('./'+str(i)+'.zip',"r") as zip_ref: zip_ref.extractall("./data/"+str(i)+'/') # Explications: # # Pour savoir comment j'ai trouvé la requête, consultez [ce lien](https://www.alexkras.com/copy-any-api-call-as-curl-request-with-chrome-developer-tools/). # On convertit la commande cURL en requête Python en passant par [ce site](https://curl.trillworks.com/#python). # J'ai supprimé les cookies car ils sont inutiles. # # On enregistre ensuite le résultat de la requête (le fichier .zip), puis on dézip l'archive. # ### Système de coordonnées # # Xavier Dupré utilise un bruteforce pour trouver le système de coordonnées utilisé dans les données. Cependant ce n'est pas nécessaire. # # En cherchant un peu sur le site de l'IREP, on trouve [ceci](http://www.georisques.gouv.fr/dossiers/irep/form-etablissement/details/2534#/) qui indique que le système utilisé est "Lambert II Etendu". # # A partir de là on trouve ce [pdf](http://www.ign.fr/sites/all/files/geodesie_projections.pdf) et cette [page](https://spatialreference.org/ref/epsg/ntf-paris-lambert-zone-ii/). # # Donc dans `pyproj` on doit utiliser `epsg:27572`. # ### Déchets # # Pour une entreprise et une année données, les données répertorient l'émission de plusieurs types de déchets. On vérifie facilement que toutes les quantités de déchets sont exprimées dans la même unité. Pour une entreprise et une année donnée on **doit donc faire la somme des émissions** (sinon on se retrouve avec des cercles concentriques à la fin). En réalité il faudrait pondérer la somme par la toxicité de chaque type de déchet. # In[3]: p1 = Proj(init='epsg:4326') # longitude / latitude p2 = Proj(init='epsg:27572') # Lambert II étendu #Initialisation du jeu de données df = pd.read_csv("./data/2003/Prod_dechets_dangereux.csv") df2 = pd.read_csv("./data/2003/etablissements.csv") long, lat = transform(p2, p1, df2.Coordonnees_X.values, df2.Coordonnees_Y.values) #Passage en coordonnées GPS df2['LLX'] = long df2['LLY'] = lat df = df.merge(df2, on="Identifiant") df = df.groupby(['Identifiant']).agg( {'Quantite':'sum', 'Nom_Etablissement_x':'first', #C'est ici qu'on somme les déchets de chaque entreprise 'LLX':'first', 'LLY':'first'}).reset_index() #pour l'année df = df.rename({'Quantite': 'Quantite2003'}, axis='columns') #On renomme la colonne Quantite en Quantite2003 for i in tqdm(range(2004,2018)): #Ajouts successifs des années suivantes df_temp = pd.read_csv("./data/"+str(i)+"/Prod_dechets_dangereux.csv") df2_temp = pd.read_csv("./data/"+str(i)+"/etablissements.csv") long, lat = transform(p2, p1, df2_temp.Coordonnees_X.values, df2_temp.Coordonnees_Y.values) df2_temp['LLX'] = long df2_temp['LLY'] = lat df_temp = df_temp.merge(df2_temp, on="Identifiant") df_temp = df_temp.groupby(['Identifiant']).agg( {'Quantite':'sum', 'Nom_Etablissement_x':'first', 'LLX':'first', 'LLY':'first'}).reset_index() df_temp = df_temp.rename({'Quantite': 'Quantite'+str(i)}, axis='columns') df = df.merge(df_temp, on=["Identifiant","Nom_Etablissement_x","LLX","LLY"], how="outer") # C'est ici que se fait #l'ajout des données de l'année i dans le DataFrame final df = df.fillna('0') #Le merge outer fait apparaitre des valeurs manquantes (NaN), on les remplace par 0 # In[4]: df.sort_values(by=['Identifiant'], axis=0).head() # ### Une manière de se limiter à la métropole # drawing # ![]() # # En longitude la métropole est visiblement entre $-5$ et $10$, et en latitude elle est entre $41$ et $52$. # # Par ailleurs on trouve des bornes plus précises sur [cette page](https://spatialreference.org/ref/epsg/ntf-paris-lambert-zone-ii/) dans la section "WGS84 Bounds". # In[5]: lim_metropole = [-5, 10, 41, 52] df_metro = df[((df.LLX >= lim_metropole[0]) & (df.LLX <= lim_metropole[1]) & (df.LLY >= lim_metropole[2]) & (df.LLY <= lim_metropole[3]))] # Ca y est le jeu de données est créé, on peut passer au sujet du devoir. # # # II. Votre DM: visualisation et package # ### Installation de Cartopy # # Facile avec conda, périlleuse sinon (passez au bureau 3028 si vous voulez que je vous explique comment j'ai fait pour OSX). # ### Consignes # # Je vous demande de créer un package nommé `ensae2019` contenant la fonction `plot_geo_time_value()` (décrite ci dessous). Le dossier contenant le package doit avoir au minimum la structure suivante: # ``` # README.md # setup.py # /ensae2019 # __init__.py # main.py # /tests # tests.py # ``` # Un exemple de `setup.py` est [ici](https://pastebin.com/raw/vANMVAB3). # Le fichier `__init__.py` doit simplement contenir la ligne suivante: # `from .main import plot_geo_time_value` # La fonction `plot_geo_time_value()` doit être définie dans `main.py`. # # Je vous parlerai des test unitaires un peu plus tard (c'est la toute fin du travail). # # *** # Une fois le package créé et installé via pip, la fonction doit être importée en exécutant # `from ensae2019 import plot_geo_time_value` # # Définissez une grille avec $i$ lignes et $j$ colonnes à l'aide de la commande # `fig, axs = plt.subplots(i, j, figsize=(20,20), subplot_kw={'projection': ccrs.PlateCarree()})` # Par exemple, avec `fig, axs = plt.subplots(2, 2, figsize=(20,20), subplot_kw={'projection': ccrs.PlateCarree()})` # la grille est $2\times 2$ dans ce cas (on peut donc y faire figurer l'évolution sur 4 années). # # Etant donnés: # `x`, `y` des vecteurs de coordonnées (localisation de chaque observation), # `year` un vecteur contenant les années choisies, # `value` une matrice contenant les données à ces différentes années (avec les colonnes dans le bon ordre), # `axs` défini précédemment, # **je demande** que la commande # `plot_geo_time_value(x, y, year=years, value=values, axs=axs)` # enregistre sur le disque (dans le dossier courant) un fichier `output.pdf` contenant la grille avec les cartes. # # J'ai fait une démo concrète de ce que j'attends plus bas dans ce notebook. # # **Bonus**: vous aurez des points bonus si # 1. vous implémentez d'autres méthodes de projection, comme Mercator (la projection PlateCarree n'est pas jolie). Il faut alors ajouter un paramètre `proj='mercator'` dans votre fonction. # # ou # # 2. vous implémentez une fonction `plot_gif_geo_time_value()` qui fait la même chose mais renvoie une animation (sous format .gif, .mp4 ou .webm) # In[ ]: #Complétez cette amorce #Inspirez vous de la toute fin du notebook de Xavier Dupré: #http://www.xavierdupre.fr/app/ensae_teaching_cs/helpsphinx3/notebooks/data_irep.html#cartes def plot_geo_time_value(x, y, year, value, proj='mercator', axs=None, name=[], hue='', **kwargs): """ Visualise l'évolution temporelle d'une donnée numérique géolocalisée. :param x: longitudes (vecteur) :param y: latitudes (vecteur) :param year: années (vecteur) :param value: valeurs numériques à représenter (DataFrame ou numpy array de taille n_observations * n_years) :param proj: méthode de projection (Mercator, PlateCarree,...) (string) :param axs: axes matplotlib sur lesquels tracer (vecteur ou numpay array) :param name: noms des lieux (vecteur) :param hue: sens de la valeur numérique (:math:`CO_2`, Ammoniac, ...) :param kwargs: paramètres additionnels """ # In[ ]: #Pour la fonction bonus, complétez cette amorce def plot_gif_geo_time_value(x, y, year, value, proj='mercator', method='gif', fig=None, ax=None, name=[], hue='', **kwargs): """ Visualise l'évolution temporelle d'une donnée numérique géolocalisée avec une animation :param x: longitudes (vecteur) :param y: latitudes (vecteur) :param year: années (vecteur) :param value: valeurs numériques à représenter (DataFrame ou numpy array de taille n_observations * n_years) :param proj: méthode de projection (Mercator, PlateCarree,...) (string) :param method: type d'animation (gif, mp4 ou webm) (string) :param axs: axes matplotlib sur lesquels tracer (vecteur ou numpay array) :param name: noms des lieux (vecteur) :param hue: sens de la valeur numérique (:math:`CO_2`, Ammoniac, ...) :param kwargs: paramètres additionnels """ # # Démo # In[6]: from ensae2019 import plot_geo_time_value, plot_gif_geo_time_value # ### De 2004 à 2007, en Mercator # #### Le fichier pdf généré par la cellule ci-dessous est visible à [cette adresse](https://files.catbox.moe/0lw8f8.pdf) # In[21]: fig, axs = plt.subplots(2, 2, figsize=(20,20), subplot_kw={'projection': ccrs.Mercator()}) x, y = df_metro['LLX'], df_metro['LLY'] years = range(2004, 2008) years_str = [str(year) for year in years] values = df_metro[[colname for colname in df_metro.columns.values if colname[-4:] in years_str]].astype('float') plot_geo_time_value(x, y, year=years, value=values, proj='mercator', axs=axs) # ### De 2004 à 2007, en PlateCarree # #### Le fichier pdf généré par la cellule ci-dessous est visible à [cette adresse](https://files.catbox.moe/ja2jet.pdf) # In[22]: fig, axs = plt.subplots(2, 2, figsize=(20,15), subplot_kw={'projection': ccrs.PlateCarree()}) x, y = df_metro['LLX'], df_metro['LLY'] years = range(2004, 2008) years_str = [str(year) for year in years] values = df_metro[[colname for colname in df_metro.columns.values if colname[-4:] in years_str]].astype('float') plot_geo_time_value(x, y, year=years, value=values, proj='plate', axs=axs) # ### Gif de 2003 à 2017, en Mercator # #### Le fichier gif généré par la cellule ci-dessous est visible à [cette adresse](https://files.catbox.moe/38yndo.gif) # In[23]: fig, ax = plt.subplots(figsize=(10,10), subplot_kw={'projection': ccrs.Mercator()}) x, y = df_metro['LLX'], df_metro['LLY'] years = range(2003, 2018) years_str = [str(year) for year in years] values = df_metro[[colname for colname in df_metro.columns.values if colname[-4:] in years_str]].astype('float') plot_gif_geo_time_value(x, y, value=values, year=years, fig=fig, ax=ax, proj='mercator', method='gif') # ### Webm de 2003 à 2017, en PlateCarree # #### Le fichier webm généré par la cellule ci-dessous est visible à [cette adresse](https://files.catbox.moe/v6oxmb.webm) (clic-droit, autoloop). # In[24]: fig, ax = plt.subplots(figsize=(10,10), subplot_kw={'projection': ccrs.PlateCarree()}) x, y = df_metro['LLX'], df_metro['LLY'] years = range(2003, 2018) years_str = [str(year) for year in years] values = df_metro[[colname for colname in df_metro.columns.values if colname[-4:] in years_str]].astype('float') plot_gif_geo_time_value(x, y, value=values, year=years, fig=fig, ax=ax, proj='plate', method='webm')