biodivMapR v1.9.7

## Addition - progress bar instead of messages - future: multisession instead of multiprocess
jbferet · Nov 7, 2022 · 656ac98 · 656ac98
1 parent 80b4539
commit 656ac98
Show file tree

Hide file tree

Showing 13 changed files with 269 additions and 105 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: biodivMapR
 Title: biodivMapR: an R package for a- and ß-diversity mapping using remotely-sensed images
-Version: 1.9.6.1
+Version: 1.9.7
 Authors@R: c(person(given = "Jean-Baptiste", 
 					family = "Feret", 
 					email = "jb.feret@teledetection.fr", 
@@ -42,7 +42,9 @@ Imports:
     vegan,
     zip,
     sf,
-    rgeos
+    rgeos,
+    progress,
+    progressr
 RoxygenNote: 7.2.1
 Suggests: 
     knitr,

diff --git a/NAMESPACE b/NAMESPACE
@@ -9,6 +9,7 @@ export(Write_Big_Image)
 export(Write_Image_NativeRes)
 export(ZipFile)
 export(apply_continuum_removal)
+export(apply_kmeans)
 export(build_image_from_list)
 export(center_reduce)
 export(change_resolution_HDR)
@@ -84,6 +85,7 @@ export(where_to_read)
 export(where_to_read_kernel)
 export(where_to_write_kernel)
 export(write_ENVI_header)
+export(write_PCA_raster)
 export(write_StarsStack)
 export(write_raster)
 import(stars)
@@ -96,6 +98,7 @@ importFrom(ecodist,nmds)
 importFrom(emstreeR,ComputeMST)
 importFrom(fields,rdist)
 importFrom(future,multiprocess)
+importFrom(future,multisession)
 importFrom(future,plan)
 importFrom(future,sequential)
 importFrom(future.apply,future_lapply)
@@ -109,6 +112,10 @@ importFrom(matrixStats,rowSds)
 importFrom(methods,as)
 importFrom(methods,is)
 importFrom(mmand,erode)
+importFrom(progress,progress_bar)
+importFrom(progressr,handlers)
+importFrom(progressr,progressor)
+importFrom(progressr,with_progress)
 importFrom(raster,brick)
 importFrom(raster,cellFromPolygon)
 importFrom(raster,cellFromXY)

diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,9 @@
+# biodivMapR v1.9.7
+
+## Addition
+- progress bar instead of messages
+- future: multisession instead of multiprocess
+
 # biodivMapR v1.9.6.1
 
 ## Addition

diff --git a/R/Lib_ImageProcess.R b/R/Lib_ImageProcess.R
@@ -660,7 +660,8 @@ get_BB_from_Vector <- function(path_raster,path_vector,Buffer = 0){
 #' Kind (Kernel index) and 'id' the sample ID to be used with the kernel
 #' @export
 
-get_random_subset_from_image <- function(ImPath, MaskPath, nb_partitions, Pix_Per_Partition, kernel=NULL,MaxRAM = 0.5) {
+get_random_subset_from_image <- function(ImPath, MaskPath, nb_partitions,
+                                         Pix_Per_Partition, kernel=NULL, MaxRAM = 0.5) {
 
   metarast <- raster::raster(ImPath)
   # updated raster package: do not use brick with 2D raster
@@ -1524,6 +1525,7 @@ write_ENVI_header <- function(HDR, HDRpath) {
 #' @param Image_Format list. description of data format corresponding to ENVI type
 #'
 #' @return None
+#' @importFrom progress progress_bar
 #' @export
 
 Write_Big_Image <- function(ImgWrite,ImagePath,HDR,Image_Format){
@@ -1535,8 +1537,11 @@ Write_Big_Image <- function(ImgWrite,ImagePath,HDR,Image_Format){
   )
   close(fidOUT)
   # for each piece of image
+  pb <- progress_bar$new(
+    format = 'Writing raster [:bar] :percent in :elapsedfull',
+    total = nbPieces, clear = FALSE, width= 100)
+
   for (i in 1:nbPieces) {
-    print(paste("Writing Image, piece #", i, "/", nbPieces))
     # read image and mask data
     Byte_Start <- SeqRead_Image$ReadByte_Start[i]
     Line_Start <- SeqRead_Image$Line_Start[i]
@@ -1564,6 +1569,7 @@ Write_Big_Image <- function(ImgWrite,ImagePath,HDR,Image_Format){
       writeBin(c(as.integer(ImgChunk)), fidOUT, size = Image_Format$Bytes, endian = .Platform$endian, useBytes = FALSE)
     }
     close(fidOUT)
+    pb$tick()
   }
   rm(ImgWrite)
   rm(ImgChunk)

diff --git a/R/Lib_MapAlphaDiversity.R b/R/Lib_MapAlphaDiversity.R
@@ -34,8 +34,8 @@
 map_alpha_div <- function(Input_Image_File=FALSE, Output_Dir='', window_size=10,
                           TypePCA = "SPCA", nbclusters = 50,
                           MinSun = 0.25, pcelim = 0.02,
-                          Index_Alpha = "Shannon", FullRes = TRUE,
-                          LowRes = FALSE, MapSTD = FALSE,
+                          Index_Alpha = "Shannon", FullRes = FALSE,
+                          LowRes = TRUE, MapSTD = TRUE,
                           nbCPU = FALSE, MaxRAM = FALSE,
                           ClassifMap = FALSE) {
 
@@ -198,8 +198,9 @@ compute_ALPHA_FromPlot <- function(SpectralSpecies_Plot,pcelim = 0.02){
 # @param MinSun numeric. minimum proportion of sunlit pixels required to consider plot
 #
 # @return list of mean and SD of alpha diversity metrics
-#' @importFrom future plan multiprocess sequential
+#' @importFrom future plan multiprocess multisession sequential
 #' @importFrom future.apply future_lapply
+#' @importFrom progressr progressor handlers with_progress
 #' @importFrom stats sd
 compute_alpha_metrics <- function(Spectral_Species_Path, window_size, nbclusters,
                                   MinSun, pcelim, nbCPU = FALSE, MaxRAM = FALSE, Index_Alpha = "Shannon") {
@@ -209,14 +210,18 @@ compute_alpha_metrics <- function(Spectral_Species_Path, window_size, nbclusters
   if (MaxRAM == FALSE) {
     MaxRAM <- 0.25
   }
-  nbPieces_Min <- split_image(HDR_SS, MaxRAM)
   if (nbCPU == FALSE) {
     nbCPU <- 1
   }
-  if (nbPieces_Min < nbCPU) {
-    nbPieces_Min <- nbCPU
-  }
-  SeqRead.SS <- where_to_read_kernel(HDR_SS, nbPieces_Min, window_size)
+  nbPieces <- split_image(HDR_SS, MaxRAM)
+  nbPieces <- nbCPU*(1+(nbPieces-1)%/%nbCPU)
+  # if (nbPieces < nbCPU) {
+  #   nbPieces <- nbCPU
+  # }
+  # if (nbPieces<10){
+  #   nbPieces <- 10
+  # }
+  SeqRead.SS <- where_to_read_kernel(HDR_SS, nbPieces, window_size)
 
   ##          prepare SS distribution map and corresponding sunlit map        ##
   # prepare SS distribution map
@@ -240,7 +245,7 @@ compute_alpha_metrics <- function(Spectral_Species_Path, window_size, nbclusters
   close(fidSSD)
   headerFpath <- paste(SSD_Path, ".hdr", sep = "")
   write_ENVI_header(HDR_SSD, headerFpath)
-  SeqWrite.SSD <- where_to_write_kernel(HDR_SS, HDR_SSD, nbPieces_Min, window_size)
+  SeqWrite.SSD <- where_to_write_kernel(HDR_SS, HDR_SSD, nbPieces, window_size)
 
   # prepare proportion of sunlit pixels from each spatial unit
   Sunlit_Path <- paste(SSD_Path, "_Sunlit", sep = "")
@@ -257,11 +262,11 @@ compute_alpha_metrics <- function(Spectral_Species_Path, window_size, nbclusters
   close(fidSunlit)
   headerFpath <- paste(Sunlit_Path, ".hdr", sep = "")
   write_ENVI_header(HDR_Sunlit, headerFpath)
-  SeqWrite.Sunlit <- where_to_write_kernel(HDR_SS, HDR_Sunlit, nbPieces_Min, window_size)
+  SeqWrite.Sunlit <- where_to_write_kernel(HDR_SS, HDR_Sunlit, nbPieces, window_size)
 
   # for each piece of image
   ReadWrite <- list()
-  for (i in 1:nbPieces_Min) {
+  for (i in 1:nbPieces) {
     ReadWrite[[i]] <- list()
     ReadWrite[[i]]$RW_SS <- ReadWrite[[i]]$RW_SSD <- ReadWrite[[i]]$RW_Sunlit <- list()
     ReadWrite[[i]]$RW_SS$Byte_Start <- SeqRead.SS$ReadByte_Start[i]
@@ -283,22 +288,35 @@ compute_alpha_metrics <- function(Spectral_Species_Path, window_size, nbclusters
 
   # 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)
-    ALPHA <- future_lapply(ReadWrite,
-                           FUN = convert_PCA_to_SSD, Spectral_Species_Path = Spectral_Species_Path,
-                           HDR_SS = HDR_SS, HDR_SSD = HDR_SSD, SS_Format = SS_Format, SSD_Format = SSD_Format,
-                           ImgFormat = ImgFormat, window_size = window_size, nbclusters = nbclusters, MinSun = MinSun,
-                           pcelim = pcelim, Index_Alpha = Index_Alpha, SSD_Path = SSD_Path, Sunlit_Path = Sunlit_Path,
-                           Sunlit_Format = Sunlit_Format, future.scheduling = Schedule_Per_Thread
-    )
+    plan(multisession, workers = nbCPU) ## Parallelize using four cores
+    # plan(multiprocess, workers = nbCPU) ## Parallelize using four cores
+    Schedule_Per_Thread <- ceiling(nbPieces / nbCPU)
+    handlers(global = TRUE)
+    handlers("progress")
+    with_progress({
+      p <- progressr::progressor(steps = nbPieces)
+      ALPHA <- future_lapply(ReadWrite,
+                             FUN = convert_PCA_to_SSD, Spectral_Species_Path = Spectral_Species_Path,
+                             HDR_SS = HDR_SS, HDR_SSD = HDR_SSD, SS_Format = SS_Format, SSD_Format = SSD_Format,
+                             ImgFormat = ImgFormat, window_size = window_size, nbclusters = nbclusters, MinSun = MinSun,
+                             pcelim = pcelim, Index_Alpha = Index_Alpha, SSD_Path = SSD_Path, Sunlit_Path = Sunlit_Path,
+                             Sunlit_Format = Sunlit_Format, p = p)
+      })
+
+
+    # ALPHA <- future_lapply(ReadWrite,
+    #                        FUN = convert_PCA_to_SSD, Spectral_Species_Path = Spectral_Species_Path,
+    #                        HDR_SS = HDR_SS, HDR_SSD = HDR_SSD, SS_Format = SS_Format, SSD_Format = SSD_Format,
+    #                        ImgFormat = ImgFormat, window_size = window_size, nbclusters = nbclusters, MinSun = MinSun,
+    #                        pcelim = pcelim, Index_Alpha = Index_Alpha, SSD_Path = SSD_Path, Sunlit_Path = Sunlit_Path,
+    #                        Sunlit_Format = Sunlit_Format, future.scheduling = Schedule_Per_Thread)
     plan(sequential)
   } else {
     ALPHA <- lapply(ReadWrite, FUN = convert_PCA_to_SSD, Spectral_Species_Path = Spectral_Species_Path,
                     HDR_SS = HDR_SS, HDR_SSD = HDR_SSD, SS_Format = SS_Format, SSD_Format = SSD_Format,
                     ImgFormat = ImgFormat, window_size = window_size, nbclusters = nbclusters, MinSun = MinSun,
                     pcelim = pcelim, Index_Alpha = Index_Alpha, SSD_Path = SSD_Path, Sunlit_Path = Sunlit_Path,
-                    Sunlit_Format = Sunlit_Format)
+                    Sunlit_Format = Sunlit_Format,p= NULL)
   }
 
   # create ful map from chunks
@@ -348,13 +366,15 @@ compute_alpha_metrics <- function(Spectral_Species_Path, window_size, nbclusters
 # @param SSD_Path
 # @param Sunlit_Path
 # @param Sunlit_Format
+# @param p
 #
 # @param
 # @param
 # @return
 convert_PCA_to_SSD <- function(ReadWrite, Spectral_Species_Path, HDR_SS, HDR_SSD,
                                SS_Format, SSD_Format, ImgFormat, window_size, nbclusters,
-                               MinSun, pcelim, Index_Alpha, SSD_Path, Sunlit_Path, Sunlit_Format) {
+                               MinSun, pcelim, Index_Alpha, SSD_Path, Sunlit_Path, Sunlit_Format,
+                               p = NULL) {
   SS_Chunk <- read_BIL_image_subset(
     Spectral_Species_Path, HDR_SS,
     ReadWrite$RW_SS$Byte_Start, ReadWrite$RW_SS$lenBin,
@@ -398,6 +418,9 @@ convert_PCA_to_SSD <- function(ReadWrite, Spectral_Species_Path, HDR_SS, HDR_SSD
   Simpson_SD_Chunk <- apply(SSD_Alpha$Simpson, 1:2, sd)
   rm(SSD_Alpha)
   gc()
+  if (!is.null(p)){
+    p()
+  }
   my_list <- list(
     "Shannon" = Shannon_Mean_Chunk, "Fisher" = Fisher_Mean_Chunk, "Simpson" = Simpson_Mean_Chunk,
     "Shannon.SD" = Shannon_SD_Chunk, "Fisher.SD" = Fisher_SD_Chunk, "Simpson.SD" = Simpson_SD_Chunk

diff --git a/R/Lib_MapBetaDiversity.R b/R/Lib_MapBetaDiversity.R
@@ -38,8 +38,8 @@
 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', FullRes = TRUE,
-                         LowRes = FALSE, nbCPU = 1, MaxRAM = 0.25,
+                         pcelim = 0.02, scaling = 'PCO', FullRes = FALSE,
+                         LowRes = TRUE, nbCPU = 1, MaxRAM = 0.25,
                          ClassifMap = FALSE, dimMDS=3) {
 
   if (ClassifMap == FALSE){
@@ -110,7 +110,7 @@ map_beta_div <- function(Input_Image_File=FALSE, Output_Dir='', window_size=10,
 #' @param dimMDS numeric. number of dimensions of the NMDS
 #
 #' @return BetaNMDS_sel
-#' @importFrom future plan multiprocess sequential
+#' @importFrom future plan multiprocess multisession sequential
 #' @importFrom future.apply future_lapply
 #' @importFrom ecodist nmds
 #' @importFrom utils find
@@ -124,7 +124,8 @@ compute_NMDS <- function(MatBCdist,dimMDS=3) {
     nbCoresNMDS <- 4
   }
   # multiprocess of spectral species distribution and alpha diversity metrics
-  plan(multiprocess, workers = nbCoresNMDS) ## Parallelize using four cores
+  # 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
@@ -158,37 +159,46 @@ compute_NMDS <- function(MatBCdist,dimMDS=3) {
 #'
 #' @return Ordination_est coordinates of each spatial unit in ordination space
 #' @importFrom snow splitRows
-#' @importFrom future plan multiprocess sequential
+#' @importFrom future plan multiprocess multisession sequential
 #' @importFrom future.apply future_lapply
+#' @importFrom progressr progressor handlers with_progress
 #' @export
 
 ordination_to_NN <- function(Beta_Ordination_sel, SSD_Path, Sample_Sel, coordTotSort,
                              nb_partitions, nbclusters, pcelim, nbCPU = 1) {
   nb_Sunlit <- dim(coordTotSort)[1]
-  # define number of samples to be sampled each time during paralle processing
+  # define number of samples to be sampled each time during parallel processing
   nb_samples_per_sub <- round(1e7 / dim(Sample_Sel)[1])
   # number of paralle processes to run
   nb.sub <- round(nb_Sunlit / nb_samples_per_sub)
   if (nb.sub == 0) nb.sub <- 1
-  id.sub <- splitRows(as.matrix(seq(1, nb_Sunlit, by = 1), ncol = 1), ncl = nb.sub)
+  id.sub <- snow::splitRows(as.matrix(seq(1, nb_Sunlit, by = 1), ncol = 1), ncl = nb.sub)
   # compute ordination coordinates from each subpart
   Nb_Units_Ordin <- dim(Sample_Sel)[1]
   if (nbCPU>1){
-    plan(multiprocess, workers = nbCPU) ## Parallelize using four cores
-    Schedule_Per_Thread <- ceiling(nb.sub / nbCPU)
-    OutPut <- future_lapply(id.sub,
-                            FUN = ordination_parallel, coordTotSort = coordTotSort, SSD_Path = SSD_Path,
-                            Sample_Sel = Sample_Sel, Beta_Ordination_sel = Beta_Ordination_sel, Nb_Units_Ordin = Nb_Units_Ordin,
-                            nb_partitions = nb_partitions, nbclusters = nbclusters, pcelim = pcelim,
-                            future.scheduling = Schedule_Per_Thread,
-                            future.packages = c("vegan", "dissUtils", "R.utils", "tools", "snow", "matlab")
-    )
+    # plan(multiprocess, workers = nbCPU) ## Parallelize using four cores
+    plan(multisession, workers = nbCPU) ## Parallelize using four cores
+    handlers(global = TRUE)
+    handlers("progress")
+    with_progress({
+      p <- progressr::progressor(steps = nb.sub)
+      OutPut <- future_lapply(id.sub,
+                              FUN = ordination_parallel, coordTotSort = coordTotSort, SSD_Path = SSD_Path,
+                              Sample_Sel = Sample_Sel, Beta_Ordination_sel = Beta_Ordination_sel, Nb_Units_Ordin = Nb_Units_Ordin,
+                              nb_partitions = nb_partitions, nbclusters = nbclusters, pcelim = pcelim, p,
+                              future.packages = c("vegan", "dissUtils", "R.utils", "tools", "snow", "matlab"))
+      })
     plan(sequential)
   } else {
-    OutPut <- lapply(id.sub, FUN = ordination_parallel, coordTotSort = coordTotSort,
-                     SSD_Path = SSD_Path, Sample_Sel = Sample_Sel,
-                     Beta_Ordination_sel = Beta_Ordination_sel, Nb_Units_Ordin = Nb_Units_Ordin,
-                     nb_partitions = nb_partitions, nbclusters = nbclusters, pcelim = pcelim)
+    handlers(global = TRUE)
+    handlers("progress")
+    with_progress({
+      p <- progressr::progressor(steps = nb.sub)
+      OutPut <- lapply(id.sub, FUN = ordination_parallel, coordTotSort = coordTotSort,
+                       SSD_Path = SSD_Path, Sample_Sel = Sample_Sel,
+                       Beta_Ordination_sel = Beta_Ordination_sel, Nb_Units_Ordin = Nb_Units_Ordin,
+                       nb_partitions = nb_partitions, nbclusters = nbclusters, pcelim = pcelim, p)
+    })
   }
   Ordination_est <- do.call("rbind", OutPut)
   gc()
@@ -206,12 +216,13 @@ ordination_to_NN <- function(Beta_Ordination_sel, SSD_Path, Sample_Sel, coordTot
 #' @param nb_partitions numeric. Number of partitions (repetitions) to be computed then averaged.
 #' @param nbclusters numeric. Number of clusters defined in k-Means.
 #' @param pcelim numeric. Minimum contribution in percents required for a spectral species
+#' @param p function.
 #
 #' @return OutPut list of mean and SD of alpha diversity metrics
 #' @export
 
 ordination_parallel <- function(id.sub, coordTotSort, SSD_Path, Sample_Sel, Beta_Ordination_sel,
-                                Nb_Units_Ordin, nb_partitions, nbclusters, pcelim) {
+                                Nb_Units_Ordin, nb_partitions, nbclusters, pcelim, p = NULL) {
 
   # get Spectral species distribution
   coordPix <- coordTotSort[id.sub, ]
@@ -231,6 +242,8 @@ ordination_parallel <- function(id.sub, coordTotSort, SSD_Path, Sample_Sel, Beta
   # get the knn closest neighbors from each kernel
   knn <- 3
   OutPut <- compute_NN_from_ordination(MatBCtmp, knn, Beta_Ordination_sel)
+
+  p()
   return(OutPut)
 }