Le 4 mars 2022, un feu de forêt apparu dans les collines dominant la ville d’Uljin en Corée du Sud a été stoppé à quelques centaines de mètres d’une des plus grandes centrales nucléaires de la planète. Pendant l’été 2021, c’est une centrale à charbon située dans le sud de la Turquie qui avait été sauvée in-extremis d’un violent incendie.
Pour un site industriel, être situé à proximité d’une forêt représente un risque. Surtout avec le changement climatique : l’augmentation de la température donc de l’évapotranspiration et la modification des précipitations contribuent à fragiliser la végétation et à aggraver les incendies.
Mais à quel point est-il fréquent que des sites industriels dangereux soient situés en forêt ? C’est ce que nous allons déterminer dans cette nouvelle démonstration du ClimateLab.
Données sur les installations à risque
La liste des installation classées pour la protection de l’environnement est disponible via le site Georisques du Ministère de l’écologie. Elle contient la localisation de chaque installation ainsi que de nombreuses informations : nom, adresse, activité, régime de classement ICPE, classement Seveso, etc.
import geopandas as gpd
import os
#Chemin vers le fichier
path = "ICPE"
file = "InstallationsClassees.shp"
#Affichage des 10 premières lignes
gdf.head(10)
Au total, un peu plus de 100 000 sites sont classés, il y a des raffineries, des papeteries, des dépôts de produits chimiques mais aussi des entrepôts de services municipaux ou des supermarchés... Beaucoup ne nous intéressent pas, nous n’allons garder que les sites qui présentent un risque majeur, c’est-à-dire les sites Seveso seuil haut.
Par simplicité, nous allons aussi exclure les sites situés dans des départements d’outre-mer.
#Sélection des installations Seveso seuil haut
seveso = gdf[gdf.lib_seveso == 'Seveso seuil haut']
#Exclusion des départements d'outre-mer
seveso = seveso[~seveso.num_dep.isin(['971', '972', '973', '974'])]
Données d’occupation des sols
Pour connaitre le type d’environnement et de végétation présent autour de chacun de ces sites, nous allons utiliser la base de données européenne Corine land Cover.
Basée sur des observations satellites, ces données ont une résolution spatiale de 100 m. Chaque pixel est classé dans une catégorie parmi 44, par exemple “tissu urbain continu”, “landes et broussailles”, “plages, dunes et sables” et bien sûr forêts.
Ces données sont disponibles sous plusieurs format, nous utiliserons la version .tif que l’on peut ouvrir avec Xarray (comme nous l’avons vu dans un tuto précédent). Pour limiter le volume de données à garder en mémoire, nous allons commencer par sélectionner grossièrement la zone qui nous intéresse : la France métropolitaine.
import xarray as xr
#Chemin vers les données Corine
path = "u2018_clc2018_v2020_20u1_raster100m\DATA"
file = "U2018_CLC2018_V2020_20u1.tif"
#Ouverture des données
da = da.sel(band = 1, drop = True)
#Sélection de la métropole
da = da.sel(y = slice(3200000, 2000000), x = slice(3200000, 4300000))
da.plot(vmin = 0, add_colorbar = False)
L’affichage permet de vérifier que l’on ne s’est pas trompé :
Définir la zone d’intérêt autour de chaque site
Les géométries fournies avec le geodataframe des sites Seveso ne sont en réalité que des points. La zone qui nous intéresse est plus étendue : nous allons donc les remplacer par des polygones quasi-circulaires de rayon 1000 mètres centrés sur chaque site.
Le système de coordonnées du geodataframe est le WGS84. Pour pouvoir exprimer le rayon de chaque cercle en mètres nous allons convertir vers une projections azimutale locale, créer le polygone puis reconvertir en WGS84.
Par exemple pour la première ligne :
from functools import partial
import pyproj
from shapely.ops import transform
#Rayon du cercle à tracer autour du site
r = 1000
#sélection du premier site dans le géodataframe
site = seveso.iloc[0]
#Latitude et longitude du site
latitude = site.geometry.y
longitude = site.geometry.x
#Création de la projection azimutale
projection_azimutale_locale = f"+proj=aeqd +R=6371000 +units=m +lat_0={latitude} +lon_0={longitude}"
#Transformation de WGS84 vers la projection locale
wgs84_to_projlocale = partial(
pyproj.transform,
pyproj.Proj("+proj=longlat +datum=WGS84 +no_defs"),
pyproj.Proj(projection_azimutale_locale),
)
#Transformation de la projections locale vers WGS84
projlocale_to_wgs84 = partial(
pyproj.transform,
pyproj.Proj(projection_azimutale_locale),
pyproj.Proj("+proj=longlat +datum=WGS84 +no_defs"),
)
#Conversion du point vers la projection locale
centre = transform(wgs84_to_projlocale, site.geometry)
#Création d'un buffer de rayon r autour de ce point
cercle = centre.buffer(r)
#Conversion de ce polygone vers la projection originale
On peut créer rapidement une carte folium pour s’assurer que tout va bien :
import folium
#Paramètres de la carte
location = [latitude, longitude]
zoom = 14
tiles = 'cartodbpositron'
#Création de la carte
zoom_start = zoom,
tiles = tiles)
#Affichage du point
point = gpd.GeoSeries(site.geometry)
point_j.add_to(Carte)
#Affichage du polygone
polygone = gpd.GeoSeries(cercle, crs = 'WGS84')
polygone_j.add_to(Carte)
display(Carte)
Le résultat correspond bien à ce qui est souhaité :
Il ne reste donc plus qu’à faire la même chose pour chaque point et à remplacer les géométries contenues dans le geodataframe :
for i, site in seveso.iterrows():
#Latitude et longitude du site
latitude = site.geometry.y
longitude = site.geometry.x
#Création de la projection azimutale
projection_azimutale_locale = f"+proj=aeqd +R=6371000 +units=m +lat_0={latitude} +lon_0={longitude}"
#Transformation de WGS84 vers la projection locale
wgs84_to_projlocale = partial(
pyproj.transform,
pyproj.Proj("+proj=longlat +datum=WGS84 +no_defs"),
pyproj.Proj(projection_azimutale_locale),
)
#Transformation de la projections locale vers WGS84
projlocale_to_wgs84 = partial(
pyproj.transform,
pyproj.Proj(projection_azimutale_locale),
pyproj.Proj("+proj=longlat +datum=WGS84 +no_defs"),
)
#Conversion du point vers la projection locale
centre = transform(wgs84_to_projlocale, site.geometry)
#Création d'un buffer de rayon r autour de ce point
cercle = centre.buffer(r)
#Conversion de ce polygone vers la projection originale
cercle = transform(projlocale_to_wgs84, cercle)
seveso.loc[i, "geometry"] = cercle
Croisement des deux jeux de données
Pour chaque site Seveso, nous avons créé un polygone définissant la zone qui nous intéresse. L’objectif est maintenant d’extraire les données Corine dans ce polygone et de voir si on y trouve des forêts.
Comme les données Corine utilise une projection différente (ETRS89, une projection azimutale équivalente de Lambert), il faut avant tout convertir le geodataframe vers le crs correspondant :
#Transformation vers le système de coordonnées du datarray
Commençons encore une fois par travailler sur le premier site de la liste. Nous allons d’abord récupérer les coordonnées maximales et minimales du polygone afin de ne garder que la partie du datarray qui nous intéresse :
#sélection du premier site dans le géodataframe
site = seveso.iloc[0]
#Récupération la bounding box de la zone
xmin, ymin, xmax, ymax = site.geometry.bounds
#Réduction du datarray à cette zone
da_reduit = da.sel(x = slice(xmin, xmax), y = slice(ymax, ymin))
Ensuite, nous allons masquer les points du datarray restant qui se trouvent en dehors du polygone (la fonction utilisée est inspirée d’un tuto précédent) :
import numpy as np
from shapely.geometry import Point
def creation_masque(da, region):
'''
Fonction prenant :
da (datarray) : un dataarray a deux dimensions latitude et longitude
region (polygone) : la geométrie utilisée pour créer le masque
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]):
x = float(da.x[n])
y = float(da.y[m])
point = Point(x, y)
masque[m][n] = point.within(region)
return masque
#Création du masque
masque = creation_masque(da_reduit, site.geometry)
#Application du masque et affichage du résultat
da_reduit.where(masque).plot()
On obtient le résultat suivant :
Dans les données Corine, les forêts correspondent aux codes 23, 24 et 25, respectivement forêts de feuillus, de conifères et mixtes. Y en a-t-il autour de ce site ? On peut le calculer de la façon suivante :
codes_foret = [23, 24, 25]
codes = da_reduit.where(masque).data.flatten()
codes = codes[~np.isnan(codes)]
pixels_foret = codes[np.isin(codes, codes_forets)]
print(f"Pourcentage de forêt : {round(len(pixels_foret)/len(codes) * 100, 2)}%")
Pas l’ombre d’une forêt ici. Selon ces données, l’environnement dans un rayon de 1 kilomètre autour du site est composé à 56% de zone portuaire, à 18% de tissu urbain discontinu, 17% de voies d’eau, 4% de plans d’eau et 4% de terres arables non-irriguées. Cela ressemble effectivement à une description de l’avant-port de Strasbourg où se trouve ce site :
Synthèse et résultats
Il ne reste plus qu’à reproduire ce calcul pour tout les sites :
#Initialisation d'une colonne pour le résultat
seveso["prct_foret_1km"] = np.nan
for i, site in seveso.iterrows():
#Récupération la bounding box de la zone
xmin, ymin, xmax, ymax = site.geometry.bounds
#Réduction du datarray à cette zone
da_reduit = da.sel(x = slice(xmin, xmax), y = slice(ymax, ymin))
#Création du masque
masque = creation_masque(da_reduit, site.geometry)
#Récupération des données dans la zone et calcul du pourcentage de foret
codes = da_reduit.where(masque).data.flatten()
codes = codes[~np.isnan(codes)]
pixels_foret = codes[np.isin(codes, codes_forets)]
prct_foret = len(pixels_foret)/len(codes) * 100
seveso.loc[i, "prct_foret_1km"] = prct_foret
#Liste des sites situés à proximité de forêt
seveso[seveso.prct_foret_1km > 0]
Et voilà notre résultat : sur 616 sites industriels classés Seveso seuil haut, 276 sont situés à proximité d’une forêt. L’environnement est composé à 10% ou plus de forêt pour 157 sites et à 20% ou plus pour 90.
A propos : Callendar est une start-up spécialisée dans le développement de solutions innovantes pour l’évaluation des risques climatiques. Le ClimateLab permet de mobiliser rapidement nos outils et notre expertise pour répondre à une question précise.