vignettes/biodivMapR_13.Rmd
biodivMapR_13.RmdThis tutorial illustrates an alternative approach to compute spectral
diversity with biodivMapR, This approach uses a selection
of spectral indices instead of reflectance data, as described in the standard
mode. If these spectral indices are appropriately selected in order
to ensure minimum intercorrelation, they can be used directly for the
computation of the spectral species, bypassing the preprocessing and
dimensionality reduction with PCA.
The computation of the spectral indices is based on the R package spinR, which
needs to be installed first:
devtools::install_gitlab('jbferet/spinR')
library(spinR)The data used here correspond to the data obtained from the application of this tutorial.
# libraries used in this script
library(biodivMapR)
library(stars)
library(raster)
################################################################################
## define data path and general parameterization for biodivMapR
################################################################################
# define data location
Data_Path <- '../03_RESULTS'
S2_Acq <- 'L2A_T21NYF_A018869_20201016T141050'
Input_Image_File <- file.path(Data_Path, S2_Acq, 'L2A_T21NYF_A018869_20201016T141050_Refl')
Input_Mask_File <- file.path(Data_Path, S2_Acq, 'CloudMask_Binary')
biodivMapR_Dir <- '../03_RESULTS/biodivMapR'
dir.create(path = biodivMapR_Dir,showWarnings = FALSE,recursive = TRUE)
# set parameters
Continuum_Removal <- FALSE
NDVI_Thresh <- 0.8
Blue_Thresh <- 500
NIR_Thresh <- 1500
TypePCA <- 'SI'
FilterPCA <- FALSE
window_size <- 10
nbCPU <- 4
MaxRAM <- 0.1
nbclusters <- 50The radiometric filtering is important when running biodivMapR, for the same reason as explained in the previous tutorial
Input_Mask_File <- perform_radiometric_filtering(Image_Path = Input_Image_File,
Mask_Path = Input_Mask_File,
Output_Dir = biodivMapR_Dir, TypePCA = TypePCA,
NDVI_Thresh = NDVI_Thresh, Blue_Thresh = Blue_Thresh,
NIR_Thresh = NIR_Thresh)Once radiometric thresholding is done, spectral indices are directly
computed from Sentinel-2 images with the function
ComputeSpectralIndices_Raster included in
spinR. Here, we selected four spectral indices.
# 1. read raster and corresponding spectral bands
ImBrick <- raster::brick(Input_Image_File)
HDR <- read_ENVI_header(get_HDR_name(Input_Image_File))
SensorBands <- HDR$wavelength
# 2. compute a list of spectral indices using spinR package
Sel_Indices <- c('EVI','CR_SWIR','NDVI')
Spectral_Indices <- spinR::compute_S2SI_Raster(Refl = ImBrick,
SensorBands = SensorBands,
Sel_Indices = Sel_Indices,
StackOut = T,
ReflFactor = 10000)The spectral indices obtained in the previous section are then used to refine the mask: extreme values which were not filtered out with the cloud mask and radiometric filtering may still remain. In order to reduce the risk of outliers to bias the analysis, a filtering based on interquartile range statistics is performed.
The updated shade mask combines the cloud mask produced with sen2Cor, the radiometric filter, and the outlier detection for each spectral index. It is written in the biodvMapR output directory.
The spectral indices raster stack is written in the same location as the reflectance raster.
# adjust mask to discard extreme values from spectral indices
Mask <- 0*Spectral_Indices$SpectralIndices[[1]]+1
# 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)
Mask[rast<IQRminmax[1] | rast>IQRminmax[2]] <- NA
}
# combine with mask produced from radiometric filtering
radiometricMask <- raster::raster(Input_Mask_File)
radiometricMask <- Mask*radiometricMask
radiometricMask[radiometricMask==0] <- NA
# apply mask on stack and convert as stars object
StarsObj <- st_as_stars(Spectral_Indices$SpectralIndices*radiometricMask)
# convert NA to 0 as expected in biodivMapR
Mask <- radiometricMask
Mask[is.na(Mask)] <- 0
################################################################################
## Write spectral indices and mask
################################################################################
# define path where to store spectral indices
NameStack <- paste(S2_Acq,'_SI',sep = '')
Input_SI_File <- file.path(Data_Path,S2_Acq,NameStack)
# save mask in raster file
Input_Mask_File <- paste(Input_Mask_File,'_Update',sep = '')
stars::write_stars(st_as_stars(Mask), dsn=Input_Mask_File, driver = "ENVI",type='Byte')
# save Stack of spectral indices in raster file with spectral indices defining band names
write_StarsStack(StarsObj = StarsObj, dsn = Input_SI_File, BandNames =Sel_Indices, datatype='Float32')The spectral species are then defined based on the selected components. These spectral species result from the application of a K-Means clustering of the selected components. The number of clusters may impact the resulting diversity maps. Therefore, in the perspective of operational application, the optimization of the number of clusters based on ground information is critical.
# define SpectralSpace_Output, a list expected by map_spectral_species
SpectralSpace_Output <- list('PCA_Files' = Input_SI_File,
'TypePCA' = TypePCA)
# use all spectral indices to compute diversity
SelectedPCs <- seq(1,length(Sel_Indices))
Kmeans_info <- map_spectral_species(Input_Image_File = Input_Image_File,
Input_Mask_File = Input_Mask_File,
Output_Dir = biodivMapR_Dir,
SpectralSpace_Output = SpectralSpace_Output,
nbclusters = nbclusters,
nbCPU = nbCPU, MaxRAM = MaxRAM,
SelectedPCs = SelectedPCs)Alpha diversity maps can then be computed from the spectral species distribution.
Several alpha diversity metrics are available (currently Richness, Shannon, Simpson).
print("MAP ALPHA DIVERSITY")
Index_Alpha <- c('Shannon')
map_alpha_div(Input_Image_File = Input_Image_File,
Output_Dir = biodivMapR_Dir,
TypePCA = TypePCA,
window_size = window_size,
nbCPU = nbCPU,
MaxRAM = MaxRAM,
Index_Alpha = Index_Alpha,
nbclusters = nbclusters)Beta diversity maps can then be computed from the spectral species distribution.
Beta diversity is based on a 2-steps (or 3-steps depending on image size) process:
The Bray-Curtis dissimilarity is computed between each pair of windows (defined by window_size) in the image
Then an ordination technique (default is Principal Coordinate Analysis, PCoA) is applied on the dissimilarity matrix to transform it into a lower dimension space (3-dimensions is default value), and a 3-bands raster image is produced from the coordinates of each window in this 3D space.
The explicit computation of the matrix becomes unmanageable for large amount of windows: in this case, the dissimilarity matrix is computed on a random subset of windows, followed by PCoA on this subset matrix. Then the dissimilarity between each window and this random subset is computed, and the coordinates of each window in this PCoA space is defined based on a weighted nearest neighbors rule.
print("MAP BETA DIVERSITY")
map_beta_div(Input_Image_File = Input_Image_File,
Output_Dir = biodivMapR_Dir,
TypePCA = TypePCA,
window_size = window_size,
nb_partitions = nb_partitions,
nbCPU = nbCPU,
MaxRAM = MaxRAM,
nbclusters = nbclusters)