Les données du réseau nivo-météorologique français sont disponibles sur le site de Météo France.

On y trouve :

• un fichier json qui contient les positions des stations,
• des fichiers csv qui contiennent les mesures de ces stations, archivées par mois.

La principale bibliothèque que nous allons utiliser dans cet article est Pandas et son extension géospatiale, GeoPandas.

Pandas est une bibliothèque Python très (très) pratique pour l'analyse et la visualisation de données. Si vous ne la connaissez pas, allez vite la découvrir ! Elle permet de travailler sur des données stockées dans des dataframes.

Geopandas est son extension géospatiale : concrètement elle permet d'ajouter une colonne de géométries dans les dataframes Pandas. Elle combine de nombreux outils de l'écosystème géospatial de Python qui permettent de travailler avec les géométries : shapely, pyproj, fiona, rtree et d'autres.

Préparation des données géographiques

Le fichier des stations contient pour chacune d'entre elles : un identifiant, un nom, sa latitude, sa longitude, son altitude, une géométrie de type POINT.

Nous pouvons rapidement visualiser les positions des stations sur une carte dynamique à l'aide de la bibliothèque Folium. Cette bibliothèque permet de créer facilement des cartes générées par la bibliothèque JavaScript Leaflet.

import folium

# Create a map centered at the given latitude and longitude
stations_map = folium.Map(location=[45,1], zoom_start=5)
# Add data from a geojson file
stations_map.add_child(folium.GeoJson('postesNivo.json'))
# Display the map
display(stations_map)

Pour travailler avec ces données, nous les chargons dans un objet GeoDataFrame à l'aide de la méthode read_file de la bibliothèque Geopandas.

import geopandas as gpd   

# Read the geolocalised data
stations = gpd.read_file('postesNivo.json')
# Convert the IDs to int
stations.ID = stations.ID.astype(int)

Nous pouvons directement visualiser le tableau des données. La méthode head disponible sur tous les dataframes Pandas permet la visualisation des premières lignes, la méthode tail des dernières.

# Display the first rows (5 by default)
stations.head()

Ou leur représentation graphique avec la méthode plot qui affiche les géométries d'un geodataframe (les longitudes sont en abscisses, les latitudes en ordonnées). Cette méthode supporte l'ensemble des options qui peuvent être passées à pyplot dans la bibliothèque graphique Matplotlib.

# Plot the geometries
stations.plot()

Le système de coordonnées est accessible via l'attribut crs de notre geoDataFrame, ici nos positions sont en WGS84 (EPSG:4326). Il s'agit du système géodésique le plus fréquent lorsque l'on travaille avec des coordonnées géographiques, typiquement des positions GPS.

stations.crs

Out : {'init': 'epsg:4326'}

En plus des méthodes standard de la bibliothèque Pandas, GeoPandas fournit une indexation basée sur les coordonnées avec l'indexeur cx. Il retourne les géométries intersectant la zone de sélection. La sélection d'un extrait du dataframe n'est plus réalisée par numéro ou nom de lignes et/ou colonnes, mais directement par les coordonnées de la zone d'intérêt.

Nous utilisons cet indexeur pour extraire les stations situées dans les Alpes.

# Select a region
alpes = stations.cx[4:8,44:47]

Les stations sont représentées par des points. Pour cartographier le manteau neigeux à l'échelle de la commune, nous récupérons les emprises des communes des Alpes grâce au projet France GeoJSON qui s'appuie sur les données de l'IGN et l'INSEE. Nous créons à partir de ces fichiers un dataframe unique contenant l'emprise des communes avec lesquelles nous travaillons.

# Read geoJSON files
com73 = gpd.read_file('communes-73-savoie.geojson')
com74 = gpd.read_file('communes-74-haute-savoie.geojson')
com38 = gpd.read_file('communes-38-isere.geojson')
com05 = gpd.read_file('communes-05-hautes-alpes.geojson')
# Concatenate the geodataframes
communes = pd.concat([com73, com74, com38, com05])

Geopandas permet de réaliser des jointures spatiales avec sa méthode sjoin. Avec l'opérateur within, nous associons les stations aux communes en intersectant leurs géométries.

# Spatial join between the station's positions and the city's surfaces
map_df = gpd.sjoin(alpes, communes, how='right', op='within')
# Display the first rows of the dataframe 
map_df.head()

Préparation des données statistiques

Les hauteurs de neige mesurées par les différentes stations du réseau sont archivées dans des fichiers csv.

Nous enregistrons le contenu du fichier contenant les données du mois de décembre 2018 dans un dataframe Pandas en utilisant la méthode read_csv.

