Tuto : calculer l’ensoleillement moyen dans une région (ou comment masquer un datarray avec un shapefile en Python)

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 ?

extraction de données régionales d’un dataarray à partir d’une forme shapefile

Dans ce tuto, nous allons voir comment “découper” un datarray à partir d’un fichier shapefile contenant les limites des régions françaises, plus précisément nous allons :

  1. Télécharger un historique d’ensoleillement et l’ouvrir dans un dataarray,
  2. 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,
  3. Créer un masque à partir de cette forme,
  4. 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 à étudier
debut = 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 :

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 :

datarray contenant l’ensoleillement moyen en France
datarray contenant l’ensoleillement moyen en France

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")regions

Le 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étropole
regions = regions[regions.reg_area_co == 'FXX']
#Classement par ordre alphabétique
regions = regions.sort_values(by = 'reg_name_up')
#Réinitialisation de l'index
regions = 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 :

Fichier shapefile : région de France métropolitaine
Fichier shapefile : région de France métropolitaine

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 moyenne
moyenne = float(da.mean())
#sélection des points où la valeur est supérieure
da.where(da.values > moyenne)

Vérifions avec .plot() :

data array masqué avec da.where()
data array masqué avec da.where()

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 :

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 latitude
lon = float(da.longitude[n]) #n_ieme longitude
#Création du point de coordonnées lat,lon
point = shapely.geometry.Point(lon, lat)
#Teste si le point est dans la région
point.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 :

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 :

dataarray masqué avec la forme d’une région française
dataarray masqué avec la forme d’une région française
Forme de la région (à gauche) et dataarray après l’application du masque (à droite)

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 :

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) :

Ensoleillement moyen (kWh/m²) par région française
Ensoleillement moyen (kWh/m²) par région française

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.

Start-up spécialisée dans l’exploitation des données climatiques, Callendar vous aide à prendre les bonnes décisions partout où le climat actuel et futur compte

Start-up spécialisée dans l’exploitation des données climatiques, Callendar vous aide à prendre les bonnes décisions partout où le climat actuel et futur compte