vignettes/biodivMapR_19.Rmd
biodivMapR_19.RmdCe tutoriel accompagne la suite de la séquence 6 de la section Cartographie de la biodiversité par imagerie satellite du module de formation COPERNICUS et biodiversité de la Copernicus Academy.
Il permet d’appliquer biodivMapR à des données issues d’imagerie Sentinel-2 (précédemment téléchargées et pré-traitées) et définies par l’utilisateur. Dans le cadre de cette séquence, des indices spectraux sont d’abord calculés à partir des données Sentinel-2, pour être ensuite utilises avec biodivMapR.
La principale différence avec le tutoriel precedent est le calcul préalable d’indices spectraux à partir des données de réflectance Sentinel-2.
Le package spinR est un package permettant de calculer des indices spectraux à partir de données d’imagerie au format raster.
biodivMapR
Les variables définies lors du précédent tutoriel sont rappelées ici.
# load biodivMapR & spinR
library(biodivMapR)
library(spinR)
# load additional packages
library(raster)
library(stars)
Input_Image_File <- Refl_path
Input_Mask_File <- cloudmasks$BinaryMask
Output_Dir <- file.path(result_path,'biodivMapR')
# TypePCA sets the name for the directory name where to save results.
# Here, as no spectral transformation (PCA, SPCA, MNF) is applied, user can define any other name
TypePCA <- 'SI'
# number of clusters (spectral species)
nbclusters <- 50
# window size: RS diversity will be computed for this window size over the raster
window_size <- 10
# computational parameters
nbCPU <- 4
MaxRAM <- 0.2Cette étape permet de réaliser un filtre pour eliminer les pixels ennuagés, ombragés ou non végétalisés. C’est un filtre assez grossier car basé sur des règles de décision très simplifiées:
le seuillage du NDVI permet d’éliminer les pixels non végétalisés.
le seuillage dans le domaine visible correspondant au bleu permet d’éliminer les pixels présentant une réflectance particulièrement forte pour de la végétation, laissant supposer la présence de brumes ou de nuages non identifiés par le masque nuage produit par la méthode de correction atmosphérique.
le seuillage dans le domaine proche infrarouge permet d’éliminer les pixels présentant une réflectance particulièrement faible pour de la végétation, laissant supposer que le pixel est situé à l’ombre.
A l’issue de cette étape, le masque défini par
Input_Mask_File est mis à jour.
# Define levels for radiometric filtering
NDVI_Thresh <- 0.8
Blue_Thresh <- 500
NIR_Thresh <- 1500
Input_Mask_File <- biodivMapR::perform_radiometric_filtering(Image_Path = Input_Image_File,
Mask_Path = Input_Mask_File,
Output_Dir = Output_Dir,
TypePCA = TypePCA,
NDVI_Thresh = NDVI_Thresh,
Blue_Thresh = Blue_Thresh,
NIR_Thresh = NIR_Thresh)Cette partie fait appel au package spinR pour produire les indices spectraux. Un crtain nombre d’indices spectraux est déja implementé. L’utilisateur peut aussi utiliser un interpreteur permettant de calculer un indice qui ne serait pas implementé.
Ici, nous sélectionnons trois indices spectraux:
Le CCCI est supposé être sensible à la teneur en chlorophylle de la végétation. Le CR_SWIR est supposé être sensible à la teneur en eau de la végétation. Le LAI est supposé être sensible au Leaf Area Index.
La robustesse de ces relations entre indices spectraux et propriétés de végétation n’est pas l’objet de ce tutoriel. L’objectif ici est de permettre à l’utilisateur de sélectionner des variables préalablement dérivées de données de télédetection, jugées pertinentes.
# get spectral bands corresponding to S2 image
HDR <- read_ENVI_header(get_HDR_name(Input_Image_File))
SensorBands <- HDR$wavelength
# get raster stack from Input_Image_File
Refl <- raster::stack(Input_Image_File)
# Define a list of spectral indices
Sel_Indices <- c('CCCI','CR_SWIR','LAI_SAVI')
# then compute the spectral indices from raster data, using spectral bands as
# close as possible from Sentinel-2 spectral bands used for the spectral indices
Spectral_Indices <- spinR::compute_S2SI_Raster(Refl = Refl,
SensorBands = SensorBands,
Sel_Indices = Sel_Indices,
StackOut = T,
ReflFactor = 10000)Une fois les indices spectraux calculés, un masque peut être mis a jour afin d’éliminer les pixels correspondant à des outliers. De plus, les indices spectraux doivent être écrits sur le disque sous la forme d’un stack correspondant au meme format que le stack produit lorsqu’une transformation spectrale est appliquée.
# initial mask values
Mask <- raster::raster(Input_Mask_File)
raster::values(Mask)[raster::values(Mask)==0] <- NA
# remove outliers from spectral indices
for (idx in Sel_Indices){
rast <- Spectral_Indices$SpectralIndices[[idx]]
IQRminmax <- biodivMapR::IQR_outliers(DistVal = raster::values(rast),weightIRQ = 3)
raster::values(Mask)[raster::values(rast)<IQRminmax[1] | raster::values(rast)>IQRminmax[2]] <- NA
}
# apply mask on stack
for (i in 1:dim(Spectral_Indices$SpectralIndices)[3]){
Spectral_Indices$SpectralIndices[[i]] <- Spectral_Indices$SpectralIndices[[i]]*raster::values(Mask)
}
# convert into stars object
StarsObj <- stars::st_as_stars(Spectral_Indices$SpectralIndices)
# write SI stack
SI_Path <- file.path(rootdir,'SpectralIndices')
dir.create(SI_Path,showWarnings = F, recursive = T)
SI_file <- file.path(SI_Path,'SI_Stack')
Input_Mask_File_update <- file.path(SI_Path,'SI_Mask')
# save mask in raster file
stars::write_stars(st_as_stars(Mask),
dsn = Input_Mask_File_update,
driver = "ENVI",
type = 'Byte')
# save Stack of spectral indices in raster file with spectral indices defining band names
biodivMapR::write_StarsStack(StarsObj = StarsObj,
dsn = SI_file,
BandNames = Sel_Indices,
datatype='Float32')Une fois les indices spectraux sélectionnés, la carte d’espèces spectrales peut être calculée. Ces espèces spectrales sont définies à l’aide d’un algorithme de k-means clustering.
SpectralSpace_Output <- list('PCA_Files' = SI_file,
'TypePCA' = TypePCA)
Kmeans_info <- biodivMapR::map_spectral_species(Input_Image_File = Input_Image_File,
Input_Mask_File = Input_Mask_File_update,
Output_Dir = Output_Dir,
SpectralSpace_Output = SpectralSpace_Output,
nbclusters = nbclusters,
SelectedPCs = seq(1,dim(raster::stack(SI_file))[3]),
nbCPU = nbCPU, MaxRAM = MaxRAM)Une carte de diversité \(\alpha\) (ici correspondant à l’indice de Shannon) peut alors être calculée.
Index_Alpha <- c('Shannon')
map_alpha_div(Input_Image_File = Input_Image_File,
Input_Mask_File = Input_Mask_File_update,
Output_Dir = Output_Dir,
TypePCA = TypePCA,
window_size = window_size,
nbCPU = nbCPU,
MaxRAM = MaxRAM,
Index_Alpha = Index_Alpha,
nbclusters = nbclusters)De la même manière, une carte de diversité \(\beta\) est aussi calculée.
map_beta_div(Input_Image_File = Input_Image_File,
Output_Dir = Output_Dir,
TypePCA = TypePCA,
window_size = window_size,
nbCPU = nbCPU,
MaxRAM = MaxRAM,
nbclusters = nbclusters)