import pandas as pd

# Read the csv file containing the measurements
snow2018 = pd.read_csv('nivo.201812.csv', delimiter=';', parse_dates=['date'])
# Remove rows without a snow height measurement
snow2018 = snow2018[snow2018.ht_neige != 'mq']

Nous ne conservons que les données qui nous intéressent pour notre analyse, à savoir: l'identifiant de la station, la date de la mesure, l'altitude de la station, la hauteur de neige.

# Index the data with the dates of mesurement
snow = snow2018[['numer_sta','date','haut_sta','ht_neige']].set_index('date')

Nous regroupons ensuite les mesures disponibles par station, puis nous en calculons la moyenne mensuelle.

# Calculate the mean snow height by station
data_for_map = snow.groupby('numer_sta').resample('M').mean()

Association des données

La première étape pour associer deux dataframes est de leur attribuer le même index.

# Set the station IDs as index in the dataframes
map_df.set_index('ID', inplace=True)
data_for_map.set_index('numer_sta', inplace=True)

Nous les aggrégeons ensuite grâce à la méthode join.

# Merge the two dataframes
merged = map_df.join(data_for_map)

Comme certaines communes possèdent plusieurs stations de mesure, nous calculons une hauteur moyenne de neige pour chacune d'entre elles. Pour cela, nous réindexons notre dataframe avec le code de chaque commune (colonne qui nous sert à regrouper les mesures), puis nous calculons la hauteur moyenne de neige par commune. Cette information est enregistrée dans une nouvelle colonne.

# Reset the dataframe index
merged.reset_index(level=0, inplace=True)
# Use the city code as index
merged.set_index('code',inplace=True)
# Calculate the mean snow height for each city
merged['com_ht_neige']=merged.groupby(['code']).mean().ht_neige

Affichage sur une carte

Comme vu précédemment, une carte statique est facilement accessible avec la méthode plot de Geopandas. Et cette méthode est complétement intégrable dans les objets Matplotlib.

import matplotlib.pyplot as plt

fig, ax = plt.subplots()

communes.plot(color='gray', ax=ax)
merged.plot(column='ht_neige', cmap='Blues', linewidth=0.8, edgecolor='1.5',ax=ax, legend=True)

Pour obtenir une carte dynamique, on peut à nouveau utiliser la bibliothèque Folium et sa méthode choropleth. Cette méthode prend en paramètres un fichier geoJSON qui contient les géométries à afficher et un dataframe Pandas qui contient les valeurs statistiques associées. On enregistre donc les géométries des communes dans un fichier geoJSON avec la méthode to_file.

# Save in geojson
merged[['index', 'nom','geometry']].to_file('nivo2O18.json', driver='GeoJSON')

Désormais nous avons tout pour cartographier le manteau neigeux des communes des Alpes sur un fond de cartes Open Street Map.

# Create a map
snow_map = folium.Map(location=[45,1], zoom_start=5)

# Add the chloropleth
snow_map.choropleth(
   geo_data='nivo2O18.json', # geoJSON file
   name='choropleth',
   data=merged, # Pandas dataframe
   columns=['index','ht_neige'], # key and value of interest from the dataframe
   key_on='feature.properties.index', # key to link the json file and the dataframe
   fill_color='Blues', # colormap
   fill_opacity=0.9,
   line_opacity=0.2,
   legend_name='Hauteur moyenne de neige (m)'
)

# Add the stations
for ix, row in merged.iterrows(): 
   # Create a popup tab with the station name and its mean snow height
   popup_df = pd.DataFrame(data=[['Station', str(row['Nom'])], ['Altitude','{:4d} m'.format(int(row['haut_sta']))], ['Neige', '{:01.2f} m'.format(row['ht_neige'])]])
   popup_html = popup_df.to_html(classes='table table-striped table-hover table-condensed table-responsive', index=False, header=False)
   # Create a marker on the map
   folium.CircleMarker(location = [float(row['Latitude']),float(row['Longitude'])], radius=2, popup=folium.Popup(popup_html), color='#0000FF', fill_color='#0000FF').add_to(snow_map)

# Display the map
display(snow_map)

Conclusion

Cet article est disponible sous forme de notebook Jupyter sur le GitHub de Makina Corpus.

Vous pourrez maintenant très simplement adapter ce code pour suivre l'évolution du manteau neigeux à partir des données fournies par Météo France.

Si cet article vous a intéressé, pourquoi ne pas venir parler d'analyse géospatiale et de visualisation de données lors de l'une de nos deux formations : Python pour l'analyse géospatiale ou Python Scientifique