update v1.10.0

Merge branch 'dev'

# Conflicts:
#	R/Lib_MapBetaDiversity.R

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: biodivMapR
 Title: biodivMapR: an R package for a- and ß-diversity mapping using remotely-sensed images
-Version: 1.9.11
+Version: 1.10.0
 Authors@R: c(person(given = "Jean-Baptiste", 
 					family = "Feret", 
 					email = "jb.feret@teledetection.fr", 

NAMESPACE

@@ -1,12 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
-export(Compute_ALPHA_SSD_per_window)
-export(Compute_ALPHA_SSD_per_window_list)
-export(Compute_ALPHA_per_window)
-export(ContinuumRemoval)
 export(ENVI_type2bytes)
 export(IQR_outliers)
-export(Normalize_SSD)
 export(VectorInRasterFootprint)
 export(WeightedCoordsNN)
 export(Write_Big_Image)
@@ -19,6 +14,9 @@ export(center_reduce)
 export(change_resolution_HDR)
 export(check_data)
 export(compute_ALPHA_FromPlot)
+export(compute_ALPHA_SSD_per_window)
+export(compute_ALPHA_SSD_per_window_list)
+export(compute_ALPHA_per_window)
 export(compute_BCdiss)
 export(compute_BETA_FromPlots)
 export(compute_FUNCT)
@@ -28,6 +26,7 @@ export(compute_alpha_metrics)
 export(compute_beta_metrics)
 export(compute_spectral_species)
 export(compute_spectral_species_FieldPlots)
+export(continuumRemoval)
 export(coordPix_kernel)
 export(create_hdr)
 export(define_output_directory)
@@ -69,6 +68,7 @@ export(mean_filter)
 export(minmax)
 export(mnf)
 export(noise)
+export(normalize_SSD)
 export(ordination_to_NN)
 export(ordination_to_NN_list)
 export(pca)

NEWS.md

@@ -1,6 +1,16 @@
+# biodivMapR v1.10.0
+
+## Fix
+- now uses red band instead of red edge band for the computation of NDVI
+--> changes thresholding
+
+## changes
+- simplified inputs for functions such as map_spectral_species(), init_kmeans(), map_beta_div()
+- updated vignettes
+
 # biodivMapR v1.9.11
 
-## Addition
+## changes
 - optimized codes for the computation of spectral species, alpha and beta diversity
 - addition of progress bars during the different steps of the computation of spectral species, alpha and beta diversity maps
 - addition of functions Compute_ALPHA_SSD_per_window_list, Compute_ALPHA_SSD_per_window, Compute_ALPHA_per_window, prepare_HDR_SSD, prepare_HDR_Sunlit, RW_bytes_all, RW_bytes

R/Lib_ContinuumRemoval.R

@@ -43,13 +43,13 @@ apply_continuum_removal <- function(Spectral_Data, Spectral, nbCPU = 1) {
       with_progress({
         p <- progressr::progressor(steps = nb_CR)
         Spectral_Data_tmp <- future.apply::future_lapply(Spectral_Data,
-                                                         FUN = ContinuumRemoval,
+                                                         FUN = continuumRemoval,
                                                          Spectral_Bands = Spectral$Wavelength,
                                                          p = p)
       })
       future::plan(sequential)
     } else {
-      Spectral_Data_tmp <- lapply(Spectral_Data, FUN = ContinuumRemoval,
+      Spectral_Data_tmp <- lapply(Spectral_Data, FUN = continuumRemoval,
                                   Spectral_Bands = Spectral$Wavelength)
     }
     Spectral_Data <- do.call("rbind", Spectral_Data_tmp)
@@ -75,7 +75,7 @@ apply_continuum_removal <- function(Spectral_Data, Spectral, nbCPU = 1) {
 #' @return samples from image and updated number of pixels to sampel if necessary
 #' @export
 
-ContinuumRemoval <- function(Minit, Spectral_Bands, p = NULL) {
+continuumRemoval <- function(Minit, Spectral_Bands, p = NULL) {
 
   # Filter and prepare data prior to continuum removal
   CR_data <- filter_prior_CR(Minit, Spectral_Bands)
@@ -159,6 +159,7 @@ ContinuumRemoval <- function(Minit, Spectral_Bands, p = NULL) {
   } else {
     CR_Results <- matrix(0, ncol = (nbBands - 3), nrow = nbSamples)
   }
+  if (!is.null(p)){p()}
   list <- ls()
   rm(list = list[-which(list == "CR_Results")])
   rm(list)

R/Lib_FilterData.R

@@ -25,9 +25,10 @@
 #'
 #' @return MaskPath = updated mask file
 #' @export
-perform_radiometric_filtering <- function(Image_Path, Mask_Path, Output_Dir, TypePCA = "SPCA",
-                                          NDVI_Thresh = 0.5, Blue_Thresh = 500, NIR_Thresh = 1500,
-                                          Blue = 480, Red = 700, NIR = 835) {
+perform_radiometric_filtering <- function(Image_Path, Mask_Path = FALSE, Output_Dir,
+                                          TypePCA = "SPCA",
+                                          NDVI_Thresh = 0.8, Blue_Thresh = 500, NIR_Thresh = 1500,
+                                          Blue = 480, Red = 670, NIR = 835) {
   # check if format of raster data is as expected
   check_data(Image_Path)
   if (!Mask_Path==FALSE){
@@ -42,7 +43,8 @@ perform_radiometric_filtering <- function(Image_Path, Mask_Path, Output_Dir, Typ
     print("Update mask based on NDVI, NIR and Blue threshold")
   }
   Shade_Update <- file.path(Output_Dir_Full, "ShadeMask_Update")
-  Mask_Path <- create_mask_from_threshold(ImPath = Image_Path, MaskPath = Mask_Path, MaskPath_Update = Shade_Update,
+  Mask_Path <- create_mask_from_threshold(ImPath = Image_Path,
+                                          MaskPath = Mask_Path, MaskPath_Update = Shade_Update,
                                           NDVI_Thresh = NDVI_Thresh, Blue_Thresh = Blue_Thresh, NIR_Thresh = NIR_Thresh,
                                           Blue = Blue, Red = Red, NIR = NIR)
   return(Mask_Path)
@@ -65,7 +67,8 @@ perform_radiometric_filtering <- function(Image_Path, Mask_Path, Output_Dir, Typ
 #' @param NIR numeric. spectral band corresponding to the NIR channel (in nanometers)
 #
 # @return MaskPath path for the updated shademask produced
-create_mask_from_threshold <- function(ImPath, MaskPath, MaskPath_Update, NDVI_Thresh, Blue_Thresh, NIR_Thresh,
+create_mask_from_threshold <- function(ImPath, MaskPath, MaskPath_Update,
+                                       NDVI_Thresh, Blue_Thresh, NIR_Thresh,
                                        Blue = 480, Red = 690, NIR = 835) {
   # define wavelength corresponding to the spectral domains Blue, Red and NIR
   Spectral_Bands <- c(Blue, Red, NIR)

R/Lib_MapAlphaDiversity.R

@@ -415,7 +415,7 @@ compute_SSD <- function(Image_Chunk, window_size, nbclusters,
   }
 
   # tictoc::tic()
-  # alphaSSD <- pbapply::pblapply(listAlpha, Compute_ALPHA_SSD_per_window, MinSun, nb_partitions, alphaIdx)
+  # alphaSSD <- pbapply::pblapply(listAlpha, compute_ALPHA_SSD_per_window, MinSun, nb_partitions, alphaIdx)
   # tictoc::toc()
 
   if (nbCPU > 1){
@@ -428,15 +428,15 @@ compute_SSD <- function(Image_Chunk, window_size, nbclusters,
       with_progress({
         p <- progressr::progressor(steps = nbCPU)
         alphaSSD <- future.apply::future_lapply(X = listAlpha,
-                                                FUN = Compute_ALPHA_SSD_per_window_list,
+                                                FUN = compute_ALPHA_SSD_per_window_list,
                                                 nb_partitions = nb_partitions,
                                                 nbclusters = nbclusters,
                                                 alphaIdx = alphaIdx,
                                                 MinSun = MinSun, pcelim = pcelim, p = p)
       })
     } else {
       alphaSSD <- future.apply::future_lapply(X = listAlpha,
-                                              FUN = Compute_ALPHA_SSD_per_window_list,
+                                              FUN = compute_ALPHA_SSD_per_window_list,
                                               nb_partitions = nb_partitions,
                                               nbclusters = nbclusters,
                                               alphaIdx = alphaIdx,
@@ -454,7 +454,7 @@ compute_SSD <- function(Image_Chunk, window_size, nbclusters,
   } else {
 
     alphaSSD <- lapply(X = listAlpha,
-                       FUN = Compute_ALPHA_SSD_per_window,
+                       FUN = compute_ALPHA_SSD_per_window,
                        nb_partitions = nb_partitions, nbclusters = nbclusters,
                        alphaIdx = alphaIdx, MinSun = MinSun, pcelim = pcelim)
 
@@ -595,11 +595,11 @@ get_Simpson <- function(Distrib) {
 #' corresponding to the list of windows
 #' @export
 
-Compute_ALPHA_SSD_per_window_list <- function(listAlpha, nb_partitions, nbclusters,
+compute_ALPHA_SSD_per_window_list <- function(listAlpha, nb_partitions, nbclusters,
                                               alphaIdx, MinSun = 0.25, pcelim = 0.02,
                                               p = NULL) {
   alphaSSD <- lapply(X = listAlpha,
-                     FUN = Compute_ALPHA_SSD_per_window,
+                     FUN = compute_ALPHA_SSD_per_window,
                      nb_partitions = nb_partitions, nbclusters = nbclusters,
                      alphaIdx = alphaIdx, MinSun = MinSun, pcelim = pcelim)
 
@@ -628,7 +628,7 @@ Compute_ALPHA_SSD_per_window_list <- function(listAlpha, nb_partitions, nbcluste
 #' corresponding to the window
 #' @export
 
-Compute_ALPHA_SSD_per_window <- function(listAlpha, nb_partitions, nbclusters,
+compute_ALPHA_SSD_per_window <- function(listAlpha, nb_partitions, nbclusters,
                                          alphaIdx, MinSun = 0.25, pcelim = 0.02) {
 
   if (listAlpha$PCsun > MinSun) {
@@ -639,7 +639,7 @@ Compute_ALPHA_SSD_per_window <- function(listAlpha, nb_partitions, nbclusters,
     if (!typeof(SSD)=='list'){
       SSD <- lapply(X = snow::splitRows(listAlpha$data, ncl = nrow(listAlpha$data)), FUN = table)
     }
-    alphawin <- lapply(SSD, FUN = Compute_ALPHA_per_window,
+    alphawin <- lapply(SSD, FUN = compute_ALPHA_per_window,
                        nbPix_Sunlit = listAlpha$nbPix_Sunlit,
                        alphaIdx = alphaIdx,
                        nbclusters = nbclusters,
@@ -672,7 +672,7 @@ Compute_ALPHA_SSD_per_window <- function(listAlpha, nb_partitions, nbclusters,
 #' @return shannon, Simpson, Fisher, SSDMap
 #' @export
 
-Compute_ALPHA_per_window <- function(SSD, nbPix_Sunlit, alphaIdx, nbclusters,
+compute_ALPHA_per_window <- function(SSD, nbPix_Sunlit, alphaIdx, nbclusters,
                                      pcelim = 0.02){
 
   ClusterID <- as.numeric(names(SSD))

R/Lib_MapBetaDiversity.R

@@ -16,7 +16,6 @@
 #' @param Output_Dir character. Output directory.
 #' @param window_size numeric. Dimensions of the spatial unit.
 #' @param TypePCA character. Type of PCA (PCA, SPCA, NLPCA...).
-#' @param nb_partitions numeric. Number of partitions (repetitions) to be computed then averaged.
 #' @param nbclusters numeric. Number of clusters defined in k-Means.
 #' @param Nb_Units_Ordin numeric. Maximum number of spatial units to be processed in NMDS.
 #' --> 1000 will be fast but may not capture important patterns if large area
@@ -39,7 +38,6 @@ map_beta_div <- function(Input_Image_File = FALSE,
                          Output_Dir = '',
                          window_size = 10,
                          TypePCA = 'SPCA',
-                         nb_partitions = 20,
                          nbclusters = 50,
                          Nb_Units_Ordin = 2000,
                          MinSun = 0.25, pcelim = 0.02, scaling = 'PCO', dimMDS = 3,
@@ -52,12 +50,13 @@ map_beta_div <- function(Input_Image_File = FALSE,
   Spectral_Species_Path <- SSDpathlist$Spectral_Species_Path
   Input_Image_File <- SSDpathlist$Input_Image_File
   nbclusters <- SSDpathlist$nbclusters
+  HDR_SSD <- read_ENVI_header(get_HDR_name(Spectral_Species_Path))
 
   # 2- compute beta diversity
   Beta <- compute_beta_metrics(ClusterMap_Path = Spectral_Species_Path,
                                MinSun = MinSun,
                                Nb_Units_Ordin = Nb_Units_Ordin,
-                               nb_partitions = nb_partitions,
+                               nb_partitions = HDR_SSD$bands,
                                nbclusters = nbclusters, pcelim = pcelim,
                                scaling = scaling, dimMDS = dimMDS,
                                nbCPU = nbCPU, MaxRAM = MaxRAM)
@@ -117,6 +116,51 @@ map_beta_div <- function(Input_Image_File = FALSE,
 #'   return(BetaNMDS_sel)
 #' }
 
+#' #' computes NMDS
+#' #
+#' #' @param MatBCdist BC dissimilarity matrix
+#' #' @param dimMDS numeric. number of dimensions of the NMDS
+#' #
+#' #' @return BetaNMDS_sel
+#' #' @importFrom future plan multiprocess multisession sequential
+#' #' @importFrom future.apply future_lapply
+#' #' @importFrom ecodist nmds
+#' #' @importFrom utils find
+#' #' @export
+#'
+#' compute_NMDS <- function(MatBCdist,dimMDS=3) {
+#'   nbiterNMDS <- 4
+#'   if (Sys.info()["sysname"] == "Windows") {
+#'     nbCoresNMDS <- 2
+#'   } else if (Sys.info()["sysname"] == "Linux") {
+#'     nbCoresNMDS <- 4
+#'   }
+#'   # multiprocess of spectral species distribution and alpha diversity metrics
+#'   # plan(multiprocess, workers = nbCoresNMDS) ## Parallelize using four cores
+#'   plan(multisession, workers = nbCoresNMDS) ## Parallelize using four cores
+#'   BetaNMDS <- future_lapply(MatBCdist,
+#'                             FUN = nmds,
+#'                             mindim = dimMDS, maxdim = dimMDS, nits = 1,
+#'                             future.packages = c("ecodist"))
+#'   plan(sequential)
+#'   # find iteration with minimum stress
+#'   Stress <- vector(length = nbiterNMDS)
+#'   for (i in 1:nbiterNMDS) {
+#'     Stress[i] <- BetaNMDS[[i]]$stress
+#'   }
+#'   print("Stress obtained for NMDS iterations:")
+#'   print(Stress)
+#'   print("Rule of thumb")
+#'   print("stress < 0.05 provides an excellent represention in reduced dimensions")
+#'   print("stress < 0.1 is great")
+#'   print("stress < 0.2 is good")
+#'   print("stress > 0.3 provides a poor representation")
+#'   MinStress <- find(Stress == min(Stress))
+#'   BetaNMDS_sel <- BetaNMDS[[MinStress]]$conf
+#'   BetaNMDS_sel <- data.frame(BetaNMDS_sel[[1]])
+#'   return(BetaNMDS_sel)
+#' }
+
 #' Identifies ordination coordinates based on nearest neighbors
 #'
 #' @param SSD_subset numeric. matrix corresponding to Spectral species distribution for a set of windows
@@ -372,7 +416,7 @@ compute_beta_metrics <- function(ClusterMap_Path,
 
   # create a Bray curtis dissimilarity matrix for each iteration
   print("compute BC dissimilarity for selected kernels")
-  Sample_Sel_list <- snow::splitCols(x = Sample_Sel,ncl = nbclusters)
+  Sample_Sel_list <- snow::splitCols(x = Sample_Sel,ncl = nb_partitions)
   plan(multisession, workers = nbCPU)
   handlers(global = TRUE)
   handlers("cli")
@@ -488,7 +532,7 @@ compute_BCdiss <- function(SSDList, pcelim=0.02) {
   #   SSD[[i]][which(SSD[[i]] < pcelim)] <- 0
   # }
 
-  SSD <- lapply(SSDList,FUN = Normalize_SSD, pcelim = pcelim)
+  SSD <- lapply(SSDList,FUN = normalize_SSD, pcelim = pcelim)
   # matrix of bray curtis dissimilarity (size = nb kernels x nb kernels)
   # Here use the package "dissUtils" to compute dissimilarity matrices sequentially
   MatBC <- diss(SSD[[1]], SSD[[2]], method = 'braycurtis')
@@ -616,7 +660,7 @@ getBCdiss <- function(Mat1, pcelim = 0.02, p = NULL){
 #'
 #' @return SSD numeric. matrix corresponding to normalized spectral species distribution
 #' @export
-Normalize_SSD <- function(SSDList, pcelim = 0.02){
+normalize_SSD <- function(SSDList, pcelim = 0.02){
 
   SumSpecies <- rowSums(SSDList)
   elim <- which(SumSpecies == 0)

R/Lib_MapFunctionalDiversity.R

@@ -96,18 +96,12 @@ map_functional_div <- function(Original_Image_File,Functional_File = FALSE,
 #' @importFrom stats sd
 #' @importFrom raster brick values nbands
 compute_Functional_metrics <- function(Functional_File, Functional_Map_Path, Selected_Features,
-                                       window_size, MinSun, nbCPU = FALSE, MaxRAM = FALSE) {
+                                       window_size, MinSun, nbCPU = 1, MaxRAM = 0.25) {
 
   ## read Functional_File and write Functional_Map_Path
   HDRname <- get_HDR_name(Functional_File)
   HDR <- read_ENVI_header(HDRname)
-  if (MaxRAM == FALSE) {
-    MaxRAM <- 0.25
-  }
   nbPieces_Min <- split_image(HDR, MaxRAM)
-  if (nbCPU == FALSE) {
-    nbCPU <- 1
-  }
   if (nbPieces_Min < nbCPU) {
     nbPieces_Min <- nbCPU
   }
@@ -168,15 +162,13 @@ compute_Functional_metrics <- function(Functional_File, Functional_Map_Path, Sel
 
   # multiprocess of spectral species distribution and alpha diversity metrics
   if (nbCPU>1){
-    plan(multiprocess, workers = nbCPU) ## Parallelize using four cores
-    Schedule_Per_Thread <- ceiling(nbPieces_Min / nbCPU)
+    plan(multisession, workers = nbCPU) ## Parallelize using four cores
     FUNCT_DIV <- future_lapply(ReadWrite,
                                FUN = Get_FunctionalMetrics_From_Traits, Functional_File = Functional_File, Selected_Features = Selected_Features,
                                MinMaxRaster = MinMaxRaster, HDR = HDR, HDR_Funct = HDR_Funct,
                                FunctIN_Format = FunctIN_Format, FunctOUT_Format = FunctOUT_Format,
                                ImgFormat = ImgFormat, window_size = window_size, MinSun = MinSun,
-                               Functional_Map_Path = Functional_Map_Path,  future.scheduling = Schedule_Per_Thread
-    )
+                               Functional_Map_Path = Functional_Map_Path)
     plan(sequential)
   } else {
     FUNCT_DIV <- lapply(ReadWrite, FUN = Get_FunctionalMetrics_From_Traits, Functional_File = Functional_File,