Nous avons déjà vu plusieurs fois comment manipuler et exploiter des données météorologiques avec Python et Xarray, mais pour l’instant nous avons toujours travaillé soit sur un seul point soit sur l’ensemble du jeu de données. Comment faire lorsqu’on souhaite extraire plusieurs points appartenant à une forme complexe, par exemple des limites administratives (région, pays…) ou un bassin versant ?
Dans ce tuto, nous allons voir comment “découper” un dataArray à partir d’un fichier shapefile contenant les limites des régions françaises, plus précisément nous allons :
Télécharger un historique d’ensoleillement et l’ouvrir dans un dataArray,
Télécharger un fichier shapefile contenant les limites des régions françaises, l’ouvrir avec Geopandas et extraire la forme d’une région,
Créer un masque à partir de cette forme,
Appliquer ce masque sur le datarray puis calculer l’ensoleillement moyen de la région.
Téléchargement et prétraitement des données météo
Pour commencer, téléchargez les données d’ensoleillement (Surface solar radiation downward) mensuelles d’ERA5-Land pour la France au format NetCDF. Nous avons déjà vu en détails de quoi il s’agit et comment vous pouvez y accéder via le Climate Data Store de Copernicus ou via l’API.
Une fois que vous avez téléchargé ces données, il faut les ouvrir et calculer la moyenne d’ensoleillement sur la période que l’on souhaite étudier (nous allons prendre 1981–2020). Là aussi, un autre tutoriel explique déjà comment faire mais reprenons rapidement…
D’abord, nous allons ouvrir le fichier NetCDF avec Xarray. Dans ce cas, le dataset ne contient qu’une seule variable ssrd, il faut cependant la sélectionner pour obtenir un dataarray :
import xarray as xrpath = "data.nc" #Chemin vers les donnéesDS = xr.open_dataset(path)da = DS["ssrd"]On va ensuite sélectionner la période à étudier :
from datetime import date#Dates de début et de fin de la période à étudierdebut = date(1991, 1, 1)fin = date(2020, 12, 31)da = da.sel(time = slice(debut, fin))Il ne reste plus qu’à calculer la moyenne de ce datarray sur la dimension “time” :
da = da.mean(dim = 'time')Voici le code pour l’ensemble :
import xarray as xr
from datetime import date
path = "data.nc"
debut = date(1991, 1, 1)
fin = date(2020, 12, 31)
var = "ssrd"
def open_data(path, var, debut, fin):
'''
fonction prenant :
path (str) : chemin vers la réanalyse
var (str) : nom de la variable étudiée
debut (datetime.date) : date de début de la période à étudier
fin (datetime.date) : date de fin de la période à étudier
Et renvoyant un datarray contenant la moyenne sur la période à étudier
'''
da = DS[var]
da = da.sel(time = slice(debut, fin))
da = da.mean(dim = 'time')
return da
da = open_data(path, var, debut, fin)
Il n’est pas indispensable de définir une fonction pour ça, mais cela nous assure que toutes les données intermédiaires seront supprimées à la fin des calculs puisque ce sont des variables locales… Compte-tenu de la taille des jeux de données utilisés, les oublier en mémoire pourrait vite poser problème !
Il est possible de vérifier rapidement le résultat avec un simple :
da.plot()Voilà le graphique obtenu :
On est bien en France. Les valeurs sont plus élevées dans le sud que dans le nord. Tout semble en ordre…
Téléchargement et ouverture des limites de région
Pour calculer l’ensoleillement moyen par région, il faut déjà connaitre les limites des régions… Celles-ci peuvent être téléchargées au format shapefile par exemple ici.
Vous obtenez un dossier compressé qui contient 4 fichiers portant le même nom mais avec des extensions différentes. Même si nous n’allons utiliser directement que le fichier .shp, les autres doivent être conservés dans le même dossier.
Pour ouvrir ce fichier, nous allons utiliser la librairie Geopandas, la cousine de Pandas spécialisée dans les données spatiales :
import geopandas as gpdregions = gpd.read_file("georef-france-region-millesime.shp")regionsLe résultat obtenu ressemble beaucoup à un dataframe classique, la différence réside dans le colonne “geometry” qui contient des objets géométriques (ou plus exactement des objets shapely). Dans notre cas, ces géométries sont des polygones représentant la forme des régions françaises.
On peut manipuler ce geodataframe comme on le ferait avec un dataframe. Par exemple pour filtrer les régions d’outre-mer et classer les régions par ordre alphabétique :
#Sélection des régions de métropoleregions = regions[regions.reg_area_co == 'FXX']#Classement par ordre alphabétiqueregions = regions.sort_values(by = 'reg_name_up')#Réinitialisation de l'indexregions = regions.reset_index(drop = True)Il est possible de visualiser simplement les formes sélectionnées avec :
regions.plot()Encore une fois tout semble en ordre :
Appliquer un filtre à un dataarray à partir d’un shapefile
Maintenant que nous avons toutes données nécessaires, il nous faut sélectionner les points du datarray appartenant à une région. Comment faire ?
Nous allons utiliser la méthode .where. Cette méthode s’utilise de la façon suivante :
da.where(condition)Les valeurs contenues dans le datarray sont laissées inchangées aux points où la condition est vraie. Si la condition est fausse, elles sont supprimées (c’est-à-dire remplacées par un NaN).
Par exemple, si on veut garder seulement les points du datarray où l’ensoleillement est supérieur à la moyenne, on peut le faire de la façon suivante :
#calcul de la moyennemoyenne = float(da.mean())#sélection des points où la valeur est supérieureda.where(da.values > moyenne)Vérifions avec .plot() :
Une application utile de .where est de trouver les coordonnées d’une valeur précise dans un dataarray. Par exemple si vous souhaitez savoir quel est le point du datarray avec l’ensoleillement maximal, vous pouvez le faire de la façon suivante :
#Calcul du maximum du datarray
maximum = float(da.max())
#On masque toutes les valeurs différentes du maximum
da_ = da.where(da.values == maximum)
#On supprime les latitudes où toutes les valeurs sont NaN
da_ = da_.dropna(dim = "latitude", how = 'all')
#Pareil avec les longitudes
da_ = da_.dropna(dim = 'longitude', how = 'all')
#Le dataarray ne contient plus qu'un point dont on peu extraire les coordonnées
print("Latitude du point d'ensoleillement maximal :", float(da_.latitude))
print("Longitude du point d'ensoleillement maximal :", float(da_.longitude))
Evidemment sélectionner les points se trouvant dans une région nécessite une condition plus complexe. Nous allons créer un masque, c’est-à-dire un tableau dont la valeur sera False (ou 0) pour les points se trouvant en dehors de la région et True (1) pour ceux se trouvant à l’interieur.
Commençons par récupérer la première région du geodataframe regions, il s’agit d’Auvergne-Rhône-Alpes puisque nous l’avons classé par ordre alphabétique :
region = regions.iloc[0]Si on a les coordonnées d’un point, il est facile de tester s’il se situe dans cette région, pour cela il faut d’abord créer le point correspondant avec shapely :
import shapelypoint = shapely.geometry.Point(lon, lat)Ensuite on peut simplement vérifier s’il se trouve dans la forme de la région de la façon suivante :
point.within(region.geometry)On obtient True si le point est dans la geometry, False sinon.
Il ne reste plus qu’à faire la même chose pour tous les points du dataarray. Commençons par initialiser notre masque en créant un Numpy array de mêmes dimensions que le dataarray :
import numpy as npmasque = np.zeros(da.shape)Ce masque est en fait un tableau dont chaque case correspond à un couple de latitude, longitude du dataarray : la valeur masque[m][n] correspond à la m-iène latitude dans da.latitude et à la n-ième longitude dans da.longitude. On peut savoir s’il se trouve dans la forme de la région de la façon suivante :
lat = float(da.latitude[m]) #m-ieme latitudelon = float(da.longitude[n]) #n_ieme longitude#Création du point de coordonnées lat,lonpoint = shapely.geometry.Point(lon, lat)#Teste si le point est dans la régionpoint.within(region.geometry)On modifie ensuite le masque en fonction :
masque[m][n] = point.within(region)Pour obtenir le masque complet, il suffit de parcourir l’ensemble de points d’assembler ces différentes étapes et de répéter l’opération pour l’ensemble des points. On peut le faire avec la fonction suivante :
import shapely
import numpy as np
def creation_masque(da, region):
'''
Fonction prenant :
da (datarray) : un dataarray a deux dimensions latitude et longitude
region (Pandas serie) : une ligne de geodataframe contenant la geometry à masquer
Et renvoyant :
masque (Numpy array) : masque pour le dataarray correspondant à la géometrie de la région
'''
#Initialisation du masque
masque = np.zeros(da.shape)
for m in range(masque.shape[0]):
for n in range(masque.shape[1]):
lat = float(da.latitude[m])
lon = float(da.longitude[n])
point = shapely.geometry.Point(lon, lat)
masque[m][n] = point.within(region.geometry)
return masque
Appliquer le masque et calculer la moyenne
Le plus facile reste à faire : appliquer le masque que l’on vient de créer avec .where :
da_regional = da.where(masque)Une petite vérification au cas où :
da_regional.plot()On a bien sélectionné les points du dataarray se situant en Auverge-Rhône-Alpes :
On peut maintenant calculer l’ensoleillement moyen dans la région en prenant la moyenne du dataarray :
float(da_regional.mean())On obtient 13 558 111 J/m², soit 3.766kWh/m². C’est un ordre de grandeur cohérent.
Pour terminer mettons l’ensemble du code en ordre afin de renouveler facilement le calcul pour chaque région :
import pandas as pd
#Initialisation d'un dataframe pour stocker les résultats
df = pd.DataFrame(columns = ["Région", "Ensoleillement"])
for _, region in regions.iterrows():
#Récupération du nom de la région
Nom_region = region.reg_name
#Création du masque
masque = creation_masque(da, region)
#Application du masque
da_regional = da.where(masque)
#Calcul de la moyenne
moyenne = float(da_regional.mean())
#Ajout du résultat au dataframe
df.loc[len(df)] = [Nom_region, moyenne]
#Pour un résultat plus lisible, on calcule l'ensoleillement
#par rapport à la région la plus favorisée en pourcent
df["Ensoleillement_relatif"] = (df.Ensoleillement / df.Ensoleillement.max()) * 100
df.Ensoleillement_relatif = df.Ensoleillement_relatif.astype('int')
#Puis on ordonne le dataframe en fonction de cette valeur
df = df.sort_values(by = 'Ensoleillement_relatif', ascending = False)
La région la plus ensoleillée est la Corse, Provence-Alpes-Côte d’Azur est presque ex-aequo avec juste 3% de soleil en moins, suivie plus loin par l’Occitanie (-12%). Les régions les plus défavorisées sont l’Île-de-France (qui reçoit en moyenne 28% d’énergie solaire en moins que la Corse) et les Hauts-de-France (30% moins) :
A propos : Callendar est une start-up spécialisée dans le développement de solutions innovantes pour l’évaluation des risques climatiques. Conscients du défi que représente l’adaptation au changement climatique, nous nous efforçons de partager notre expertise au travers d’outils gratuits ou de tutoriels comme celui-ci.