--- a/.shed.yml	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,14 +0,0 @@
-categories:
-    - Ecology
-owner: ecology
-remote_repository_url: https://github.com/Marie59/Sentinel_2A/srs_tools
-long_description: |
-  Compute biodiversity indicators for remote sensing data from Sentinel 2
-type: unrestricted
-auto_tool_repositories:
-  name_template: "{{ tool_id }}"
-  description_template: "Wrapper for srs biodiversity indicators tool: {{ tool_name }}."
-suite:
-  name: "sentinel_2"
-  description: "A suite of biodiversity indicators for Sentinel 2 data"
-  type: unrestricted

--- a/Lib_preprocess_S2.r	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,1335 +0,0 @@
-# == == == == == == == == == == == == == == == == == == == == == == == == == == ==
-# preprocS2
-# Lib_preprocess_S2.R
-# == == == == == == == == == == == == == == == == == == == == == == == == == == ==
-# PROGRAMMERS:
-# Jean-Baptiste FERET <jb.feret@teledetection.fr>
-# Copyright 2021/08 Jean-Baptiste FERET
-# == == == == == == == == == == == == == == == == == == == == == == == == == == ==
-# This Library contains functions to preprocess Sentinel-2 images downloaded from
-# different data hubs, such as THEIA, PEPS or SCIHUB
-# == == == == == == == == == == == == == == == == == == == == == == == == == == ==
-
-#" This function adjusts information from ENVI header
-#"
-#" @param dsn character. path where to store the stack
-#" @param bands list. should include "bandname", and if possible "wavelength"
-#" @param sensor character. Name of the sensor used to acquire the image
-#" @param stretch boolean. Set TRUE to get 10% stretching at display for reflectance, mentioned in hdr only
-#"
-#" @return None
-#" @importFrom utils read.table
-#" @importFrom raster hdr raster
-#" @export
-adjust_envi_hdr <- function(dsn, bands, sensor = "Unknown", stretch = FALSE) {
-
-  # Edit hdr file to add metadata
-  hdr <- read_envi_header(get_hdr_name(dsn))
-  hdr$`band names` <- bands$bandname
-  if (length(bands$wavelength) == length(bands$bandname)) {
-    hdr$wavelength <- bands$wavelength
-  }else {
-    hdr$wavelength <- NULL
-  }
-  if (stretch == TRUE) {
-    hdr$`default stretch` <- "0.000000 1000.000000 linear"
-  }
-  hdr$`z plot range` <- NULL
-  hdr$`data ignore value` <- "-Inf"
-  hdr$`sensor type` <- sensor
-  write_envi_header(hdr = hdr, hdrpath = get_hdr_name(dsn))
-
-  # remove unnecessary files
-  file2remove <- paste(dsn, ".aux.xml", sep = "")
-  if (file.exists(file2remove)) file.remove(file2remove)
-  file2remove <- paste(dsn, ".prj", sep = "")
-  if (file.exists(file2remove)) file.remove(file2remove)
-  file2remove <- paste(dsn, ".stx", sep = "")
-  if (file.exists(file2remove)) file.remove(file2remove)
-  return(invisible())
-}
-
-#" This function saves reflectance files
-#"
-#" @param s2sat character. Sentinel-2 mission ("2A" or "2B")
-#" @param tile_s2 character. S2 tile name (2 numbers + 3 letters)
-#" @param dateacq_s2 double. date of acquisition
-#"
-#" @return s2mission character. name of the S2 mission (2A or 2B)
-#" @importFrom sen2r safe_getMetadata check_scihub_connection s2_list
-#" @export
-check_s2mission <- function(s2sat, tile_s2, dateacq_s2) {
-
-  # is mission already defined by user?
-  if (!is.null(s2sat)) {
-    if (s2sat == "2A") {
-      s2mission <- "2A"
-    }else if (s2sat == "2B") {
-      s2mission <- "2B"
-    }else {
-      message("Could not identify if image from Sentinel-2A or -2B")
-      message("Defining central wavelength of spectral bands based on S2A")
-      s2mission <- "2A"
-    }
-  }else {
-    message("Could not identify if image from Sentinel-2A or -2B")
-    message("Defining central wavelength of spectral bands based on S2A")
-    s2mission <- "2A"
-  }
-  return(s2mission)
-}
-
-#" this function aims at computing directory size
-#" @param path character. path for directory
-#" @param recursive boolean . set T if recursive
-#"
-#" @return size_files numeric. size in bytes
-#" - image stack
-#" - path for individual band files corresponding to the stack
-#" - path for vector (reprojected if needed)
-#"
-#" @importFrom raster raster
-#" @importFrom tools file_path_sans_ext file_ext
-#" @export
-dir_size <- function(path, recursive = TRUE) {
-  stopifnot(is.character(path))
-  files <- list.files(path, full.names = TRUE, recursive = recursive)
-  vect_size <- sapply(files, function(x) file.size(x))
-  size_files <- sum(vect_size)
-  return(size_files)
-}
-
-#" This function reads S2 data from L2A directories downloaded from
-#" various data hubs including THEIA, PEPS & SCIHUB (SAFE format & LaSRC)
-#" @param path_dir_s2 character. path for S2 directory
-#" @param path_vector character. path for vector file
-#" @param s2source character. type of directory format (depends on atmospheric correction: SAFE produced from Sen2Cor)
-#" @param resolution numeric. buffer applied to vector file (in meters)
-#" @param interpolation character. method for resampling. default = "bilinear"
-#" @param fre_sre character. SRE or FRE products from THEIA
-#"
-#" @return listout list.
-#" - image stack
-#" - path for individual band files corresponding to the stack
-#" - path for vector (reprojected if needed)
-#"
-#" @importFrom raster raster
-#" @importFrom tools file_path_sans_ext file_ext
-#" @export
-extract_from_s2_l2a <- function(path_dir_s2, path_vector = NULL, s2source = "SAFE",
-                                resolution = 10, interpolation = "bilinear", fre_sre = "FRE") {
-  # Get list of paths corresponding to S2 bands and depending on S2 directory
-  s2_bands <- get_s2_bands(path_dir_s2 = path_dir_s2,
-                           s2source = s2source,
-                           resolution = resolution,
-                           fre_sre = fre_sre)
-
-  if (length(s2_bands$s2bands_10m) > 0) {
-    rastmp <- raster::raster(s2_bands$s2bands_10m[[1]])
-  } else if (length(s2_bands$s2bands_20m) > 0) {
-    rastmp <- raster::raster(s2_bands$s2bands_20m[[1]])
-  }
-  # check if vector and raster share the same projection. if not, re-project vector
-  if (!is.null(path_vector)) {
-    raster_proj <- raster::projection(rastmp)
-    path_vector_reproj <- paste(tools::file_path_sans_ext(path_vector), "_reprojected.shp", sep = "")
-    path_vector <- reproject_shp(path_vector_init = path_vector,
-                                 newprojection = raster_proj,
-                                 path_vector_reproj = path_vector_reproj)
-  }
-  # Extract data corresponding to the vector footprint (if provided) & resample data if needed
-  if (length(s2_bands$s2bands_10m) > 0) {
-    stack_10m <- read_s2bands(s2_bands = s2_bands$s2bands_10m, path_vector = path_vector,
-                              resampling = 1, interpolation = interpolation)
-  }
-  if (length(s2_bands$s2bands_20m) > 0) {
-    if (resolution == 10 && s2source != "LaSRC") {
-      resampling <- 2
-    }else {
-      resampling <- 1
-    }
-    stack_20m <- read_s2bands(s2_bands = s2_bands$s2bands_20m, path_vector = path_vector,
-                              resampling = resampling, interpolation = interpolation)
-  }
-  # get full stack including 10m and 20m spatial resolution
-  if (length(s2_bands$s2bands_10m) > 0 && length(s2_bands$s2bands_20m) > 0) {
-    diffxstart <- attributes(stack_10m)$dimensions[[1]]$from - attributes(stack_20m)$dimensions[[1]]$from
-    diffxstop <- attributes(stack_10m)$dimensions[[1]]$to - attributes(stack_20m)$dimensions[[1]]$to
-    diffystart <- attributes(stack_10m)$dimensions[[2]]$from - attributes(stack_20m)$dimensions[[2]]$from
-    diffystop <- attributes(stack_10m)$dimensions[[2]]$to - attributes(stack_20m)$dimensions[[2]]$to
-    if (!diffxstop == 0) {
-      # size of 20m > size of 10m --> reduce 20m
-      # size of 10m > size of 20m --> reduce 10m
-      if (diffxstop > 0) {
-        stack_10m <- stack_10m[, 1:(dim(stack_10m)[1] - diffxstop), , ]
-      }else if (diffxstop < 0) {
-        stack_20m <- stack_20m[, 1:(dim(stack_20m)[1] + diffxstop), , ]
-      }
-    }
-    if (!diffystop == 0) {
-      if (diffystop > 0) {
-        stack_10m <- stack_10m[, , 1:(dim(stack_10m)[2] - diffystop), ]
-      }else if (diffystop < 0) {
-        stack_20m <- stack_20m[, , 1:(dim(stack_20m)[2] + diffystop), ]
-      }
-    }
-    if (!diffxstart == 0) {
-      if (diffxstart > 0) {
-        stack_20m <- stack_20m[, (1 + diffxstart):dim(stack_20m)[1], , ]
-      }else if (diffxstart < 0) {
-        stack_10m <- stack_10m[, (1 - diffxstart):dim(stack_10m)[1], , ]
-      }
-    }
-    if (!diffystart == 0) {
-      if (diffystart > 0) {
-        stack_20m <- stack_20m[, , (1 + diffystart):dim(stack_20m)[2], ]
-      }else if (diffystart < 0) {
-        stack_10m <- stack_10m[, , (1 - diffystart):dim(stack_10m)[2], ]
-      }
-    }
-    # reorder bands with increasing wavelength
-    s2bands <- c("B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B11", "B12", "Cloud")
-    namebands <- c(names(s2_bands$s2bands_10m), names(s2_bands$s2bands_20m))
-    reorder_bands <- match(s2bands, namebands)
-    namebands <- namebands[reorder_bands]
-    listfiles <- c(stack_10m$attr, stack_20m$attr)[reorder_bands]
-
-    # adjust size to initial vector footprint without buffer
-    # --> buffer is needed in order to ensure that extraction following
-    # footprint of vector matches for images of different spatial resolution
-    # get bounding box corresponding to footprint of image or image subset
-    bb_xycoords <- get_bb(path_raster = listfiles[1],
-                          path_vector = path_vector, buffer = 0)
-
-    # prepare reading data for extent defined by bounding box
-    nxoff <- bb_xycoords$UL$col
-    nyoff <- bb_xycoords$UL$row
-    nxsize <- bb_xycoords$UR$col - bb_xycoords$UL$col + 1
-    nysize <- bb_xycoords$LR$row - bb_xycoords$UR$row + 1
-    nbufxsize <- nxsize
-    nbufysize <- nysize
-    s2_stack <- stars::read_stars(listfiles, along = "band",
-                                  RasterIO = list(nXOff = nxoff, nYOff = nyoff,
-                                                  nXSize = nxsize, nYSize = nysize,
-                                                  nBufXSize = nbufxsize, nBufYSize = nbufysize,
-                                                  resample = "nearest_neighbour"), proxy = TRUE)
-
-
-    names(s2_stack$attr) <- namebands
-  }else if (length(s2_bands$s2bands_10m) > 0) {
-    s2_stack <- stack_10m
-    namebands <- names(s2_bands$s2bands_10m)
-    names(s2_stack$attr) <- namebands
-  }else if (length(s2_bands$s2bands_20m) > 0) {
-    s2_stack <- stack_20m
-    namebands <- names(s2_bands$s2bands_20m)
-    names(s2_stack$attr) <- namebands
-  }
-
-  listout <- list("s2_stack" = s2_stack, "s2_bands" = s2_bands, "path_vector" = path_vector,
-                  "namebands" = namebands)
-  return(listout)
-}
-
-#" This function gets coordinates of a bounding box defined by a vector (optional) and a raster
-#"
-#" @param path_raster character. path for raster file
-#" @param path_vector character. path for vector file
-#" @param buffer numeric. buffer applied to vector file (in meters)
-#"
-#" @return bb_xycoords list. Coordinates (in pixels) of the upper/lower right/left corners of bounding box
-#" @export
-get_bb <- function(path_raster, path_vector = NULL, buffer = 0) {
-
-  if (!is.null(path_vector)) {
-    # get bounding box with a 50m buffer in order to allow for interpolation
-    bb_xycoords <- get_bb_from_vector(path_raster = path_raster,
-                                      path_vector = path_vector,
-                                      buffer = buffer)
-  }else if (is.null(path_vector)) {
-    bb_xycoords <- get_bb_from_fullimage(path_raster)
-  }
-  return(bb_xycoords)
-}
-
-#" This function gets extreme coordinates of a bounding box corresponding to a full image
-#"
-#" @param path_raster character. path for raster file
-#"
-#" @return bb_xycoords list. Coordinates (in pixels) of the upper/lower right/left corners of bounding box
-#" @importFrom raster raster
-#" @export
-get_bb_from_fullimage <- function(path_raster) {
-  # get raster coordinates corresponding to Full image
-  rasterobj <- raster::raster(path_raster)
-  bb_xycoords <- list()
-  bb_xycoords[["UL"]] <- data.frame("row" = 1, "col" = 1)
-  bb_xycoords[["UR"]] <- data.frame("row" = 1, "col" = dim(rasterobj)[2])
-  bb_xycoords[["LL"]] <- data.frame("row" = dim(rasterobj)[1], "col" = 1)
-  bb_xycoords[["LR"]] <- data.frame("row" = dim(rasterobj)[1], "col" = dim(rasterobj)[2])
-  return(bb_xycoords)
-}
-
-#" This gets bounding box corresponding to a vector from a raster (UL, UR, LL, LR corners)
-#"
-#" @param path_raster character. path for raster file
-#" @param path_vector character. path for vector file
-#" @param buffer numeric. buffer applied to vector file (in meters)
-#"
-#" @return bb_xycoords list. Coordinates (in pixels) of the upper/lower right/left corners of bounding box
-#" @importFrom sf st_read st_bbox st_crop
-#" @importFrom rgeos gbuffer bbox2SP
-#" @importFrom sp SpatialPoints bbox
-#" @importFrom raster projection extract extent raster
-#" @importFrom methods as
-#" @export
-get_bb_from_vector <- function(path_raster, path_vector, buffer = 0) {
-
-  data_raster <- raster::raster(path_raster)
-  # extract BB coordinates from vector
-  bb_vector <- rgeos::gbuffer(spgeom = as(sf::st_read(dsn = path_vector, quiet = TRUE), "Spatial"),
-                              width = buffer, byid = TRUE)
-  # extract BB coordinates from raster
-  bb_raster <- rgeos::bbox2SP(bbox = bbox(data_raster))
-  # compute intersection
-  intersect <- rgeos::gIntersection(bb_vector, bb_raster)
-  bbext <- raster::extent(intersect)
-  xmin <- bbext[1]
-  xmax <- bbext[2]
-  ymin <- bbext[3]
-  ymax <- bbext[4]
-  # get coordinates of bounding box corresponding to vector
-  corners <- list()
-  corners[["UR"]] <- sp::SpatialPoints(coords = cbind(xmax, ymax))
-  corners[["LR"]] <- sp::SpatialPoints(coords = cbind(xmax, ymin))
-  corners[["UL"]] <- sp::SpatialPoints(coords = cbind(xmin, ymax))
-  corners[["LL"]] <- sp::SpatialPoints(coords = cbind(xmin, ymin))
-  raster::projection(corners[["UL"]]) <- raster::projection(corners[["UR"]]) <-
-    raster::projection(corners[["LL"]]) <- raster::projection(corners[["LR"]]) <-
-    raster::projection(sf::st_read(dsn = path_vector, quiet = TRUE))
-  # get coordinates for corners of bounding box
-  bb_xycoords <- list()
-  for (corner in names(corners)) {
-    ex_df <- as.data.frame(raster::extract(data_raster, corners[[corner]], cellnumbers = TRUE))
-    colrow <- ind2sub(data_raster, ex_df$cell)
-    bb_xycoords[[corner]] <- data.frame("row" = colrow$row, "col" = colrow$col)
-  }
-  return(bb_xycoords)
-}
-
-#" get hdr name from image file name, assuming it is BIL format
-#"
-#" @param impath path of the image
-#"
-#" @return corresponding hdr
-#" @import tools
-#" @export
-get_hdr_name <- function(impath) {
-  if (tools::file_ext(impath) == "") {
-    impathhdr <- paste(impath, ".hdr", sep = "")
-  }else if (tools::file_ext(impath) == "bil") {
-    impathhdr <- gsub(".bil", ".hdr", impath)
-  }else if (tools::file_ext(impath) == "zip") {
-    impathhdr <- gsub(".zip", ".hdr", impath)
-  }else {
-    impathhdr <- paste(tools::file_path_sans_ext(impath), ".hdr", sep = "")
-  }
-
-  if (!file.exists(impathhdr)) {
-    message("WARNING : COULD NOT FIND hdr FILE")
-    print(impathhdr)
-    message("Process may stop")
-  }
-  return(impathhdr)
-}
-
-#" This function returns path for the spectral bands to be used
-#"
-#" @param path_dir_s2 character. Path for the directory containing S2 data. either L2A .SAFE S2 file or THEIA directory
-#" @param s2source character. defines if data comes from SciHub as SAFE directory, from THEIA or from LaSRC
-#" @param resolution numeric. spatial resolution of the final image: 10m or 20m
-#" @param fre_sre character. SRE or FRE products from THEIA
-#"
-#" @return listbands list. contains path for spectral bands corresponding to 10m and 20m resolution
-#" @export
-get_s2_bands <- function(path_dir_s2, s2source = "SAFE", resolution = 10, fre_sre = "FRE") {
-
-  if (s2source == "SAFE" || s2source == "Sen2Cor") {
-    listbands <- get_s2_bands_from_sen2cor(path_dir_s2 = path_dir_s2, resolution = resolution)
-  }else if (s2source == "THEIA") {
-    listbands <- get_s2_bands_from_theia(path_dir_s2 = path_dir_s2, resolution = resolution,
-                                         fre_sre = fre_sre)
-  }else if (s2source == "LaSRC") {
-    listbands <- get_s2_bands_from_lasrc(path_dir_s2 = path_dir_s2, resolution = resolution)
-  }else {
-    message("The data source (Atmospheric correction) for Sentinel-2 image is unknown")
-    message("Please provide S2 images from one of the following data sources:")
-    message("- LaSRC (atmospheric correction: LaSRC)")
-    message("- THEIA (atmospheric correction: MAJA)")
-    message("- SAFE (atmospheric correction: Sen2Cor)")
-    s2bands_10m <- s2bands_20m <- granule <- mtdfile <- metadata_msi <- metadata_lasrc <- NULL
-    listbands <- list("s2bands_10m" = s2bands_10m, "s2bands_20m" = s2bands_20m, "GRANULE" = granule,
-                      "metadata" = mtdfile, "metadata_MSI" = metadata_msi,
-                      "metadata_lasrc" = metadata_lasrc)
-  }
-  return(listbands)
-}
-
-#" This function returns path for the spectral bands in SAFE / sen2Cor directory
-#"
-#" @param path_dir_s2 character. Path for the SAFE directory containing S2 data
-#" @param resolution numeric. spatial resolution of the final image: 10m or 20m
-#"
-#" @return listbands list. contains path for spectral bands corresponding to 10m and 20m resolution, as well name of as granule
-#" @export
-get_s2_bands_from_sen2cor <- function(path_dir_s2, resolution = 10) {
-  # build path for all bands
-  if (resolution == 10) {
-    b10m <- c("B02", "B03", "B04", "B08")
-    b20m <- c("B05", "B06", "B07", "B8A", "B11", "B12")
-  }else {
-    b10m <- c()
-    b20m <- c("B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B11", "B12")
-  }
-  # get granule directory & path for corresponding metadata XML file
-  granule <- list.dirs(list.dirs(path_dir_s2, recursive = FALSE)[grep(pattern = "GRANULE",
-                                                                     x = list.dirs(path_dir_s2, recursive = FALSE))], recursive = FALSE)
-  mtdfile <- file.path(granule, "MTD_TL.xml")
-  if (file.exists(file.path(path_dir_s2, "MTD_MSIL2A.xml"))) {
-    mtd_msi_file <- file.path(path_dir_s2, "MTD_MSIL2A.xml")
-  } else {
-    mtd_msi_file <- NULL
-  }
-
-  # Define path for bands
-  s2bands_20m_dir <- file.path(granule, "IMG_DATA", "R20m")
-  s2bands_10m_dir <- file.path(granule, "IMG_DATA", "R10m")
-  s2bands_10m <- s2bands_20m <- list()
-  for (band in b20m) {
-    s2bands_20m[[band]] <- file.path(s2bands_20m_dir, list.files(s2bands_20m_dir, pattern = band))
-  }
-  for (band in b10m) {
-    s2bands_10m[[band]] <- file.path(s2bands_10m_dir, list.files(s2bands_10m_dir, pattern = band))
-  }
-  # get cloud mask
-  cloud <- "MSK_CLDPRB_20m"
-  cloud_20m_dir <- file.path(granule, "QI_DATA")
-  s2bands_20m[["Cloud"]] <- file.path(cloud_20m_dir, list.files(cloud_20m_dir, pattern = cloud))
-  listbands <- list("s2bands_10m" = s2bands_10m,
-                    "s2bands_20m" = s2bands_20m,
-                    "GRANULE" = granule,
-                    "metadata" = mtdfile,
-                    "metadata_MSI" = mtd_msi_file,
-                    "metadata_lasrc" = NULL)
-  return(listbands)
-}
-
-#" This function returns path for the spectral bands in LaSRC directory
-#"
-#" @param path_dir_s2 character. Path for the SAFE directory containing S2 data
-#" @param resolution numeric. spatial resolution of the final image: 10m or 20m
-#"
-#" @return listbands list. contains path for spectral bands corresponding to 10m and 20m resolution, as well name of as granule
-#" @importFrom stringr str_subset
-#" @export
-get_s2_bands_from_lasrc <- function(path_dir_s2, resolution = 10) {
-
-  # get granule directory & path for corresponding metadata XML file
-  granule <- path_dir_s2
-  mtdfile <- file.path(granule, "MTD_TL.xml")
-  if (file.exists(file.path(path_dir_s2, "MTD_MSIL1C.xml"))) {
-    mtd_msi_file <- file.path(path_dir_s2, "MTD_MSIL1C.xml")
-  } else {
-    mtd_msi_file <- NULL
-  }
-
-  # build path for all bands
-  b10m <- c("band2", "band3", "band4", "band5", "band6", "band7", "band8", "band8a", "band11", "band12")
-  b10m_standard <- c("B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B11", "B12")
-  # Define path for bands
-  s2bands_10m <- s2bands_20m <- list()
-  for (i in 1:seq_along(b10m)) {
-    s2bands_10m[[b10m_standard[i]]] <- file.path(path_dir_s2,
-                                                 list.files(path_dir_s2,
-                                                            pattern = paste(b10m[i], ".tif", sep = "")))
-  }
-
-  # get metadata file containing offset
-  mtd_lasrc <- str_subset(list.files(path_dir_s2, pattern = "S2"), ".xml$")
-  if (file.exists(file.path(path_dir_s2, mtd_lasrc))) {
-    metadata_lasrc <- file.path(path_dir_s2, mtd_lasrc)
-  } else {
-    metadata_lasrc <- NULL
-  }
-  # get cloud mask
-  cloud <- "CLM"
-  s2bands_10m[["Cloud"]] <- file.path(path_dir_s2, list.files(path_dir_s2, pattern = cloud))
-  listbands <- list("s2bands_10m" = s2bands_10m,
-                    "s2bands_20m" = s2bands_20m,
-                    "GRANULE" = granule,
-                    "metadata" = mtdfile,
-                    "metadata_MSI" = mtd_msi_file,
-                    "metadata_lasrc" = metadata_lasrc)
-  return(listbands)
-}
-
-#" This function returns path for the spectral bands in THEIA directory
-#"
-#" @param path_dir_s2 character. Path for the SAFE directory containing S2 data
-#" @param resolution numeric. spatial resolution of the final image: 10m or 20m
-#" @param fre_sre character. SRE or FRE products from THEIA
-#"
-#" @return listbands list. contains path for spectral bands corresponding to 10m and 20m resolution, as well name of as granule
-#" @export
-get_s2_bands_from_theia <- function(path_dir_s2, resolution = 10, fre_sre = "FRE") {
-
-  # build path for all bands
-  if (resolution == 10) {
-    b10m <- c("B02", "B03", "B04", "B08")
-    b20m <- c("B05", "B06", "B07", "B8A", "B11", "B12")
-  } else {
-    b10m <- c()
-    b20m <- c("B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B11", "B12")
-  }
-
-  # get path_tile_s2 directory & path for corresponding metadata XML file
-  path_tile_s2 <- list.dirs(path_dir_s2, recursive = FALSE)
-  files_tile_s2 <- list.files(path_tile_s2, recursive = FALSE)
-  mtdfile <- file.path(path_tile_s2, files_tile_s2[grep(pattern = "MTD_ALL.xml", x = files_tile_s2)])
-
-  # Define path for bands
-  s2bands_10m_dir <- s2bands_20m_dir <- path_tile_s2
-  s2bands_10m <- s2bands_20m <- list()
-  for (band in b20m) {
-    band_20m_pattern <- paste0(gsub("0", "", band), ".tif") # for THEAI band 2 is "B2" ("B02" for SAFE)
-    list_files_20m <- list.files(s2bands_20m_dir, pattern = band_20m_pattern)
-    s2bands_20m[[band]] <- file.path(s2bands_20m_dir, list_files_20m)[grep(pattern = fre_sre,
-                                                                           x = file.path(s2bands_20m_dir, list_files_20m))]
-  }
-  for (band in b10m) {
-    band_10m_pattern <- paste0(gsub("0", "", band), ".tif") # for THEAI band 2 is "B2" ("B02" for SAFE)
-    list_files_10m <- list.files(s2bands_10m_dir, pattern = band_10m_pattern)
-    s2bands_10m[[band]] <- file.path(s2bands_10m_dir, list_files_10m)[grep(pattern = fre_sre,
-                                                                           x = file.path(s2bands_10m_dir, list_files_10m))]
-  }
-
-  # get cloud mask 10m
-  cloud_10m <- "CLM_R1"
-  cloud_10m_dir <- file.path(path_tile_s2, "MASKS")
-  s2bands_10m[["Cloud"]] <- file.path(cloud_10m_dir, list.files(cloud_10m_dir, pattern = cloud_10m))
-
-  # get cloud mask 20m
-  cloud_20m <- "CLM_R2"
-  cloud_20m_dir <- file.path(path_tile_s2, "MASKS")
-  s2bands_20m[["Cloud"]] <- file.path(cloud_20m_dir, list.files(cloud_20m_dir, pattern = cloud_20m))
-
-  # return list bands
-  listbands <- list("s2bands_10m" = s2bands_10m,
-                    "s2bands_20m" = s2bands_20m,
-                    "path_tile_s2" = path_tile_s2,
-                    "metadata" = mtdfile)
-  return(listbands)
-}
-
-#" This function check S2 data level:
-#" - L2A: already atmospherically corrected
-#" - L1C: requires atmospheric corrections with sen2cor
-#"
-#" @param prodname character. original name for the S2 image
-#"
-#" @return s2level character. S2 level: L1C or L2A
-#" @export
-get_s2_level <- function(prodname) {
-  prodname <- basename(prodname)
-  if (length(grep(pattern = "L1C_", x = prodname)) == 1) {
-    s2level <- "L1C"
-  } else if (length(grep(pattern = "L2A_", x = prodname)) == 1) {
-    s2level <- "L2A"
-  }
-  return(s2level)
-}
-
-#" This function gets tile from S2 image
-#"
-#" @param prodname character. original name for the S2 image
-#"
-#" @return tilename character
-#" @importFrom tools file_path_sans_ext
-#" @export
-get_tile <- function(prodname) {
-  prodname <- basename(prodname)
-  tilename <- tools::file_path_sans_ext(gsub("_.*", "", gsub(".*_T", "", prodname)))
-  return(tilename)
-}
-
-#" This function gets acquisition date from S2 image
-#"
-#" @param prodname character. original name for the S2 image
-#"
-#" @return dateacq character
-#" @export
-get_date <- function(prodname) {
-  prodname <- basename(prodname)
-  dateacq <- as.Date(gsub("T.*", "", gsub(".*_20", "20", prodname)), format = "%Y%m%d")
-  return(dateacq)
-}
-
-#" download S2 L1C data from Copernicus hub or Google cloud
-#"
-#" @param list_safe safe object. produced with sen2r::s2_list
-#" @param l1c_path character. path for storage of L1C image
-#" @param path_vector path for a vector file
-#" @param time_interval dates. time interval for S2 query
-#" @param googlecloud boolean. set to TRUE if google cloud SDK is installed and
-#" @param forcegoogle boolean. set to TRUE if only google requested
-#" sen2r configured as an alternative hub for S2 download
-#"
-#" @return prodname character. S2 Product name
-#" @importFrom sen2r safe_is_online s2_list s2_download s2_order check_gcloud
-#" @export
-get_s2_l1c_image <- function(list_safe, l1c_path, path_vector, time_interval,
-                             googlecloud = FALSE, forcegoogle = FALSE) {
-  # Check if available from Copernicus hub first
-  copernicus_avail <- sen2r::safe_is_online(list_safe)
-  # if available: download
-  prodname <- attr(list_safe, which = "name")
-  if (file.exists(file.path(l1c_path, prodname))) {
-    message("L1C file already downloaded")
-    message(file.path(l1c_path, prodname))
-  } else {
-    if (copernicus_avail == TRUE && forcegoogle == FALSE) {
-      sen2r::s2_download(list_safe, outdir = l1c_path)
-    } else if (copernicus_avail == FALSE || forcegoogle == TRUE) {
-      # if not available and googlecloud==TRUE
-      if (googlecloud == TRUE) {
-        # check if google cloud SDK available from this computer
-        ggc <- sen2r::check_gcloud()
-        if (ggc == TRUE) {
-          message("downloading from Google cloud")
-          list_safe_ggc <- sen2r::s2_list(spatial_extent = sf::st_read(dsn = path_vector),
-                                          time_interval = time_interval,
-                                          server = "gcloud")
-          prodname <- attr(list_safe_ggc, which = "name")
-          if (file.exists(file.path(l1c_path, prodname))) {
-            message("L1C file already downloaded")
-            message(file.path(l1c_path, prodname))
-          } else {
-            sen2r::s2_download(list_safe_ggc, outdir = l1c_path)
-            # check if QI_DATA exists in DATASTRIP, and create it if not the case
-            datastrip_path <- file.path(l1c_path, prodname, "DATASTRIP")
-            dsdir <- list.dirs(datastrip_path, recursive = FALSE)
-            if (length(match(list.dirs(dsdir, recursive = FALSE, full.names = FALSE), "QI_DATA")) == 0) {
-              dir.create(file.path(dsdir, "QI_DATA"))
-            }
-          }
-        } else if (ggc == FALSE) {
-          message("googlecloud set to TRUE but missing")
-          message("Please install Google cloud SDK")
-          message("https://cloud.google.com/sdk/docs/install")
-          message("and/or set configuration of sen2r following instructions")
-          message("https://www.r-bloggers.com/2021/06/downloading-sentinel-2-archives-from-google-cloud-with-sen2r/")
-        }
-      }
-    }
-    if (copernicus_avail == FALSE && googlecloud == FALSE) {
-      message("S2 image in Long Term Archive (LTA)")
-      message("Ordering image from  LTA")
-      message("This may take 1 day, please run your script later")
-      orders2 <- sen2r::s2_order(list_safe)
-      message("An alternative is possible with Google cloud SDK")
-      message("https://cloud.google.com/sdk/docs/install")
-      message("and/or set configuration of sen2r following instructions")
-      message("https://www.r-bloggers.com/2021/06/downloading-sentinel-2-archives-from-google-cloud-with-sen2r/")
-    }
-  }
-  return(prodname)
-}
-
-#" download S2 L2A data from Copernicus hub or convert L1C to L2A
-#"
-#" @param l2a_path character. path for storage of L2A image
-#" @param spatial_extent path for a vector file
-#" @param dateacq character. date of acquisition
-#" @param deletel1c Boolean. set TRUE to delete L1C images
-#" @param Sen2Cor Boolean. set TRUE to automatically perform atmospheric corrections using sen2Cor
-#" @param googlecloud boolean. set to TRUE if google cloud SDK is installed and
-#" sen2r configured as an alternative hub for S2 download
-#"
-#" @return pathl2a character. Path for L2A image
-#" @importFrom sen2r s2_list s2_download
-#" @importFrom R.utils getAbsolutePath
-
-#" @export
-get_s2_l2a_image <- function(l2a_path, spatial_extent, dateacq,
-                             deletel1c = FALSE, sen2cor = TRUE,
-                             googlecloud = FALSE) {
-
-  # Needs to be updated: define path for L1c data
-  l1c_path <- l2a_path
-  # define time interval
-  time_interval <- as.Date(c(dateacq, dateacq))
-  # get list S2 products corresponding to study area and date of interest using sen2r package
-  if (googlecloud == TRUE) {
-    server <- c("scihub", "gcloud")
-  } else if (googlecloud == FALSE) {
-    server <- "scihub"
-  }
-  list_safe <- sen2r::s2_list(spatial_extent = sf::st_read(dsn = spatial_extent),
-                              time_interval = time_interval,
-                              server = server, availability = "check")
-  # download products
-  sen2r::s2_download(list_safe, outdir = l2a_path)
-  # name all products
-  prodname <- attr(list_safe, which = "name")
-  prodfullpath <- file.path(l2a_path, prodname)
-  if (sen2cor == TRUE) {
-    for (imgname in prodname) {
-      s2level <- get_s2_level(imgname)
-      if (s2level == "L1C") {
-        datepattern <- gsub(pattern = "-", replacement = "", x = dateacq)
-        pathl2a <- s2_from_l1c_to_l2a(prodname = imgname, l1c_path = l2a_path, l2a_path = l2a_path,
-                                      datepattern = datepattern, tmp_path = NULL)
-        if (deletel1c == TRUE) {
-          unlink(x = R.utils::getAbsolutePath(file.path(l1c_path, prodname)),
-                 recursive = TRUE, force = TRUE)
-          # delete from full path and add atmospherically corrected
-          whichimg <- grep(x = prodfullpath, pattern = imgname)
-          dateacq <- get_date(imgname)
-          tilename <- get_tile(imgname)
-          pathl2a <- list.files(path = l2a_path, pattern = tilename, full.names = TRUE)
-          pathl2a <- pathl2a[grep(x = pathl2a, pattern = dateacq)]
-          pathl2a <- pathl2a[grep(x = basename(pathl2a), pattern = "L2A")]
-          prodfullpath[whichimg] <- pathl2a
-        }
-      }
-    }
-  }
-
-  return(prodfullpath)
-}
-
-#" convert image coordinates from index to X-Y
-#"
-#" @param Raster image raster object
-#" @param image_index coordinates corresponding to the raster
-ind2sub <- function(data_raster, image_index) {
-  c <- ((image_index - 1) %% data_raster@ncols) + 1
-  r <- floor((image_index - 1) / data_raster@ncols) + 1
-  my_list <- list("col" = c, "row" = r)
-  return(my_list)
-}
-
-#" mosaicing a set of rasters
-#"
-#" @param list_rasters character. list of paths corresponding to rasters to mosaic
-#" @param dst_mosaic character. path and name of mosaic produced
-#" @param stretch boolean. Set TRUE to get 10% stretching at display for reflectance, mentioned in hdr only
-#"
-#" @return None
-#" @importFrom gdalUtils mosaic_rasters
-#" @importFrom raster hdr raster
-#" @export
-mosaic_rasters <- function(list_rasters, dst_mosaic, stretch = FALSE) {
-
-  # produce mosaic
-  gdalUtils::mosaic_rasters(gdalfile = list_rasters, dst_dataset = dst_mosaic,
-                            separate = FALSE, of = "Ehdr", verbose = TRUE)
-
-  # convert hdr to ENVI format
-  raster::hdr(raster(dst_mosaic), format = "ENVI")
-  # add info to hdr based on initial rasters
-  hdr_init <- read_envi_header(get_hdr_name(list_rasters[1]))
-  hdr <- read_envi_header(get_hdr_name(dst_mosaic))
-  hdr$`band names` <- hdr_init$`band names`
-  hdr$wavelength <- hdr_init$wavelength
-  if (stretch == TRUE) {
-    hdr$`default stretch` <- "0.000000 1000.000000 linear"
-  }
-  hdr$`z plot range` <- NULL
-  hdr$`data ignore value` <- "-Inf"
-  hdr$`sensor type` <- hdr_init$`sensor type`
-  hdr$`coordinate system string` <- read.table(paste(file_path_sans_ext(dst_mosaic), ".prj", sep = ""))
-  write_envi_header(hdr = hdr, hdrpath = get_hdr_name(dst_mosaic))
-  return(invisible())
-}
-
-#" Reads ENVI hdr file
-#"
-#" @param hdrpath Path of the hdr file
-#"
-#" @return list of the content of the hdr file
-#" @export
-read_envi_header <- function(hdrpath) {
-  if (!grepl(".hdr$", hdrpath)) {
-    stop("File extension should be .hdr")
-  }
-  hdr <- readLines(hdrpath)
-  ## check ENVI at beginning of file
-  if (!grepl("ENVI", hdr[1])) {
-    stop("Not an ENVI header (ENVI keyword missing)")
-  } else {
-    hdr <- hdr [-1]
-  }
-  ## remove curly braces and put multi-line key-value-pairs into one line
-  hdr <- gsub("\\{([^}]*)\\}", "\\1", hdr)
-  l <- grep("\\{", hdr)
-  r <- grep("\\}", hdr)
-
-  if (length(l) != length(r)) {
-    stop("Error matching curly braces in header (differing numbers).")
-  }
-
-  if (any(r <= l)) {
-    stop("Mismatch of curly braces in header.")
-  }
-
-  hdr[l] <- sub("\\{", "", hdr[l])
-  hdr[r] <- sub("\\}", "", hdr[r])
-
-  for (i in rev(seq_along(l))) {
-    hdr <- c(
-      hdr [seq_len(l [i] - 1)],
-      paste(hdr [l [i]:r [i]], collapse = "\n"),
-      hdr [-seq_len(r [i])]
-    )
-  }
-
-  ## split key = value constructs into list with keys as names
-  hdr <- sapply(hdr, split_line, "=", USE.NAMES = FALSE)
-  names(hdr) <- tolower(names(hdr))
-
-  ## process numeric values
-  tmp <- names(hdr) %in% c(
-    "samples", "lines", "bands", "header offset", "data type",
-    "byte order", "default bands", "data ignore value",
-    "wavelength", "fwhm", "data gain values"
-  )
-  hdr [tmp] <- lapply(hdr [tmp], function(x) {
-    as.numeric(unlist(strsplit(x, ", ")))
-  })
-
-  return(hdr)
-}
-
-#" This function reads a list of files corresponding to S2 bands
-#" S2 bands are expected to have uniform spatial resolution and footprint
-#" @param s2_bands list. list of S2 bands obtained from get_s2_bands
-#" @param path_vector path for a vector file
-#" @param resampling numeric. resampling factor (default = 1, set to resampling = 2 to convert 20m into 10m resolution)
-#" @param interpolation character. method for resampling. default = "bilinear"
-#"
-#" @return stack_s2 list. contains stack of S2 bands
-#"
-#" @importFrom stars read_stars
-#" @importFrom sf st_bbox st_read st_crop
-#" @export
-
-read_s2bands <- function(s2_bands, path_vector = NULL,
-                         resampling = 1, interpolation = "bilinear") {
-  # get bounding box corresponding to footprint of image or image subset
-  bb_xycoords <- get_bb(path_raster = s2_bands[[1]],
-                        path_vector = path_vector, buffer = 50)
-
-  # prepare reading data for extent defined by bounding box
-  nxoff <- bb_xycoords$UL$col
-  nyoff <- bb_xycoords$UL$row
-  nxsize <- bb_xycoords$UR$col - bb_xycoords$UL$col + 1
-  nysize <- bb_xycoords$LR$row - bb_xycoords$UR$row + 1
-  nbufxsize <- resampling * nxsize
-  nbufysize <- resampling * nysize
-  if (resampling == 1) {
-    interpolation <- "nearest_neighbour"
-  }
-  # write interpolated individual bands in temp directory
-  tmpdir <- tempdir()
-  tmpfile <- list()
-  for (band in names(s2_bands)) {
-    stack_s2_tmp <- stars::read_stars(s2_bands[[band]], along = "band",
-                                      RasterIO = list(nXOff = nxoff, nYOff = nyoff,
-                                                      nXSize = nxsize, nYSize = nysize,
-                                                      nBufXSize = nbufxsize, nBufYSize = nbufysize,
-                                                      resample = interpolation), proxy = FALSE)
-    if (!is.null(path_vector)) {
-      stack_s2_tmp <- sf::st_crop(x = stack_s2_tmp, y = st_bbox(st_read(dsn = path_vector, quiet = TRUE)))
-    }
-    tmpfile[[band]] <- file.path(tmpdir, tools::file_path_sans_ext(basename(s2_bands[[band]])))
-    if (band == "Cloud") {
-      stars::write_stars(obj = stack_s2_tmp, dsn = tmpfile[[band]],
-                         driver =  "ENVI", type = "Byte", overwrite = TRUE)
-    } else {
-      stars::write_stars(obj = stack_s2_tmp, dsn = tmpfile[[band]],
-                         driver =  "ENVI", type = "Int16", overwrite = TRUE)
-    }
-    gc()
-  }
-
-  stack_s2 <- stars::read_stars(tmpfile, along = "band", proxy = TRUE)
-  return(stack_s2)
-}
-
-#" This function reads a raster stack, and gets footprint as pixel coordinates or vector file as input
-#" @param path_raster character. path for raster file
-#" @param path_vector character. path for vector file
-#" @param bbpix list. coordinates of pixels corresponding to a bounding box
-#"
-#" @return starsobj stars object corresponding to raster or raster subset
-#"
-#" @importFrom stars read_stars
-#" @importFrom sf st_bbox st_read st_crop
-#" @export
-read_raster <- function(path_raster, path_vector = NULL, bbpix = NULL) {
-  # get bounding box corresponding to footprint of image or image subset
-  if (is.null(bbpix)) {
-    bb_xycoords <- get_bb(path_raster = path_raster,
-                          path_vector = path_vector, buffer = 0)
-  } else {
-    bb_xycoords <- bbpix
-  }
-  # prepare reading data for extent defined by bounding box
-  nxoff <- bb_xycoords$UL$col
-  nyoff <- bb_xycoords$UL$row
-  nxsize <- bb_xycoords$UR$col - bb_xycoords$UL$col + 1
-  nysize <- bb_xycoords$LR$row - bb_xycoords$UR$row + 1
-  nbufxsize <- nxsize
-  nbufysize <- nysize
-  starsobj <- stars::read_stars(path_raster, along = "band",
-                                RasterIO = list(nXOff = nxoff, nYOff = nyoff,
-                                                  nXSize = nxsize, nYSize = nysize,
-                                                  nBufXSize = nbufxsize, nBufYSize = nbufysize),
-                                proxy = FALSE)
-  return(starsobj)
-}
-
-#" This function reprojects a shapefile and saves reprojected shapefile
-#"
-#" @param path_vector_init character. path for a shapefile to be reprojected
-#" @param newprojection character. projection to be applied to path_vector_init
-#" @param path_vector_reproj character. path for the reprojected shapefile
-#"
-#" @return path_vector character. path of the shapefile
-#" - path_vector_init if the vector did not need reprojection
-#" - path_vector_reproj if the vector needed reprojection
-#"
-#" @importFrom rgdal readOGR writeOGR
-#" @importFrom sp spTransform
-#" @importFrom raster projection
-#" @export
-reproject_shp <- function(path_vector_init, newprojection, path_vector_reproj) {
-
-  dir_vector_init <- dirname(path_vector_init)
-  # shapefile extension
-  fileext <- file_ext(basename(path_vector_init))
-  if (fileext == "shp") {
-    name_vector_init <- file_path_sans_ext(basename(path_vector_init))
-    vector_init_ogr <- rgdal::readOGR(dir_vector_init, name_vector_init, verbose = FALSE)
-  } else if (fileext == "kml") {
-    vector_init_ogr <- rgdal::readOGR(path_vector_init, verbose = FALSE)
-  }
-  vector_init_proj <- raster::projection(vector_init_ogr)
-
-  if (!vector_init_proj == newprojection) {
-    dir_vector_reproj <- dirname(path_vector_reproj)
-    name_vector_reproj <- file_path_sans_ext(basename(path_vector_reproj))
-    vector_reproj <- sp::spTransform(vector_init_ogr, newprojection)
-    rgdal::writeOGR(obj = vector_reproj, dsn = dir_vector_reproj, layer = name_vector_reproj,
-                    driver = "ESRI Shapefile", overwrite_layer = TRUE)
-    path_vector <- path_vector_reproj
-  } else {
-    path_vector <- path_vector_init
-  }
-  return(path_vector)
-}
-
-
-#" perform atmospheric corrections to convert L1C to L2A data with Sen2cor
-#"
-#" @param prodname character. produced with sen2r::s2_list
-#" @param l1c_path character. path of directory where L1C image is stored
-#" @param l2a_path character. path of directory where L2A image is stored
-#" @param datepattern character. pattern corresponding to date of acquisition to identify L2A directory
-#" @param tmp_path character. path of temporary directory where L2A image is stored
-#" sen2r configured as an alternative hub for S2 download
-#"
-#" @return pathl2a character. S2 Product name
-#" @importFrom sen2r safe_is_online s2_list s2_download s2_order
-#" @importFrom R.utils getAbsolutePath
-#"
-#" @export
-s2_from_l1c_to_l2a <- function(prodname, l1c_path, l2a_path, datepattern, tmp_path = NULL) {
-
-  # define path for tmp directory
-  if (is.null(tmp_path)) {
-    tmp_path <- tempdir(check = TRUE)
-  }
-  tmp_prodlist <- prodname
-  # perform Sen2Cor atmospheric corrections
-  binpath <- sen2r::load_binpaths()
-  # 2- open a command prompt and directly run sen2cor with following command line
-  cmd <- paste(binpath$sen2cor,
-               "--output_dir", R.utils::getAbsolutePath(l2a_path),
-               R.utils::getAbsolutePath(file.path(l1c_path, prodname)), sep = " ")
-  system(cmd)
-  pathl2a <- list.files(path = l2a_path, pattern = datepattern, full.names = TRUE)
-
-  return(pathl2a)
-}
-
-#" This function saves cloud masks.
-#" "cloudMask_Binary" is default binary mask with 0 where clouds are detected and 1 for clean pixels
-#" "cloudMask_RAW" is the original cloud layer produced by atmospheric correction algorithm
-#" --> may be useful to refine cloud mask
-#"
-#" @param s2_stars list. stars object containing raster data. Can be produced with function extract_from_s2_l2a
-#" @param cloud_path character.
-#" @param s2source character.
-#" @param footprint character. path for vector file defining footprint of interest in the image
-#" @param saveraw boolean. should the original cloud mask layer be saved?
-#" @param maxchunk numeric. Size of individual chunks to be written (in Mb)
-#"
-#" @return list of cloudmasks (binary mask, and raw mask if required)
-#" @importFrom sf st_read
-#" @importFrom stars write_stars
-#" @importFrom raster raster
-#" @export
-save_cloud_s2 <- function(s2_stars, cloud_path, s2source = "SAFE",
-                          footprint = NULL, saveraw = FALSE, maxchunk = 256) {
-
-  whichcloud <- which(names(s2_stars$attr) == "Cloud")
-  # Save cloud mask
-  if (saveraw == TRUE) {
-    cloudraw <- file.path(cloud_path, "CloudMask_RAW")
-    obj <- stars::read_stars(s2_stars$attr[whichcloud], proxy = TRUE)
-    sizeobj <- dim(obj)[1] * dim(obj)[2] / (1024**2)
-    nbchunks <- ceiling(sizeobj / maxchunk)
-    stars::write_stars(obj,
-                       dsn = cloudraw,
-                       driver =  "ENVI",
-                       type = "Byte",
-                       chunk_size = c(dim(obj)[1], dim(obj)[2] / nbchunks),
-                       progress = TRUE)
-  } else {
-    cloudraw <- NULL
-  }
-  # Save cloud mask as in biodivMapR (0 = clouds, 1 = pixel ok)
-  cloudmask <- stars::read_stars(s2_stars$attr[whichcloud], proxy = FALSE)
-  if (s2source == "SAFE" || s2source == "THEIA") {
-    cloudy <- which(cloudmask[[1]] > 0)
-    sunny <- which(cloudmask[[1]] == 0)
-  } else if (s2source == "LaSRC") {
-    cloudy <- which(is.na(cloudmask[[1]]))
-    sunny <- which(cloudmask[[1]] == 1)
-  }
-
-  cloudmask[[1]][cloudy] <- 0
-  cloudmask[[1]][sunny] <- 1
-  cloudbin <- file.path(cloud_path, "CloudMask_Binary")
-  stars::write_stars(cloudmask, dsn = cloudbin, driver =  "ENVI", type = "Byte", overwrite = TRUE)
-  cloudmasks <- list("BinaryMask" = cloudbin, "RawMask" = cloudraw)
-  # delete temporary file
-  file.remove(s2_stars$attr[whichcloud])
-  if (file.exists(paste(s2_stars$attr[whichcloud], ".hdr", sep = ""))) file.remove(paste(s2_stars$attr[whichcloud], ".hdr", sep = ""))
-  gc()
-  return(cloudmasks)
-}
-
-#" This function saves reflectance files
-#"
-#" @param s2_stars list. stars object containing raster data. Can be produced with function extract_from_s2_l2a
-#" @param refl_path character. path for reflectance file to be stored
-#" @param format character. file format for reflectance data
-#" @param datatype character. data type (integer, float, 16bits, 32bits...)
-#" @param s2sat character. Sentinel-2 mission ("2A" or "2B")
-#" @param tile_s2 character. S2 tile name (2 numbers + 3 letters)
-#" @param dateacq_s2 double. date of acquisition
-#" @param MTD character. path for metadata file
-#" @param MTD_MSI character. path for metadata MSI file
-#" @param mtd_lasrc character. path for metadata LaSRC file
-#" @param maxchunk numeric. Size of individual chunks to be written (in Mb)
-#"
-#" @return None
-#" @importFrom stars write_stars st_apply
-#" @importFrom XML xml
-#" @export
-save_reflectance_s2 <- function(s2_stars, refl_path, format = "ENVI", datatype = "Int16",
-                                s2sat = NULL, tile_s2 = NULL, dateacq_s2 = NULL,
-                                mtd = NULL, mtd_msi = NULL, mtd_lasrc = NULL,
-                                maxchunk = 256) {
-  # identify if S2A or S2B, if possible
-  s2mission <- check_s2mission(s2sat = s2sat, tile_s2 = tile_s2, dateacq_s2 = dateacq_s2)
-
-  # define central wavelength corresponding to each spectral band
-  if (s2mission == "2A") {
-    wl_s2 <- list("B02" = 496.6, "B03" = 560.0, "B04" = 664.5,
-                  "B05" = 703.9, "B06" = 740.2, "B07" = 782.5, "B08" = 835.1,
-                  "B8A" = 864.8, "B11" = 1613.7, "B12" = 2202.4)
-  } else if (s2mission == "2B") {
-    wl_s2 <- list("B02" = 492.1, "B03" = 559.0, "B04" = 665.0,
-                  "B05" = 703.8, "B06" = 739.1, "B07" = 779.7, "B08" = 833.0,
-                  "B8A" = 864.0, "B11" = 1610.4, "B12" = 2185.7)
-  }
-  if (s2mission == "2A") {
-    sensor <- "Sentinel_2A"
-  } else if (s2mission == "2B") {
-    sensor <- "Sentinel_2B"
-  }
-
-  # apply offset when necessary
-  listbands_bis <- c("B2", "B3", "B4", "B5", "B6", "B7", "B8", "B8A", "B11", "B12")
-  if (!is.null(mtd_msi) && is.null(mtd_lasrc)) {
-    # read XML file containing info about geometry of acquisition
-    s2xml <- XML::xmlToList(mtd_msi)
-    xml_offset <- s2xml$General_Info$Product_Image_Characteristics$BOA_ADD_offset_VALUES_LIST
-    bands <- lapply(s2xml$General_Info$Product_Image_Characteristics$Spectral_Information_List, "[[", 4)
-    if (!is.null(xml_offset) && !is.null(bands)) {
-      bandid  <- lapply(bands, "[[", 1)
-      bandname  <- lapply(bands, "[[", 2)
-      offset <- data.frame("bandname" = unlist(bandname),
-                           "bandid" = unlist(bandid),
-                           "offset" = unlist(lapply(xml_offset, "[[", 1)))
-      selbands <- match(listbands_bis, offset$bandname)
-      offset <- offset[selbands, ]
-      boa_quantval <- as.numeric(s2xml$General_Info$Product_Image_Characteristics$QUANTIFICATION_VALUES_LIST$BOA_QUANTIFICATION_VALUE[1])
-    } else {
-      offset <- data.frame("bandname" = listbands_bis,
-                           "bandid" = c(1, 2, 3, 4, 5, 6, 7, 8, 11, 12),
-                           "offset" = 0)
-      boa_quantval <- 10000
-    }
-  } else if (!is.null(mtd_lasrc)) {
-    # read XML file containing info about geometry of acquisition
-    s2xml <- XML::xmlToList(mtd_lasrc)
-    attributes_lasrc <- s2xml$bands[[14]]$.attrs
-    attributes_lasrc_df <- data.frame(attributes_lasrc)
-    if (match("add_offset", rownames(attributes_lasrc_df)) > 0 && match("scale_factor", rownames(attributes_lasrc_df)) > 0) {
-      xml_offset <- as.numeric(attributes_lasrc[["add_offset"]])
-      boa_quantval <- 1 / as.numeric(attributes_lasrc[["scale_factor"]])
-      offset <- data.frame("bandname" = listbands_bis,
-                           "bandid" = c(1, 2, 3, 4, 5, 6, 7, 8, 11, 12),
-                           "offset" = xml_offset)
-    } else {
-      offset <- data.frame("bandname" = listbands_bis,
-                           "bandid" = c(1, 2, 3, 4, 5, 6, 7, 8, 11, 12),
-                           "offset" = 0)
-      boa_quantval <- 10000
-    }
-  } else {
-    offset <- data.frame("bandname" = listbands_bis,
-                         "bandid" = c(1, 2, 3, 4, 5, 6, 7, 8, 11, 12),
-                         "offset" = 0)
-    boa_quantval <- 10000
-  }
-
-  # identify where spectral bands are in the stars object
-  stars_spectral <- list()
-  starsnames <- names(s2_stars$attr)
-  stars_spectral$bandname <- starsnames[which(!starsnames == "Cloud")]
-  stars_spectral$wavelength <- wl_s2[stars_spectral$bandname]
-
-  sortedwl <- names(wl_s2)
-  reorder <- match(sortedwl, stars_spectral$bandname)
-  elim <- which(is.na(reorder))
-  if (length(elim) > 0) {
-    reorder <- reorder[-elim]
-  }
-  pathr <- s2_stars$attr[reorder]
-
-  names(pathr) <- NULL
-  s2_stars2 <- stars::read_stars(pathr, along = "band", proxy = TRUE)
-  stars_spectral$bandname <- stars_spectral$bandname[reorder]
-  stars_spectral$wavelength <- stars_spectral$wavelength[reorder]
-
-  uniqueoffset <- as.numeric(unique(offset$offset))
-  if (length(uniqueoffset) > 1) {
-    message("Warning: BOA offset differs between bands.")
-    message("offset will not be applied to the final S2 reflectance raster")
-    message("check metadata file to identify the offset applied on each band")
-    print(mtd_msi)
-  } else {
-    message("applying offset to reflectance data")
-    if (is.null(mtd_lasrc) || uniqueoffset == 0) {
-      offsets2 <- function(x) (round(x + uniqueoffset) * (10000 / boa_quantval))
-      s2_stars2 <- stars::st_apply(X = s2_stars2, MARGIN = "band", FUN = offsets2)
-    } else {
-      offsets2 <- function(x) (round(10000 * ((x + uniqueoffset * boa_quantval) / boa_quantval)))
-      s2_stars2 <- stars::st_apply(X = s2_stars2, MARGIN = "band", FUN = offsets2)
-    }
-  }
-  write_stack_s2(stars_s2 = s2_stars2, stars_spectral = stars_spectral, refl_path = refl_path,
-                 format = format, datatype = datatype, sensor = sensor, maxchunk = maxchunk)
-  # save metadata file as well if available
-  if (!is.null(mtd)) {
-    if (file.exists(mtd)) {
-      file.copy(from = mtd, to = file.path(dirname(refl_path), basename(mtd)), overwrite = TRUE)
-    }
-  }
-  # save metadata file as well if available
-  if (!is.null(mtd_msi)) {
-    if (file.exists(mtd_msi)) {
-      file.copy(from = mtd_msi, to = file.path(dirname(refl_path), basename(mtd_msi)), overwrite = TRUE)
-    }
-  }
-  # save LaSRC metadata file as well if available
-  if (!is.null(mtd_lasrc)) {
-    if (file.exists(mtd_lasrc)) {
-      file.copy(from = mtd_lasrc, to = file.path(dirname(refl_path), basename(mtd_lasrc)), overwrite = TRUE)
-    }
-  }
-  # delete temporary file
-  for (pathtemp in pathr) {
-    file.remove(pathtemp)
-    if (file.exists(paste(pathtemp, ".hdr", sep = ""))) file.remove(paste(pathtemp, ".hdr", sep = ""))
-  }
-  gc()
-  return(invisible())
-}
-
-#" ENVI functions
-#"
-#" based on https://github.com/cran/hyperSpec/blob/master/R/read.ENVI.R
-#" added wavelength, fwhm, ... to header reading
-#" Title
-#"
-#" @param x character.
-#" @param separator character
-#" @param trim_blank boolean.
-#"
-#" @return list.
-#" @export
-split_line <- function(x, separator, trim_blank = TRUE) {
-  tmp <- regexpr(separator, x)
-  key <- substr(x, 1, tmp - 1)
-  value <- substr(x, tmp + 1, nchar(x))
-  if (trim_blank) {
-    blank_pattern <- "^[[:blank:]]*([^[:blank:]]+.*[^[:blank:]]+)[[:blank:]]*$"
-    key <- sub(blank_pattern, "\\1", key)
-    value <- sub(blank_pattern, "\\1", value)
-  }
-  value <- as.list(value)
-  names(value) <- key
-  return(value)
-}
-
-#" save raster footprint as vector file
-#"
-#" @param path_raster character. path for a raster file
-#" @param path_vector character. path for a vector file
-#" @param driver character. driver for vector
-#"
-#" @return None
-#" @importFrom raster raster extent projection
-#" @importFrom sf st_as_sf st_write
-#" @export
-vectorize_raster_extent <- function(path_raster, path_vector, driver = "ESRI Shapefile") {
-  rast <- raster(path_raster)
-  e <- extent(rast)
-  # coerce to a SpatialPolygons object
-  p <- as(e, "SpatialPolygons")
-  projection(p) <- projection(rast)
-  p <- sf::st_as_sf(p)
-  sf::st_write(obj = p, path_vector, driver = driver)  # create to a shapefile
-  return(invisible())
-}
-
-#" writes ENVI hdr file
-#"
-#" @param hdr content to be written
-#" @param hdrpath Path of the hdr file
-#"
-#" @return None
-#" @importFrom stringr str_count
-#" @export
-write_envi_header <- function(hdr, hdrpath) {
-  h <- lapply(hdr, function(x) {
-    if (length(x) > 1 || (is.character(x) && stringr::str_count(x, "\\w+") > 1)) {
-      x <- paste0("{", paste(x, collapse = ", "), "}")
-    }
-    # convert last numerics
-    x <- as.character(x)
-  })
-  writeLines(c("ENVI", paste(names(hdr), h, sep = " = ")), con = hdrpath)
-  return(invisible())
-}
-
-#" This function writes a raster Stack object into a ENVI raster file
-#"
-#" @param stackobj list. raster stack
-#" @param stackpath character. path where to store the stack
-#" @param bands list. should include "bandname", and if possible "wavelength"
-#" @param datatype character. should be INT2S or FLT4S for example
-#" @param sensor character. Name of the sensor used to acquire the image
-#" @param stretch boolean. Set TRUE to get 10% stretching at display for reflectance, mentioned in hdr only
-#"
-#" @return None
-#" @importFrom utils read.table
-#" @export
-write_rasterstack_envi <- function(stackobj, stackpath, bands, datatype = "INT2S",
-                                   sensor = "Unknown", stretch = FALSE) {
-
-  r <- raster::writeRaster(x = stackobj, filename = stackpath, format = "Ehdr", overwrite = TRUE, datatype = datatype)
-  raster::hdr(r, format = "ENVI")
-  # Edit hdr file to add metadata
-  hdr <- read_envi_header(get_hdr_name(stackpath))
-  hdr$`band names` <- bands$bandname
-  if (length(bands$wavelength) == length(bands$bandname)) {
-    hdr$wavelength <- bands$wavelength
-  } else {
-    hdr$wavelength <- NULL
-  }
-  if (stretch == TRUE) {
-    hdr$`default stretch` <- "0.000000 1000.000000 linear"
-  }
-  hdr$`z plot range` <- NULL
-  hdr$`data ignore value` <- "-Inf"
-  hdr$`coordinate system string` <- read.table(paste(stackpath, ".prj", sep = ""))
-  proj <- strsplit(x = strsplit(x = projection(stackobj), split = " ")[[1]][1], split = "=")[[1]][2]
-  zone <- strsplit(x = strsplit(x = projection(stackobj), split = " ")[[1]][2], split = "=")[[1]][2]
-  datum <- strsplit(x = strsplit(x = projection(stackobj), split = " ")[[1]][3], split = "=")[[1]][2]
-  oldproj <- hdr$`map info`
-  newproj <- gsub(pattern = "projection", replacement = proj, x = oldproj)
-  newproj <- paste(newproj, zone, datum, sep = ", ")
-  hdr$`map info` <- newproj
-  hdr$`sensor type` <- sensor
-  write_envi_header(hdr = hdr, hdrpath = get_hdr_name(stackpath))
-
-  # remove unnecessary files
-  file2remove <- paste(stackpath, ".aux.xml", sep = "")
-  file.remove(file2remove)
-  file2remove <- paste(stackpath, ".prj", sep = "")
-  file.remove(file2remove)
-  file2remove <- paste(stackpath, ".stx", sep = "")
-  file.remove(file2remove)
-  return(invisible())
-}
-
-
-#" This function writes a stars object into a raster file
-#"
-#" @param stars_s2 list. stars object containing raster data. Can be produced with function Crop_n_resample_S2
-#" @param stars_spectral list. band name to be saved in the stack and spectral bands corresponding to the image
-#" @param refl_path character. path where to store the image
-#" @param format character. default = ENVI BSQ. otherwise use gdal drivers
-#" @param datatype character. should be Int16 or Float64 for example
-#" @param sensor character. Name of the sensor used to acquire the image
-#" @param maxchunk numeric. Size of individual chunks to be written (in Mb)
-#"
-#" @return None
-#" @export
-write_stack_s2 <- function(stars_s2, stars_spectral, refl_path, format = "ENVI",
-                           datatype = "Int16", sensor = "Unknown", maxchunk = 256) {
-
-  # write raster file from proxy using chunks
-  sizeobj <- 2 * dim(stars_s2)[1] * dim(stars_s2)[2] * dim(stars_s2)[3] / (1024**2)
-  nbchunks <- ceiling(sizeobj / maxchunk)
-  stars::write_stars(obj = stars_s2,
-                     dsn = refl_path,
-                     driver =  format,
-                     type = datatype,
-                     chunk_size = c(dim(stars_s2)[1], ceiling(dim(stars_s2)[2] / nbchunks)),
-                     progress = TRUE)
-
-  if (format == "ENVI") {
-    adjust_envi_hdr(dsn = refl_path, bands = stars_spectral,
-                    sensor = sensor, stretch = TRUE)
-  }
-  return(invisible())
-}

--- a/alpha_beta.r	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,142 +0,0 @@
-#Rscript
-
-###########################################
-##    Mapping alpha and beta diversity   ##
-###########################################
-
-#####Packages : stars
-#               utils
-#               biodivmapr
-#               raster
-#               sf
-#               mapview
-#               leafpop
-#               RColorBrewer
-#               labdsv
-#               rgdal
-#               ggplot2
-#               gridExtra
-
-#####Load arguments
-
-args <- commandArgs(trailingOnly = TRUE)
-
-#####Import the S2 data
-
-if (length(args) < 1) {
-    stop("This tool needs at least 1 argument")
-}else {
-    data_raster <- args[1]
-    rasterheader <- args[2]
-    data <- args[3]
-    # type of PCA:
-    # PCA: no rescaling of the data
-    # SPCA: rescaling of the data
-    typepca <- as.character(args[4])
-    alpha <- as.logical(args[5])
-    beta <- as.logical(args[6])
-    funct <- as.logical(args[7])
-    all <- as.logical(args[8])
-    nbcpu <- as.integer(args[9])
-    source(args[10])
-}
-
-################################################################################
-##              DEFINE PARAMETERS FOR DATASET TO BE PROCESSED                 ##
-################################################################################
-if (data_raster == "") {
-  #Create a directory where to unzip your folder of data
-  dir.create("data_dir")
-  unzip(data, exdir = "data_dir")
-  # Path to raster
-  data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl")
-  input_image_file <- file.path("data_dir/results/Reflectance", data_raster[1])
-  input_header_file <- file.path("data_dir/results/Reflectance", data_raster[2])
-
-} else {
-  input_image_file <- file.path(getwd(), data_raster, fsep = "/")
-  input_header_file <- file.path(getwd(), rasterheader, fsep = "/")
-}
-
-################################################################################
-##                              PROCESS IMAGE                                 ##
-################################################################################
-# 1- Filter data in order to discard non vegetated / shaded / cloudy pixels
-
-print("PERFORM PCA ON RASTER")
-pca_output <- biodivMapR::perform_PCA(Input_Image_File = input_image_file, Input_Mask_File = input_mask_file,
-                          Output_Dir = output_dir, TypePCA = typepca, FilterPCA = filterpca, nbCPU = nbcpu, MaxRAM = maxram)
-
-pca_files <- pca_output$PCA_Files
-pix_per_partition <- pca_output$Pix_Per_Partition
-nb_partitions <- pca_output$nb_partitions
-# path for the updated mask
-input_mask_file <- pca_output$MaskPath
-
-
-selected_pcs <- seq(1, dim(raster::stack(input_image_file))[3])
-
-selected_pcs <- all(selected_pcs)
-################################################################################
-##                      MAP ALPHA AND BETA DIVERSITY                          ##
-################################################################################
-print("MAP SPECTRAL SPECIES")
-
-kmeans_info <- biodivMapR::map_spectral_species(Input_Image_File = input_image_file, Output_Dir = output_dir, PCA_Files = pca_files, Input_Mask_File = input_mask_file, SelectedPCs = selected_pcs, Pix_Per_Partition = pix_per_partition, nb_partitions = nb_partitions, TypePCA = typepca, nbCPU = nbcpu, MaxRAM = maxram, nbclusters = nbclusters)
-
-image_name <- tools::file_path_sans_ext(basename(input_image_file))
-if (alpha == TRUE || beta == TRUE || all == TRUE) {
-## alpha
-  print("MAP ALPHA DIVERSITY")
-  index_alpha <- c("Shannon")
-  alpha_div <- biodivMapR::map_alpha_div(Input_Image_File = input_image_file, Output_Dir = output_dir, TypePCA = typepca, window_size = window_size, nbCPU = nbcpu, MaxRAM = maxram, Index_Alpha = index_alpha, nbclusters = nbclusters, FullRes = TRUE, LowRes = FALSE, MapSTD = FALSE)
-
-  alpha_zip <- file.path(output_dir, image_name, typepca, "ALPHA", "Shannon_10_Fullres.zip")
-  alpha_path <- file.path(output_dir, image_name, typepca, "ALPHA")
-  unzip(alpha_zip, exdir = alpha_path)
-  alpha_path <- file.path(output_dir, image_name, typepca, "ALPHA", "Shannon_10_Fullres")
-  alpha_raster <- raster::raster(alpha_path)
-  get_alpha <- convert_raster(alpha_raster)
-
- if (alpha == TRUE || all == TRUE) {
-   colnames(get_alpha) <- c("Alpha", "longitude", "latitude")
-   plot_indices(get_alpha, titre = "Alpha")
-
-   write.table(get_alpha, file = "alpha.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
-}
-  if (beta == TRUE || all == TRUE) {
-## beta
-  print("MAP BETA DIVERSITY")
-  beta_div <- biodivMapR::map_beta_div(Input_Image_File = input_image_file, Output_Dir = output_dir, TypePCA = typepca, window_size = window_size, nb_partitions = nb_partitions, nbCPU = nbcpu, MaxRAM = maxram, nbclusters = nbclusters)
-
-  beta_path <- file.path(output_dir, image_name, typepca, "BETA", "BetaDiversity_BCdiss_PCO_10")
-  beta_raster <- raster::raster(beta_path)
-  get_beta <- convert_raster(beta_raster)
-
-  colnames(get_beta) <- c("Beta", "longitude", "latitude")
-  plot_indices(get_beta, titre = "Beta")
-
-  write.table(get_beta, file = "beta.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
-  }
-}
-
-
-################################################################################
-##          COMPUTE ALPHA AND BETA DIVERSITY FROM FIELD PLOTS                 ##
-################################################################################
-
-if (funct == TRUE || all == TRUE) {
-  mapper <- biodivMapR::map_functional_div(Original_Image_File = input_image_file, Functional_File = pca_files,  Selected_Features = selected_pcs, Output_Dir = output_dir, window_size = window_size, nbCPU = nbcpu, MaxRAM = maxram, TypePCA = typepca, FullRes = TRUE, LowRes = FALSE, MapSTD = FALSE)
-
-  funct_zip <- file.path(output_dir, image_name, typepca, "FUNCTIONAL", "FunctionalDiversity_Map_MeanFilter_Fullres.zip")
-  funct_path <- file.path(output_dir, image_name, typepca, "FUNCTIONAL")
-  unzip(funct_zip, exdir = funct_path)
-  funct_path <- file.path(output_dir, image_name, typepca, "FUNCTIONAL", "FunctionalDiversity_Map_MeanFilter_Fullres")
-  funct_raster <- raster::raster(funct_path)
-  get_funct <- convert_raster(funct_raster)
-
-  colnames(get_funct) <- c("Functionnal", "longitude", "latitude")
-  plot_indices(get_funct, titre = "Functionnal")
-
-  write.table(get_funct, file = "Functionnal.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
-}

--- a/alpha_beta.xml	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,153 +0,0 @@
-<tool id="srs_diversity_maps" name="Map diversity" version="@VERSION@" profile="20.01">
-    <description>from remote sensing data</description>
-    <macros>
-        <import>macro.xml</import>
-    </macros>
-    <expand macro="SRS_requirements">
-    </expand>
-    <command detect_errors="exit_code"><![CDATA[ 
-        #import re 
-        #if $method.origin == 'envi_bil': 
-          #set input_raster = $method.input_raster
-          #set input_raster_identifier = re.sub('[^\s\w\-]', '_', str($input_raster.element_identifier)) 
-          #set input_header = $method.input_header
-          #set input_header_identifier = re.sub('[^\s\w\-]+[^.hdr]', '_', str($input_header.element_identifier)) 
-          cp '${input_raster}' '${input_raster_identifier}' &&
-          cp '${input_header}' '${input_header_identifier}' &&
-        #end if
-        Rscript
-            '$__tool_directory__/alpha_beta.r'
-             #if $method.origin == 'envi_bil':
-              '$input_raster_identifier' 
-              '$input_header_identifier'
-              ''
-            #else:
-              ''
-              ''
-              '$method.input'
-            #end if
-            '$typepca'
-            #if $type == 'alpha':
-              'TRUE'
-              'FALSE'
-              'FALSE'
-              'FALSE'
-            #else if $type == 'beta':
-              'FALSE'
-              'TRUE'
-              'FALSE'
-              'FALSE'
-            #else if $type == 'funct':
-              'FALSE'
-              'FALSE'
-              'TRUE'
-              'FALSE'
-            #else:
-              'FALSE'
-              'FALSE'
-              'FALSE'
-              'TRUE'
-            #end if
-            \${GALAXY_SLOTS:-1}
-            '$__tool_directory__/functions.r'
-            '$output_alpha'
-            '$output_beta'
-            '$output_funct'
-            '$plots'
-        ]]>
-    </command>
-    <inputs>
-        <conditional name="method">
-            <param name="origin" type="select" label="In which format are your data ?">
-                <option value="zipper">The data you are using are in a zip folder Reflectance</option>
-                <option value="envi_bil">Your already have the files in ENVI BIL format</option>
-            </param>
-            <when value="zipper">
-                <param name="input" type="data" format="zip" multiple="true" label="Input data"/>
-            </when>
-            <when value="envi_bil">
-                <param name="input_raster" type="data" format="bil" label="Input raster" help="It can be the raw data in bil or the PCA raster layer in bil"/>
-                <param name="input_header" type="data" format="hdr" label="Input header"/>
-            </when>
-        </conditional>
-        <param name="typepca" type="select" label="Do you want to do a PCA or a SPCA ?" display="radio" help="If you choose PCA there is no rescaling of the data as oppposed as if you choose SPCA">
-            <option value="SPCA">SPCA</option>
-            <option value="PCA">PCA</option>
-        </param>
-        <param name="type" type="select" label="Alpha, beta, functional diversity and comparison plot and map" display="radio">
-            <option value="alpha">Alpha diversity map</option>
-            <option value="beta">Beta diversity map</option>
-            <option value="funct">Functional diversity map</option>
-            <option value="all">All of the above</option>
-        </param>
-    </inputs>
-    <outputs>
-        <data name="output_alpha" from_work_dir="alpha.tabular" format="tabular" label="Alpha diversity">
-            <filter> type == 'alpha' or  type == 'all'</filter>
-        </data>
-        <data name="output_beta" from_work_dir="beta.tabular" format="tabular" label="Beta diversity">
-            <filter> type == 'beta' or  type == 'all'</filter>
-        </data>
-        <data name="output_funct" from_work_dir="Functionnal.tabular" format="tabular" label="Functionnal diversity">
-            <filter> type == 'funct' or  type == 'all'</filter>
-        </data>
-        <collection type="list" name="plots" label="${type} plot">
-            <discover_datasets pattern="(?P&lt;designation&gt;.+)\.png" visible="false" format="png"/>
-        </collection>
-    </outputs>
-    <tests>
-        <test>
-            <param name="origin" value="envi_bil"/>
-            <param name="input_raster" value="S2A_Subset"/>
-            <param name="input_header" value="S2A_Subset.hdr"/>
-            <param name="type" value="alpha"/>
-            <output name="output_alpha">
-                <assert_contents>
-                    <has_n_columns n="3"/>
-                </assert_contents>
-            </output>
-            <output_collection name="plots" type="list" count="1"/>
-        </test>
-    </tests>
-    <help><![CDATA[
-========================================================================
-Process satellite remote sensing data to produce biodiversity indicators
-========================================================================
-
-
-**What it does**
-
-Féret and Asner (2014) developed a method for **tropical forest** diversity mapping based on very high spatial resolution airborne imaging spectroscopy.
-
-The goal of this tool using the package biodivMapR is to produce (spectral) diversity maps based on (optical) images.
-
-**Input description**
-
-It expects an image file as input, with a specific data format. ENVI HDR image with BIL interleave required.
-The image is an ENVI raster including :
-
-- A binary file (which has no extension here).
-
-- A header file (with .hdr extension).
-
-The header file is a text file including all necessary metadata which can be read with a text editor. It includes image dimensions, projection, and the name and central wavelength for each spectral band.
-
-In order to get such input we advise to use the tool preprocessing sentinel 2 data. 
-
-+--------------+----------+
-|      BIL     | ENVI HDR |
-+==============+==========+
-| raster stack | Metadata |
-+--------------+----------+
-|      ...     |    ...   |
-+--------------+----------+
-
-**Output**
-
-- Three tabulars : alpha, beta, functionnal each of them with 3 colomns latitude, longitude and the indice.
-
-- Three png graph for each indice 
-
-    ]]></help>
-    <expand macro="SRS_BDMRref"/>
-</tool>

--- a/biodiv_indices_global.r	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,112 +0,0 @@
-#Rscript
-
-###########################################
-##    Mapping  biodiversity indicators   ##
-###########################################
-
-#####Packages :  raster,
-#                rgdal,
-#                sp,
-#                rasterdiv,
-#                ggplot2,
-
-#####Load arguments
-
-args <- commandArgs(trailingOnly = TRUE)
-
-#####Import the S2 data
-
-if (length(args) < 1) {
-    stop("This tool needs at least 1 argument")
-}else {
-    data_raster <- args[1]
-    data_header <- args[2]
-    data <- args[3]
-    alpha <- as.numeric(args[4])
-    source(args[5])
-
-}
-
-########################################################################
-##                   COMPUTE BIODIVERSITY INDICES                     ##
-########################################################################
-
-if (data_raster == "") {
-  #Create a directory where to unzip your folder of data
-  dir.create("data_dir")
-  unzip(data, exdir = "data_dir")
-  # Path to raster
-  data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl")
-  data_raster <- file.path("data_dir/results/Reflectance", data_raster[1])
-}
-
-# Read raster
-copndvi <- raster::raster(data_raster)
-
-copndvilr <- raster::reclassify(copndvi, cbind(252, 255, NA), right = TRUE)
-
-#Resample using raster::aggregate and a linear factor of 10
-copndvilr <- raster::aggregate(copndvilr, fact = 20)
-#Set float numbers as integers to further speed up the calculation
-storage.mode(copndvilr[]) <- "integer"
-
-# Convert raster to SpatialPointsDataFrame
-r_pts <- convert_raster(copndvilr)
-
-#Shannon's Diversity
-sha <- rasterdiv::Shannon(copndvilr, window = 9, np = 1)
-sha_df <- data.frame(raster::rasterToPoints(sha, spatial = TRUE))
-sha_name <- "Shannon"
-r_pts[, sha_name] <- sha_df[, 1]
-
-#Renyi's Index
-ren <- rasterdiv::Renyi(copndvilr, window = 9, alpha = alpha, np = 1)
-ren_df <- data.frame(raster::rasterToPoints(ren[[1]]))
-ren_name <- "Renyi"
-r_pts[, ren_name] <- ren_df[, 3]
-
-#Berger-Parker's Index
-ber <- rasterdiv::BergerParker(copndvilr, window = 9, np = 1)
-ber_df <- data.frame(raster::rasterToPoints(ber, spatial = TRUE))
-ber_name <- "Berger-Parker"
-r_pts[, ber_name] <- ber_df[, 1]
-
-#Pielou's Evenness
-pie <- rasterdiv::Pielou(copndvilr, window = 9, np = 1)
-pie_df <- data.frame(raster::rasterToPoints(pie))
-if (length(pie_df[, 1]) == length(r_pts[, 1])) {
-  pie_name <- "Pielou"
-  r_pts[, pie_name] <- pie_df[, 1]
-}
-
-#Hill's numbers
-hil <- rasterdiv::Hill(copndvilr, window = 9, alpha = alpha, np = 1)
-hil_df <- data.frame(raster::rasterToPoints(hil[[1]]))
-hil_name <- "Hill"
-r_pts[, hil_name] <- hil_df[, 3]
-
-#Parametric Rao's quadratic entropy with alpha ranging from 1 to 5
-prao <- rasterdiv::paRao(copndvilr, window = 9, alpha = alpha, dist_m = "euclidean", np = 1)
-prao_df <- data.frame(raster::rasterToPoints(prao$window.9[[1]]))
-prao_name <- "Prao"
-r_pts[, prao_name] <- prao_df[, 3]
-
-#Cumulative Residual Entropy
-cre <- rasterdiv::CRE(copndvilr, window = 9, np = 1)
-cre_df <- data.frame(raster::rasterToPoints(cre))
-if (length(cre_df[, 1]) == length(r_pts[, 1])) {
-  cre_name <- "CRE"
-  r_pts[, cre_name] <- cre_df[, 1]
-}
-
-if (length(cre_df[, 1]) == length(r_pts[, 1]) || length(pie_df[, 1]) == length(r_pts[, 1])) {
-list_indice <- list("Shannon", "Renyi", "Berger-Parker", "Pielou", "Hill", "Prao", "CRE")
-} else {
-list_indice <- list("Shannon", "Renyi", "Berger-Parker", "Hill", "Prao")
-}
-## Plotting all the graph and writing a tabular
-for (indice in list_indice) {
-  plot_indices(data = r_pts, titre = indice)
-}
-
-write.table(r_pts, file = "BiodivIndex.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)

--- a/biodiv_indices_global.xml	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,123 +0,0 @@
-<tool id="srs_global_indices" name="Compute biodiversity indices" version="@VERSION@" profile = "20.01">
-    <description>from remote sensing data</description>
-    <macros>
-        <import>macro.xml</import>
-    </macros>
-    <requirements>
-        <requirement type="package" version="4.2.2">r-base</requirement>
-        <requirement type="package" version="3.5_21">r-raster</requirement>
-        <requirement type="package" version="1.5_32">r-rgdal</requirement>
-        <requirement type="package" version="1.5_1">r-sp</requirement>
-        <requirement type="package" version="0.3.1">r-rasterdiv</requirement>
-        <requirement type="package" version="3.4.0">r-ggplot2</requirement>
-        <!--requirement type="package" version="1.5_21">r-terra</requirement-->
-    </requirements>
-    <command detect_errors="exit_code"><![CDATA[ 
-        #import re 
-        #if $method.type == 'envi_bil': 
-          #set input_raster = $method.input_raster
-          #set input_raster_identifier = re.sub('[^\s\w\-]', '_', str($input_raster.element_identifier)) 
-          #set input_header = $method.input_header
-          #set input_header_identifier = re.sub('[^\s\w\-]+[^.hdr]', '_', str($input_header.element_identifier)) 
-          cp '${input_raster}' '${input_raster_identifier}' &&
-          cp '${input_header}' '${input_header_identifier}' &&
-        #end if
-        Rscript
-            '$__tool_directory__/biodiv_indices_global.r'
-             #if $method.type == 'envi_bil':
-              '$input_raster_identifier' 
-              '$input_header_identifier'
-              ''
-            #else:
-              ''
-              ''
-              '$method.input'
-            #end if
-            '$alpha'
-            '$__tool_directory__/functions.r'
-            '$output_indices'
-            '$plots_png'
-        ]]>
-    </command>
-    <inputs>
-        <conditional name="method">
-            <param name="type" type="select" label="In which format are your data ?">
-                <option value="zipper">The data you are using are in a zip folder Reflectance</option>
-                <option value="envi_bil">Your already have the files in ENVI BIL format</option>
-            </param>
-            <when value="zipper">
-                <param name="input" type="data" format="zip" multiple="true" label="Input data"/>
-            </when>
-            <when value="envi_bil">
-                <param name="input_raster" type="data" format="bil" label="Input raster" help="It can be the raw data in bil or the PCA raster layer in bil"/>
-                <param name="input_header" type="data" format="hdr" label="Input header"/>
-            </when>
-        </conditional>
-        <param name="alpha" type="text" label="Write a number of the value of alpha you want to calculate Renyi, Prao and Hill indicators (between 0 and 5)"/>
-    </inputs>
-    <outputs>
-        <data name="output_indices" from_work_dir="BiodivIndex.tabular" format="tabular" label="Biodiversity indices tabular">
-        </data>
-        <collection type="list" name="plots_png" label="Biodiversity indices plots">
-            <discover_datasets pattern="(?P&lt;designation&gt;.+)\.png" visible="false" format="png"/>
-        </collection>
-    </outputs>
-    <tests>
-        <test>
-            <param name="type" value="envi_bil"/>
-            <param name="input_raster" value="S2A_Subset"/>
-            <param name="input_header" value="S2A_Subset.hdr"/>
-            <param name="alpha" value="1"/>
-            <output name="output_indices">
-                <assert_contents>
-                    <has_n_lines n="2501"/>
-                </assert_contents>
-            </output>
-            <output_collection name="plots_png" type="list" count="7"/>
-        </test>
-    </tests>
-    <help><![CDATA[
-=========================================================================
-Computes global biodiversity indices from satellite remote sensing data 
-=========================================================================
-
-
-**What it does**
-
-This tool aims to provide functions to apply Information Theory based diversity indexes on RasterLayer such as Shannon's entropy or Cumulative Residual Entropy (CRE).
-
-
-**Input description**
-
-It expects an image file as input, with a specific data format. ENVI HDR image with BIL interleave required.
-The image is an ENVI raster including :
-
-- A binary file (which has no extension here).
-
-- A header file (with .hdr extension).
-
-The header file is a text file including all necessary metadata which can be read with a text editor. It includes image dimensions, projection, and the name and central wavelength for each spectral band.
-
-In order to get such input we advise to use the tool preprocessing sentinel 2 data. If you did so you can directly enter the "Reflectance" output from this tool and thus select the otpion "The data you are using are in a zip folder Reflectance". 
-
-⚠️ If you do not use this Reflectance folder make sure that your data are respectively in bil and hdr format in the datatypes.
-
-- A number for the alpha indice which used to calculate the following indicators : Renyi, Hill and Prao.
-
-+--------------+----------+---------------+
-|      BIL     | ENVI HDR |  Number alpha |
-+==============+==========+===============+
-| raster stack | Metadata |       1       |
-+--------------+----------+---------------+
-|      ...     |    ...   |      ...      |
-+--------------+----------+---------------+
-
-**Output**
-
-- One tabular with 9 columns longitude, latitude columns and one for each indices.
-
-- Seven png graph one for each indices.
-
-    ]]>    </help>
-        <expand macro="SRS_rasterdivref"/>
-</tool>

--- a/comparison_div.r	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,190 +0,0 @@
-#Rscript
-
-###########################################
-##    Mapping alpha and beta diversity   ##
-###########################################
-
-#####Packages : stars
-#               utils
-#               biodivmapr
-#               raster
-#               sf
-#               mapview
-#               leafpop
-#               RColorBrewer
-#               labdsv
-#               rgdal
-#               ggplot2
-#               gridExtra
-##remotes::install_github("jbferet/biodivMapR")
-#####Load arguments
-
-args <- commandArgs(trailingOnly = TRUE)
-
-#####Import the S2 data
-
-if (length(args) < 1) {
-    stop("This tool needs at least 1 argument")
-}else {
-    data_raster <- args[1]
-    rasterheader <- args[2]
-    data <- args[3]
-    plots_zip <- args[4]
-    choice <- as.character(args[5])
-    source(args[6])
-    # type of PCA:
-    # PCA: no rescaling of the data
-    # SPCA: rescaling of the data
-    typepca <- as.character(args[7])
-}
-
-################################################################################
-##              DEFINE PARAMETERS FOR DATASET TO BE PROCESSED                 ##
-################################################################################
-if (data_raster == "") {
-  #Create a directory where to unzip your folder of data
-  dir.create("data_dir")
-  unzip(data, exdir = "data_dir")
-  # Path to raster
-  data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl")
-  input_image_file <- file.path("data_dir/results/Reflectance", data_raster[1])
-  input_header_file <- file.path("data_dir/results/Reflectance", data_raster[2])
-
-} else {
-  input_image_file <- file.path(getwd(), data_raster, fsep = "/")
-  input_header_file <- file.path(getwd(), rasterheader, fsep = "/")
-}
-
-################################################################################
-##                              PROCESS IMAGE                                 ##
-################################################################################
-# 1- Filter data in order to discard non vegetated / shaded / cloudy pixels
-
-print("PERFORM PCA ON RASTER")
-pca_output <- biodivMapR::perform_PCA(Input_Image_File = input_image_file, Input_Mask_File = input_mask_file,
-                          Output_Dir = output_dir, TypePCA = typepca, FilterPCA = filterpca, nbCPU = nbcpu, MaxRAM = maxram)
-
-pca_files <- pca_output$PCA_Files
-pix_per_partition <- pca_output$Pix_Per_Partition
-nb_partitions <- pca_output$nb_partitions
-# path for the updated mask
-input_mask_file <- pca_output$MaskPath
-
-# 3- Select principal components from the PCA raster
-# Select components from the PCA/SPCA/MNF raster
-sel_compo <- c("1\n", "2\n", "3\n", "4\n", "5\n", "6\n", "7\n", "8")
-image_name <- tools::file_path_sans_ext(basename(input_image_file))
-output_dir_full <- file.path(output_dir, image_name, typepca, "PCA")
-
-write.table(sel_compo, paste0(output_dir_full, "/Selected_Components.txt"))
-sel_pc <-  file.path(output_dir_full, "Selected_Components.txt")
-
-
-################################################################################
-##                      MAP ALPHA AND BETA DIVERSITY                          ##
-################################################################################
-print("MAP SPECTRAL SPECIES")
-
-kmeans_info <- biodivMapR::map_spectral_species(Input_Image_File = input_image_file, Output_Dir = output_dir, PCA_Files = pca_files, Input_Mask_File = input_mask_file, Pix_Per_Partition = pix_per_partition, nb_partitions = nb_partitions, nbCPU = nbcpu, MaxRAM = maxram, nbclusters = nbclusters, TypePCA = typepca)
-
-################################################################################
-##          COMPUTE ALPHA AND BETA DIVERSITY FROM FIELD PLOTS                 ##
-################################################################################
-## read selected features from dimensionality reduction
-
-## path for selected components
-
-# location of the directory where shapefiles used for validation are saved
-dir.create("VectorDir")
-unzip(plots_zip, exdir = "VectorDir")
-
-# list vector data
-path_vector <- biodivMapR::list_shp("VectorDir")
-name_vector <- tools::file_path_sans_ext(basename(path_vector))
-
-# location of the spectral species raster needed for validation
-path_spectralspecies <- kmeans_info$SpectralSpecies
-# get diversity indicators corresponding to shapefiles (no partitioning of spectral dibversity based on field plots so far...)
-biodiv_indicators <- biodivMapR::diversity_from_plots(Raster_SpectralSpecies = path_spectralspecies, Plots = path_vector, nbclusters = nbclusters, Raster_Functional = pca_files, Selected_Features = FALSE)
-
-shannon_rs <- c(biodiv_indicators$Shannon)[[1]]
-fric <- c(biodiv_indicators$FunctionalDiversity$FRic)
-feve <- c(biodiv_indicators$FunctionalDiversity$FEve)
-fdiv <- c(biodiv_indicators$FunctionalDiversity$FDiv)
-# if no name for plots
-biodiv_indicators$Name_Plot <- seq(1, length(biodiv_indicators$Shannon[[1]]), by = 1)
-
-
-####################################################
-# write RS indicators                              #
-####################################################
-# write a table for Shannon index
-
-# write a table for all spectral diversity indices corresponding to alpha diversity
-results <- data.frame(name_vector, biodiv_indicators$Richness, biodiv_indicators$Fisher,
-                      biodiv_indicators$Shannon, biodiv_indicators$Simpson,
-                      biodiv_indicators$FunctionalDiversity$FRic,
-                      biodiv_indicators$FunctionalDiversity$FEve,
-                      biodiv_indicators$FunctionalDiversity$FDiv)
-
-names(results) <- c("ID_Plot", "Species_Richness", "Fisher", "Shannon", "Simpson", "fric", "feve", "fdiv")
-write.table(results, file = "Diversity.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
-
-if (choice == "Y") {
-# write a table for Bray Curtis dissimilarity
-bc_mean <- biodiv_indicators$BCdiss
-bray_curtis <- data.frame(name_vector, bc_mean)
-colnames(bray_curtis) <- c("ID_Plot", bray_curtis[, 1])
-write.table(bray_curtis, file = "BrayCurtis.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
-
-####################################################
-# illustrate results
-####################################################
-# apply ordination using PCoA (same as done for map_beta_div)
-
-mat_bc_dist <- as.dist(bc_mean, diag = FALSE, upper = FALSE)
-betapco <- labdsv::pco(mat_bc_dist, k = 3)
-
-# assign a type of vegetation to each plot, assuming that the type of vegetation
-# is defined by the name of the shapefile
-
-nbsamples <- shpname <- c()
-for (i in 1:length(path_vector)) {
-  shp <- path_vector[i]
-  nbsamples[i] <- length(rgdal::readOGR(shp, verbose = FALSE))
-  shpname[i] <- tools::file_path_sans_ext(basename(shp))
-}
-
-type_vegetation <- c()
-for (i in 1: length(nbsamples)) {
-  for (j in 1:nbsamples[i]) {
-    type_vegetation <- c(type_vegetation, shpname[i])
-  }
-}
-
-#data frame including a selection of alpha diversity metrics and beta diversity expressed as coordinates in the PCoA space
-results <- data.frame("vgtype" = type_vegetation, "pco1" = betapco$points[, 1], "pco2" = betapco$points[, 2], "pco3" = betapco$points[, 3], "shannon" = shannon_rs, "fric" = fric, "feve" = feve, "fdiv" = fdiv)
-
-#plot field data in the PCoA space, with size corresponding to shannon index
-g1 <- ggplot2::ggplot(results, ggplot2::aes(x = pco1, y = pco2, color = vgtype, size = shannon)) + ggplot2::geom_point(alpha = 0.6) + ggplot2::scale_color_manual(values = c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
-
-g2 <- ggplot2::ggplot(results, ggplot2::aes(x = pco1, y = pco3, color = vgtype, size = shannon)) + ggplot2::geom_point(alpha = 0.6) + ggplot2::scale_color_manual(values = c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
-
-g3 <- ggplot2::ggplot(results, ggplot2::aes(x = pco2, y = pco3, color = vgtype, size = shannon)) + ggplot2::geom_point(alpha = 0.6) + ggplot2::scale_color_manual(values = c("#e6140a", "#e6d214", "#e68214", "#145ae6"))
-
-#extract legend
-get_legend <- function(a_gplot) {
-    tmp <- ggplot2::ggplot_gtable(ggplot2::ggplot_build(a_gplot))
-    leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
-    legend <- tmp$grobs[[leg]]
-    return(legend)
-}
-
-legend <- get_legend(g3)
-gall <- gridExtra::grid.arrange(gridExtra::arrangeGrob(g1 + ggplot2::theme(legend.position = "none"), g2 + ggplot2::theme(legend.position = "none"), g3 + ggplot2::theme(legend.position = "none"), nrow = 1), legend, nrow = 2, heights = c(3, 2))
-
-
-filename <- ggplot2::ggsave("BetaDiversity_PcoA1_vs_PcoA2_vs_PcoA3.png", gall, scale = 0.65, width = 12, height = 9, units = "in", dpi = 200, limitsize = TRUE)
-
-filename
-}

--- a/comparison_div.xml	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,138 +0,0 @@
-<tool id="srs_process_data" name="Compare diversity indicators" version="@VERSION@" profile="20.01">
-    <description>with remote sensing data</description>
-    <macros>
-        <import>macro.xml</import>
-    </macros>
-    <expand macro="SRS_requirements"/>
-    <command detect_errors="exit_code"><![CDATA[
-        #import re 
-        #if $method.type == 'envi_bil': 
-          #set input_raster = $method.input_raster
-          #set input_raster_identifier = re.sub('[^\s\w\-]', '_', str($input_raster.element_identifier)) 
-          #set input_header = $method.input_header
-          #set input_header_identifier = re.sub('[^\s\w\-]+[^.hdr]', '_', str($input_header.element_identifier)) 
-          cp '${input_raster}' '${input_raster_identifier}' &&
-          cp '${input_header}' '${input_header_identifier}' &&
-        #end if
-        Rscript
-            '$__tool_directory__/comparison_div.r'
-             #if $method.type == 'envi_bil':
-              '$input_raster_identifier' 
-              '$input_header_identifier'
-              ''
-            #else:
-              ''
-              ''
-              '$method.input'
-            #end if
-            '$input_plot'
-            '$choice'
-            '$__tool_directory__/functions.r'
-            '$typepca'
-            '$output_div'
-            '$output_beta'
-            '$plots'
-        ]]>
-    </command>
-    <inputs>
-        <conditional name="method">
-            <param name="type" type="select" label="In which format are your data ?">
-                <option value="zipper">The data you are using are in a zip folder Reflectance</option>
-                <option value="envi_bil">Your already have the files in ENVI BIL format</option>
-            </param>
-            <when value="zipper">
-                <param name="input" type="data" format="zip" multiple="true" label="Input data"/>
-            </when>
-            <when value="envi_bil">
-                <param name="input_raster" type="data" format="bil" label="Input raster" help="It can be the raw data in bil or the PCA raster layer in bil"/>
-                <param name="input_header" type="data" format="hdr" label="Input header"/>
-            </when>
-        </conditional>
-        <param name="input_plot" type="data" format="data" label="Plots folder zip"/>
-        <param name="choice" type="select" label="Do you want to compute beta diversity (needs mutliple locations) ?" display="radio">
-            <option value="Y">Yes</option>
-            <option value="N">No</option>
-        </param>
-        <param name="typepca" type="select" label="Do you want to do a PCA or a SPCA ?" display="radio" help="If you choose PCA there is no rescaling of the data as oppposed as if you choose SPCA">
-            <option value="SPCA">SPCA</option>
-            <option value="PCA">PCA</option>
-        </param>
-    </inputs>
-    <outputs>
-        <data name="output_div" from_work_dir="Diversity.tabular" format="tabular" label="Global diversity">
-        </data>
-        <data name="output_beta" from_work_dir="BrayCurtis.tabular" format="tabular" label="Bray Curtis">
-            <filter>choice == 'Y'</filter>
-        </data>
-        <collection type="list" name="plots" label="Comparison plot">
-            <discover_datasets pattern="(?P&lt;designation&gt;.+)\.png" visible="false" format="png"/>
-            <filter>choice =='Y'</filter>
-        </collection>
-    </outputs>
-    <tests>
-        <test>
-            <param name="type" value="envi_bil"/>
-            <param name="input_raster" value="S2A_Subset"/>
-            <param name="input_header" value="S2A_Subset.hdr"/>
-            <param name="input_plot" value="S2A_T33NUD_Plots.zip"/>
-            <param name="choice" value="Y"/>
-            <output name="output_div">
-                <assert_contents>
-                    <has_n_lines n="25"/>
-                </assert_contents>
-            </output>
-            <output name="output_beta">
-                <assert_contents>
-                    <has_n_lines n="25"/>
-                </assert_contents>
-            </output>
-            <output_collection name="plots" type="list" count="1"/>
-        </test>
-    </tests>
-    <help><![CDATA[
-========================================================================
-Process satellite remote sensing data to produce biodiversity indicators
-========================================================================
-
-
-**What it does**
-
-Féret and Asner (2014) developed a method for **tropical forest** diversity mapping based on very high spatial resolution airborne imaging spectroscopy.
-
-The goal of this tool using the package biodivMapR is to compute diversity indices over each spatial polygon of a shapefile of plots, if available, in order to compare field inventories with diversity indices estimated from remotely-sensed images.
-
-**Input description**
-
-It expects an image file as input, with a specific data format. ENVI HDR image with BIL interleave required.
-The image is an ENVI raster including :
-
-- A binary file (which has no extension here).
-
-- A header file (with .hdr extension).
-
-The header file is a text file including all necessary metadata which can be read with a text editor. It includes image dimensions, projection, and the name and central wavelength for each spectral band.
-
-In order to get such input we advise to use the tool preprocessing sentinel 2 data. 
-
-+--------------+----------+--------------+
-|      BIL     | ENVI HDR |  Shapefiles  |
-+==============+==========+==============+
-| raster stack | Metadata |  plots.zip   |
-+--------------+----------+--------------+
-|      ...     |    ...   |      ...     |
-+--------------+----------+--------------+
-
-**Output**
-
-- Two tabulars : 
-    - One matrix for Bray-Curtis indicator
-    - One table for the following indicators; Species richness, shannon, fisher, simpson, richness, eveness, divergence
-
-- One comparison png plot in the Pcoa space that summarizes α- and β-diversity in scatterplots and illustrates that the combination of the three components computed with PCoA allows proper differentiation among vegetation types:
-    - PCoA#1 allows differentiating medium and high diversity forests from low diversity forest and low vegetation, but does not discriminate medium and high diversity forests.
-    - PCoA#2 allows differentiating low diversity forest from medium/high diversity forests and low vegetation
-    - PCoA#3 allows differentiating medium diversity forests from high diversity forests and low vegetation.
-
-    ]]>    </help>
-        <expand macro="SRS_BDMRref"/>
-</tool>

--- a/functions.r	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,79 +0,0 @@
-#Rscript
-
-#########################
-##      Functions      ##
-#########################
-
-#####Packages : raster
-#               sp
-#               ggplot2
-
-####Set paramaters for all tools using BiodivMapR
-
-# path for the Mask raster corresponding to image to process
-# expected to be in ENVI HDR format, 1 band, integer 8bits
-# expected values in the raster: 0 = masked, 1 = selected
-# set to FALSE if no mask available
-input_mask_file <- FALSE
-
-# relative or absolute path for the Directory where results will be stored
-# For each image processed, a subdirectory will be created after its name
-output_dir <- "RESULTS"
-
-# SPATIAL RESOLUTION
-# resolution of spatial units for alpha and beta diversity maps (in pixels), relative to original image
-# if Res.Map = 10 for images with 10 m spatial resolution, then spatial units will be 10 pixels x 10m = 100m x 100m surfaces
-# rule of thumb: spatial units between 0.25 and 4 ha usually match with ground data
-# too small window_size results in low number of pixels per spatial unit, hence limited range of variation of diversity in the image
-window_size <- 10
-
-# PCA FILTERING: Set to TRUE if you want second filtering based on PCA outliers to be processed. Slower
-filterpca <- TRUE
-
-################################################################################
-##                    DEFINE PARAMETERS FOR METHOD                            ##
-################################################################################
-nbcpu <- 4
-maxram <- 0.5
-nbclusters <- 50
-
-################################################################################
-##                              PROCESS IMAGE                                 ##
-################################################################################
-# 1- Filter data in order to discard non vegetated / shaded / cloudy pixels
-ndvi_thresh <- 0.5
-blue_thresh <- 500
-nir_thresh  <- 1500
-continuum_removal <- TRUE
-
-
-
-#### Convert raster to dataframe
-
-# Convert raster to SpatialPointsDataFrame
-convert_raster <- function(data_raster) {
-r_pts <- raster::rasterToPoints(data_raster, spatial = TRUE)
-
-# reproject sp object
-geo_prj <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
-r_pts <- sp::spTransform(r_pts, sp::CRS(geo_prj))
-
-
-# Assign coordinates to @data slot, display first 6 rows of data.frame
-r_pts@data <- data.frame(r_pts@data, longitude = sp::coordinates(r_pts)[, 1],
-                         latitude = sp::coordinates(r_pts)[, 2])
-
-return(r_pts@data)
-}
-
-
-#### Potting indices
-
-plot_indices <- function(data, titre) {
-  graph_indices <- ggplot2::ggplot(data) +
-  ggplot2::geom_point(ggplot2::aes_string(x = data[, 2], y = data[, 3], color = data[, titre]), shape = "square", size = 2) + ggplot2::scale_colour_gradient(low = "blue", high = "orange", na.value = "grey50") +
-  ggplot2::xlab("Longitude") + ggplot2::ylab("Latitude") +
-  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1), plot.title = ggplot2::element_text(color = "black", size = 12, face = "bold")) + ggplot2::ggtitle(titre)
-
-ggplot2::ggsave(paste0(titre, ".png"), graph_indices, width = 12, height = 10, units = "cm")
-}

--- a/indices_spectral.r	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,105 +0,0 @@
-#Rscript
-
-###########################################
-##    Mapping alpha and beta diversity   ##
-###########################################
-
-#####Packages :  expint,
-#                pracma,
-#                R.utils,
-#                raster,
-#                sp,
-#                matrixStats,
-#                ggplot2,
-#                expandFunctions,
-#                stringr,
-#                XML,
-#                rgdal,
-#                stars,
-#####Load arguments
-
-args <- commandArgs(trailingOnly = TRUE)
-
-#####Import the S2 data
-
-if (length(args) < 1) {
-    stop("This tool needs at least 1 argument")
-}else {
-    data_raster <- args[1]
-    data_header <- args[2]
-    data <- args[3]
-    source(args[4])
-    source(args[5])
-    source(args[6])
-    indice_choice <- strsplit(args[7], ",")[[1]]
-    source(args[8])
-    output_raster <- as.character(args[9])
-
-}
-
-########################################################################
-##                  COMPUTE SPECTRAL INDEX : NDVI                     ##
-########################################################################
-
-if (data != "") {
-  #Create a directory where to unzip your folder of data
-  dir.create("data_dir")
-  unzip(data, exdir = "data_dir")
-
-  # Read raster
-  data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl")
-  data_raster <- file.path("data_dir/results/Reflectance", data_raster[1])
-  refl <- raster::brick(data_raster)
-  refl2 <- raster::raster(data_raster)
-} else {
-  # Read raster
-  refl <- raster::brick(data_raster)
-  refl2 <- raster::raster(data_raster)
-}
-# get raster band name and clean format. Expecting band name and wav
-# reflFactor = 10000 when reflectance is coded as INT16
-refl <- raster::aggregate(refl, fact = 10)
-
-# Convert raster to SpatialPointsDataFrame
-refl2 <- raster::aggregate(refl2, fact = 10)
-r_pts <- convert_raster(refl2)
-table_ind <- r_pts
-# create directory for Spectralelength to be documented in image
-hdr_refl <- read_envi_header(get_hdr_name(data_raster))
-sensorbands <- hdr_refl$wavelength
-# compute a set of spectral indices defined by indexlist from S2 data indices
-si_path <- file.path("SpectralIndices")
-dir.create(path = si_path, showWarnings = FALSE, recursive = TRUE)
-# Save spectral indices
-
-indices <- lapply(indice_choice, function(x) {
-  indices_list <- computespectralindices_raster(refl = refl, sensorbands = sensorbands,
-                                                  sel_indices = x,
-                                                  reflfactor = 10000, stackout = FALSE)
-
-  index_path <- file.path(si_path, paste(basename(data_raster), "_", x, sep = ""))
-  spec_indices <- stars::write_stars(stars::st_as_stars(indices_list$spectralindices[[1]]), dsn = index_path, driver = "ENVI", type = "Float32")
-
-  # write band name in HDR
-  hdr <- read_envi_header(get_hdr_name(index_path))
-  hdr$`band names` <- x
-  write_envi_header(hdr = hdr, hdrpath = get_hdr_name(index_path))
-  # Writting the tabular and the plot
-  r_pts[, x] <- as.data.frame(sapply(spec_indices, c))
-  plot_indices(data = r_pts, titre = x)
-  return(r_pts)
-})
-
-new_table <- as.data.frame(indices)
-new_table <- new_table[, !grepl("longitude", names(new_table))]
-new_table <- new_table[, !grepl("latitude", names(new_table))]
-new_table <- new_table[, !grepl(basename(data_raster), names(new_table))]
-
-table_ind <- cbind(table_ind, new_table)
-if (length(indice_choice) == 1) {
-  colnames(table_ind) <- c(basename(data_raster), "longitude", "latitude", indice_choice)
-}
-
-write.table(table_ind, file = "Spec_Index.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
-
-# Get the raster layer of the indice as an output

--- a/indices_spectral.xml	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,270 +0,0 @@
-<tool id="srs_spectral_indices" name="Compute spectral indices" version="@VERSION@" profile = "20.01">
-    <description>as NDVI from remote sensing data</description>
-    <macros>
-        <import>macro.xml</import>
-    </macros>
-    <requirements>
-        <requirement type="package" version="4.2.2">r-base</requirement>
-        <requirement type="package" version="2.12.2">r-r.utils</requirement>
-        <requirement type="package" version="3.5_21">r-raster</requirement>
-        <requirement type="package" version="1.5_1">r-sp</requirement>
-        <requirement type="package" version="1.5_32">r-rgdal</requirement>
-        <requirement type="package" version="0.3.1">r-rasterdiv</requirement>
-        <requirement type="package" version="0.5_5">r-stars</requirement>
-        <requirement type="package" version="1.5.0">r-stringr</requirement>
-        <requirement type="package" version="2.4.2">r-pracma</requirement>
-        <requirement type="package" version="0.1_8">r-expint</requirement>
-        <requirement type="package" version="0.63.0">r-matrixstats</requirement>
-        <requirement type="package" version="3.4.0">r-ggplot2</requirement>
-        <requirement type="package" version="3.99_0.13">r-xml</requirement>
-        <requirement type="package" version="2.2.2">r-zip</requirement>
-        <requirement type="package" version="1.0.10">r-dplyr</requirement>
-    </requirements>
-    <command detect_errors="exit_code"><![CDATA[ 
-         #import re 
-        #if $method.type == 'envi_bil': 
-          #set input_raster = $method.input_raster
-          #set input_raster_identifier = re.sub('[^\s\w\-]', '_', str($input_raster.element_identifier)) 
-          #set input_header = $method.input_header
-          #set input_header_identifier = re.sub('[^\s\w\-]+[^.hdr]', '_', str($input_header.element_identifier)) 
-          cp '${input_raster}' '${input_raster_identifier}' &&
-          cp '${input_header}' '${input_header_identifier}' &&
-        #end if 
-        Rscript
-            '$__tool_directory__/indices_spectral.r'
-             #if $method.type == 'envi_bil':
-              '$input_raster_identifier' 
-              '$input_header_identifier'
-              ''
-            #else:
-              ''
-              ''
-              '$method.input'
-            #end if
-            '$__tool_directory__/prosail-master/R/Lib_PROSAIL.R'
-            '$__tool_directory__/prosail-master/R/Lib_SpectralIndices.R'
-            '$__tool_directory__/prosail-master/R/Lib_PROSAIL_HybridInversion.R'
-            '$input_indice'
-            '$__tool_directory__/functions.r'
-            '$choice'
-            '$output_indices'
-            '$plots'
-            '$output_raster'
-        ]]>
-    </command>
-    <inputs>
-        <conditional name="method">
-            <param name="type" type="select" label="In which format are your data ?" help="You can directly take the Reflectance output from the tool 'Preprocessing sentinel 2 data'">
-                <option value="zipper">The data you are using are in a zip folder Reflectance</option>
-                <option value="envi_bil">You already have the files ENVI BIL</option>
-            </param>
-            <when value="zipper">
-                <param name="input" type="data" format="zip" multiple="true" label="Input data"/>
-            </when>
-            <when value="envi_bil">
-                <param name="input_raster" type="data" format="bil" label="Input raster"/>
-                <param name="input_header" type="data" format="hdr" label="Input raster header"/>
-            </when>
-        </conditional>
-        <param name="input_indice" type="select" label="Input the type of indice you want" multiple="true">
-            <option value="ARI1">ARI1</option>
-            <option value="ARI2">ARI2</option>
-            <option value="ARVI">ARVI</option>
-            <option value="BAI">BAI</option>
-            <option value="BAIS2">BAIS2</option>
-            <option value="CCCI">CCCI</option>
-            <option value="CHL_RE">CHL RE</option>
-            <option value="CRI1">CRI1</option>
-            <option value="CRI2">CRI2</option>
-            <option value="EVI">EVI</option>
-            <option value="EVI2">EVI2</option>
-            <option value="GRVI1">GRVI1</option>
-            <option value="GNDVI">GNDVI</option>
-            <option value="IRECI">IRECI</option>
-            <option value="LAI_SAVI">LAI SAVI</option>
-            <option value="MCARI">MCARI</option>
-            <option value="mNDVI705">mNDVI705</option>
-            <option value="MSAVI2">MSAVI2</option>
-            <option value="MSI">MSI</option>
-            <option value="mSR705">mSR705</option>
-            <option value="MTCI">MTCI</option>
-            <option value="nBR_RAW">nBR_RAW</option>
-            <option value="NDI_45">NDI_45</option>
-            <option value="NDII">NDII</option>
-            <option value="NDSI">NDSI</option>
-            <option value="NDVI">NDVI</option>
-            <option value="NDVI_G">NDVI_G</option>
-            <option value="NDVI705">NDVI705</option>
-            <option value="NDWI">NDWI</option>
-            <option value="NDWI1">NDWI1</option>
-            <option value="NDWI2">NDWI2</option>
-            <option value="PSRI">PSRI</option>
-            <option value="PSRI_NIR">PSRI_NIR</option>
-            <option value="RE_NDVI">RE_NDVI</option>
-            <option value="RE_NDWI">RE_NDWI</option>
-            <option value="S2REP">S2REP</option>
-            <option value="SAVI">SAVI</option>
-            <option value="SIPI">SIPI</option>
-            <option value="SR">SR</option>
-            <option value="CR_SWIR">CR_SWIR</option>
-        </param>
-        <param name="choice" type="select" label="Do you want the raster layer of the indice as an output ?" display="radio">
-            <option value="N">No</option>
-            <option value="Y">Yes</option>
-        </param>
-    </inputs>
-    <outputs>
-        <data name="output_indices" from_work_dir="Spec_Index.tabular" format="tabular" label="${input_indice} tabular"/>
-        <collection type="list" name="plots" label="${input_indice} plot">
-            <discover_datasets pattern="(?P&lt;designation&gt;.+)\.png" visible="false" format="png"/>
-        </collection>
-        <collection type="list" name="output_raster" label="${input_indice} raster">
-            <discover_datasets pattern="(?P&lt;designation&gt;.+)" directory="SpectralIndices/" visible="false" recurse="true" format="bil,hdr"/>
-                <filter> choice == 'Y'</filter>
-        </collection>
-    </outputs>
-    <tests>
-        <test>
-            <param name="type" value="envi_bil"/>
-            <param name="input_raster" value="S2A_Subset"/>
-            <param name="input_header" value="S2A_Subset.hdr"/>
-            <param name="input_indice" value="NDVI"/>
-            <param name="choice" value="N"/>
-            <output name="output_indices">
-                <assert_contents>
-                    <has_n_columns n="4"/>
-                    <has_n_lines n="10001"/>
-                    <has_text_matching expression="237\.25\t13.72\d*\t3.17\d*\t0.80\d*" />
-                </assert_contents>
-            </output>
-            <output_collection name="plots" type="list" count="1"/>
-        </test>
-    </tests>
-    <help><![CDATA[
-=========================================================================
-Computes biodiversity spectral indices from satellite remote sensing data 
-=========================================================================
-
-
-**What it does**
-
-This tool estimates vegetation properties (leaf and canopy) from sensor measurements. You can choose which vegetation property you want to compute.
-
-**Input description**
-
-It expects an image file as input, with a specific data format. ENVI HDR image with BIL interleave required.
-The image is an ENVI raster including :
-
-- A binary file (which has no extension here).
-
-- A header file (with .hdr extension).
-
-The header file is a text file including all necessary metadata which can be read with a text editor. It includes image dimensions, projection, and the name and central wavelength for each spectral band.
-
-In order to get such input we advise to use the tool preprocessing sentinel 2 data. If you did so you can directly enter the "Reflectance" output from this tool and thus select the otpion "The data you are using are in a zip folder Reflectance". 
-
-⚠️ If you do not use this Reflectance folder make sure that your data are respectively in bil and hdr format in the datatypes.
-
-Finally, you can choose whether or not you want to have the raster layer of your indice.
-
-+--------------+----------+---------------+
-|      BIL     |    HDR   |Spectral indice|
-+==============+==========+===============+
-| raster stack | Metadata |     NDVI      |
-+--------------+----------+---------------+
-|      ...     |    ...   |      ...      |
-+--------------+----------+---------------+
-
-**Output**
-
-- One tabular with 3 columns : longitude, latitude and the chosen indice. 
-
-- One png plot for the vizualisation of the chosen indice.
-
-- Opionnnal, the raster layer in a zip file.
-
-**Additionnal informations**
-
-This tools allows you to compute one of GEO BON EBV 'Canopy Chlorophyll Content' (https://portal.geobon.org/ebv-detail?id=13). This EBV is computed by GEO BON on the Netherlands, here you can compute it on which ever Sentinel 2 data you want by chosing to calculate the indice CCI. When you chose your data you have to select a cloud cover smaller than 10% (you can selct this character directly on the different platforms wher you can download Sentinel 2 data).
-
-**Indices description**
-
-+-------------------+---------------------------------------------------------+
-|    Abbreviation   |                    Name                                 |
-+===================+=========================================================+ 
-|  ARI              |Anthocyanin reflectance index                            |
-+-------------------+---------------------------------------------------------+
-|  ARVI             |Atmospherically Resistant Vegetation Index               |
-+-------------------+---------------------------------------------------------+
-|  BAI              |Burn Area Index                                          |
-+-------------------+---------------------------------------------------------+
-|  BAIS2            |Burned Area Index for Sentinel 2                         |
-+-------------------+---------------------------------------------------------+
-|  CCCI             |Canopy Chlorophyll Content Index                         |
-+-------------------+---------------------------------------------------------+
-|  CHL_RE           |Chlorophyll Red-Edge ?                                   |
-+-------------------+---------------------------------------------------------+
-|  EVI              |Enhanced Vegetation Index                                |
-+-------------------+---------------------------------------------------------+
-|  GRVI1            |Green Ratio Vegetation Index                             |
-+-------------------+---------------------------------------------------------+
-|  GNDVI            |Green Normalized Difference Vegetation Index             |
-+-------------------+---------------------------------------------------------+
-|  IRECI            |Inverted Red-Edge Chlorophyll Index                      |
-+-------------------+---------------------------------------------------------+
-|  LAI SAVI         |Leaf Area Index Soil Adjusted Vegetation Index ?         |
-+-------------------+---------------------------------------------------------+
-|  MCARI            |Modified Chlorophyll Absorption in Reflectance Index     |
-+-------------------+---------------------------------------------------------+
-|  mNDVI705         |Modified NDVI 705                                        |
-+-------------------+---------------------------------------------------------+
-|  MSAVI2           |Modified Soil Adjusted Vegetation Index 2                |
-+-------------------+---------------------------------------------------------+
-|  MSI              |Moisture Stress Index                                    |
-+-------------------+---------------------------------------------------------+
-|  mSR705           |Modified Simple Ratio 705                                |
-+-------------------+---------------------------------------------------------+
-|  MTCI             |MERIS Terrestrial Chlorophyll Index                      |
-+-------------------+---------------------------------------------------------+
-|  NBR_RAW          |Normalized Burn Ratio RAW                                |
-+-------------------+---------------------------------------------------------+
-|  NDI_45           |                                                         |
-+-------------------+---------------------------------------------------------+
-|  NDII             |Normalized Difference 819/1600 NDII                      |
-+-------------------+---------------------------------------------------------+
-|  NDSI             |Normalized Difference Snow Index                         |
-+-------------------+---------------------------------------------------------+
-|  NDVI             |Normalized Difference Vegetation Index                   |
-+-------------------+---------------------------------------------------------+
-|  NDVI_G           |                                                         |
-+-------------------+---------------------------------------------------------+
-|  NDVI705          |NDVI 705                                                 |
-+-------------------+---------------------------------------------------------+
-|  NDWI             |Normalized Difference Water Index                        |
-+-------------------+---------------------------------------------------------+
-|  NDWI1            |NDMI ? Normalized Difference Moisture Index ?            |
-+-------------------+---------------------------------------------------------+
-|  NDWI2            |                                                         |
-+-------------------+---------------------------------------------------------+
-|  PSRI             |Plant Senescing Reflectance Index                        |
-+-------------------+---------------------------------------------------------+
-|  PSRI_NIR         |                                                         |
-+-------------------+---------------------------------------------------------+
-|  RE_NDVI          |Red Edge NDVI ?                                          |
-+-------------------+---------------------------------------------------------+
-|  RE_NDWI          |Red Edge NDWI                                            |
-+-------------------+---------------------------------------------------------+
-|  S2REP            |Sentinel 2 Red-Edge Position                             |
-+-------------------+---------------------------------------------------------+
-|  SAVI             |Soil Adjusted Vegetation Index                           |
-+-------------------+---------------------------------------------------------+
-|  SIPI             |Structure Intensive Pigment Index 1                      |
-+-------------------+---------------------------------------------------------+
-|  SR               |Simple Ratio ?                                           |
-+-------------------+---------------------------------------------------------+
-|  CR_SWIR          |CR short wave infrared                                   |
-+-------------------+---------------------------------------------------------+
-
-    ]]>    </help>
-        <expand macro="SRS_prosailref"/>
-</tool>

--- a/macro.xml	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,102 +0,0 @@
-<macros>
-    <token name="@VERSION@">0.0.1</token>
-    <xml name="SRS_requirements">
-        <requirements>
-            <requirement type="package" version="4.2.2">r-base</requirement>
-            <requirement type="package" version="2.12.2">r-r.utils</requirement>
-            <requirement type="package" version="1.0_7">r-sf</requirement>
-            <requirement type="package" version="1.5_32">r-rgdal</requirement>
-            <requirement type="package" version="0.6_1">r-rgeos</requirement>
-            <requirement type="package" version="0.5_5">r-stars</requirement>
-            <requirement type="package" version="1.5.0">r-stringr</requirement>
-            <requirement type="package" version="1.9.8">r-biodivmapr</requirement>
-            <requirement type="package" version="0.36">r-xfun</requirement>
-            <requirement type="package" version="0.51.4">r-rastervis</requirement>
-            <requirement type="package" version="3.4.0">r-ggplot2</requirement>
-            <requirement type="package" version="2.11.0">r-mapview</requirement>
-            <requirement type="package" version="2.3">r-gridextra</requirement>
-            <requirement type="package" version="3.0.0">r-emstreer</requirement>
-            <requirement type="package" version="3.2.3">r-filesstrings</requirement>
-            <!--requirement type="package" version="9.0.1">proj</requirement-->
-            <yield/>
-        </requirements>
-    </xml>
-    <xml name="SRS_input">
-       <param name="input" type="data" format="tabular,csv,txt" label="Input table"/>
-       <param name="colnames" type="boolean" label="First line is a header line" checked="true"/>
-    </xml>
-    <xml name="SRS_bibref">
-       <citations>
-            <citation type="bibtex">
-            @Manual{,
-            title = {obisindicators: Develop marine biodiversity indicators from OBIS data},
-            author = {Ben Ben and Pieter Provoost and Tylar Murray},
-            year = {2022},
-            note = {https://marinebon.org/obisindicators,
-            https://github.com/marinebon/obisindicators},
-            }
-            </citation>
-        </citations>
-    </xml>
-    <xml name="SRS_S2ref">
-        <citations>
-            <citation type="bibtex">
-            @Manual{,
-            title = {preprocS2: preprocessing of Sentinel-2 data downloaded from various data hubs / produced from various atmospheric correction methods},
-            author = {Jean-Baptiste Feret},
-            year = {2022},
-            note = {R package version 1.1.3},
-            url = {https://gitlab.com/jbferet/preprocS2},
-            }
-            </citation>   
-        </citations>
-    </xml>
-    <xml name="SRS_BDMRref">
-        <citations>
-            <citation type="bibtex">
-            @Manual{,
-            title = {biodivMapR: biodivMapR: an R package for a- and ß-diversity mapping using remotely-sensed images},
-            author = {Jean-Baptiste Feret and Florian {de Boissieu}},
-            year = {2022},
-            note = {R package version 1.9.4},
-            url = {https://github.com/jbferet/biodivMapR},
-            }
-            </citation>   
-        </citations>
-    </xml>
-    <xml name="SRS_rasterdivref">
-        <citations>
-            <citation type="doi">doi:10.1016/j.ecolind.2016.07.039</citation> 
-            <citation type="doi">doi:10.1101/2021.01.23.427872</citation>     
-        </citations>
-    </xml>
-    <xml name="SRS_prosailref">
-        <citations>
-            <citation type="bibtex">
-            @Manual{,
-            title = {prosail: PROSAIL leaf and canopy radiative transfer model and inversion routines},
-            author = {Jean-Baptiste Feret and Florian {de Boissieu}},
-            year = {2022},
-            note = {R package version 1.1.1},
-            url = {https://gitlab.com/jbferet/prosail},
-            }
-            </citation>   
-        </citations>
-    </xml>
-    <xml name="SRS_metaref">
-        <citations>
-            <citation type="doi">doi.org/10.5281/zenodo.5907920</citation>
-        </citations>
-    </xml>
-    <xml name="SRS_vegetpheno">
-        <citations>
-            <citation type="doi">doi.org/10.1016/j.rse.2013.09.022</citation>
-        </citations>
-    </xml>
-    <xml name="topic">
-        <edam_topics>
-           <edam_topic>topic_0610</edam_topic>
-           <edam_topic>topic_3050</edam_topic>
-        </edam_topics>
-    </xml>
-</macros>

--- a/pca_raster.r	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,75 +0,0 @@
-#Rscript
-
-###########################################
-##    Getting PCA raster   ##
-###########################################
-
-#####Packages : stars
-#               utils
-#               biodivmapr
-#               raster
-#               sf
-#               mapview
-#               leafpop
-#               RColorBrewer
-#               labdsv
-#               rgdal
-#               ggplot2
-#               gridExtra
-## remotes::install_github("jbferet/biodivMapR")
-#####Load arguments
-
-args <- commandArgs(trailingOnly = TRUE)
-
-#####Import the S2 data
-
-if (length(args) < 1) {
-    stop("This tool needs at least 1 argument")
-}else {
-    data_raster <- args[1]
-    rasterheader <- args[2]
-    data <- args[3]
-    typepca <- as.character(args[4])
-    source(args[5])
-}
-
-################################################################################
-##              DEFINE PARAMETERS FOR DATASET TO BE PROCESSED                 ##
-################################################################################
-# expected to be in ENVI HDR
-
-if (data_raster == "") {
-  #Create a directory where to unzip your folder of data
-  dir.create("data_dir")
-  unzip(data, exdir = "data_dir")
-  # Path to raster
-  data_raster <- list.files("data_dir/results/Reflectance", pattern = "_Refl")
-  input_image_file <- file.path("data_dir/results/Reflectance", data_raster[1])
-  input_header_file <- file.path("data_dir/results/Reflectance", data_raster[2])
-
-} else {
-  input_image_file <- file.path(getwd(), data_raster, fsep = "/")
-  input_header_file <- file.path(getwd(), rasterheader, fsep = "/")
-}
-
-################################################################################
-##                              PROCESS IMAGE                                 ##
-################################################################################
-# 1- Filter data in order to discard non vegetated / shaded / cloudy pixels
-print("PERFORM PCA ON RASTER")
-pca_output <- biodivMapR::perform_PCA(Input_Image_File = input_image_file, Input_Mask_File = input_mask_file,
-                          Output_Dir = output_dir, TypePCA = typepca, FilterPCA = filterpca, nbCPU = nbcpu, MaxRAM = maxram)
-
-
-pca_path <- file.path(output_dir, basename(data_raster), typepca, "PCA", "OutputPCA_8_PCs")
-pca_raster <- raster::raster(pca_path)
-get_pca <- convert_raster(pca_raster)
-
-colnames(get_pca) <- c("PCA", "longitude", "latitude")
-plot_indices(get_pca, titre = "PCA")
-
-write.table(get_pca, file = "PCA.tabular", sep = "\t", dec = ".", na = " ", row.names = FALSE, col.names = TRUE, quote = FALSE)
-#### Get the raster layer files
-pca_files <- file.path("RESULTS", basename(data_raster), typepca, "PCA")
-to_dir_short <- output_dir
-file.copy(pca_files, to_dir_short) #copy files from long to short paths

--- a/pca_raster.xml	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,120 +0,0 @@
-<tool id="srs_pca" name="Compute a PCA raster" version="@VERSION@" profile = "20.01">
-    <description>from remote sensing data</description>
-    <macros>
-        <import>macro.xml</import>
-    </macros>
-    <expand macro="SRS_requirements"/>
-    <command detect_errors="exit_code"><![CDATA[ 
-        #import re 
-        #if $method.origin == 'envi_bil': 
-          #set input_raster = $method.input_raster
-          #set input_raster_identifier = re.sub('[^\s\w\-]', '_', str($input_raster.element_identifier)) 
-          #set input_header = $method.input_header
-          #set input_header_identifier = re.sub('[^\s\w\-]+[^.hdr]', '_', str($input_header.element_identifier)) 
-          cp '${input_raster}' '${input_raster_identifier}' &&
-          cp '${input_header}' '${input_header_identifier}' &&
-        #end if
-        Rscript
-            '$__tool_directory__/pca_raster.r'
-             #if $method.origin == 'envi_bil':
-              '$input_raster_identifier' 
-              '$input_header_identifier'
-              ''
-            #else:
-              ''
-              ''
-              '$method.input'
-            #end if
-            '$type'
-            '$__tool_directory__/functions.r'
-            '$output_raster'
-            '$output_pca'
-            '$plots'
-        ]]>
-    </command>
-    <inputs>
-        <conditional name="method">
-            <param name="origin" type="select" label="In which format are your data ?">
-                <option value="zipper">The data you are using are in a zip folder Reflectance</option>
-                <option value="envi_bil">Your already have the files ENVI BIL</option>
-            </param>
-            <when value="zipper">
-                <param name="input" type="data" format="zip" multiple="true" label="Input data"/>
-            </when>
-            <when value="envi_bil">
-                <param name="input_raster" type="data" format="bil" label="Input raster"/>
-                <param name="input_header" type="data" format="hdr" label="Input raster header"/>
-            </when>
-        </conditional>
-        <param name="type" type="select" label="Do you want to do a PCA or a SPCA ?" display="radio" help="If you choose PCA there is no rescaling of the data as oppposed as if you choose SPCA">
-            <option value="SPCA">SPCA</option>
-            <option value="PCA">PCA</option>
-        </param>
-    </inputs>
-    <outputs>
-        <collection type="list" name="output_raster" label="PCA raster">
-            <discover_datasets pattern="(?P&lt;designation&gt;.+)" directory="RESULTS/" visible="false" recurse="true" format="bil,hdr"/>
-        </collection>
-        <data name="output_pca" from_work_dir="PCA.tabular" format="tabular" label="PCA tabular"/>
-        <collection type="list" name="plots" label="PCA plot">
-            <discover_datasets pattern="(?P&lt;designation&gt;.+)\.png" visible="false" format="png"/>
-        </collection>
-    </outputs>
-    <tests>
-        <test>
-            <param name="origin" value="envi_bil"/>
-            <param name="input_raster" value="S2A_Subset"/>
-            <param name="input_header" value="S2A_Subset.hdr"/>
-            <param name="type" value="PCA"/>
-            <output_collection name="output_raster" type="list" count="5"/>
-            <output name="output_pca">
-                <assert_contents>
-                    <has_n_lines n="960435"/>
-                </assert_contents>
-            </output>
-            <output_collection name="plots" type="list" count="1"/>
-        </test>
-    </tests>
-    <help><![CDATA[
-========================================================================
-Process satellite remote sensing data to produce biodiversity indicators
-========================================================================
-
-
-**What it does**
-
-Féret and Asner (2014) developed a method for **tropical forest** diversity mapping based on very high spatial resolution airborne imaging spectroscopy.
-
-The goal of this tool using the package biodivMapR is to produce a Principal Components Analysis (PCA) based on (optical) images. Principal component analysis prepares for feature extraction in order to reduce noise in the image. It is performed on a random subset of the image to ensure computational efficiency for large images.
-
-**Input description**
-
-It expects an image file as input, with a specific data format. ENVI HDR image with BIL interleave required.
-The image is an ENVI raster including :
-
-- A binary file (which has no extension here).
-
-- A header file (with .hdr extension).
-
-The header file is a text file including all necessary metadata which can be read with a text editor. It includes image dimensions, projection, and the name and central wavelength for each spectral band.
-
-In order to get such input we advise to use the tool preprocessing sentinel 2 data. 
-
-+--------------+----------+
-|      BIL     | ENVI HDR |
-+==============+==========+
-| raster stack | Metadata |
-+--------------+----------+
-|      ...     |    ...   |
-+--------------+----------+
-
-**Output**
-
-- One tabular with 3 columns : longitude, latitude and pca.
-
-- One png graph.
-
-- One  ENVI BIL files with the one file for the bil raster and one file for the hdr header.
-    ]]>    </help>
-        <expand macro="SRS_BDMRref"/>
-</tool>

--- a/preprocess_S2.r	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,121 +0,0 @@
-#Rscript
-
-###########################################
-##     Preprocessing Sentinel 2 data     ##
-###########################################
-
-#####Packages : sen2r,
-#    jqr,
-#    protolite,
-#    raster,
-#    sf,
-#    rgeos,
-#    sp,
-#    raster,
-#    stars,
-#    stringr,
-#    progress,
-#    rgdal,
-#    R.utils,
-#    gdalUtils,
-#    fasterize,
-#    XML,
-#    XML2
-
-#####Load arguments
-
-args <- commandArgs(trailingOnly = TRUE)
-
-if (length(args) < 1) {
-    stop("This tool needs at least 1 argument")
-}else {
-    data <- args[1]
-    source(args[2])
-    data_source <- as.character(args[3])
-    sat_type <- as.character(args[4])
-}
-
-##____________________________________________________________________##
-##        Define where data is stored and where to write results      ##
-##--------------------------------------------------------------------##
-
-#Create a directory where to unzip your folder of data
-dir.create("data_dir")
-unzip(data, exdir = "data_dir")
-
-# Result directory
-result_path <- "results"
-dir.create(path = result_path, showWarnings = FALSE, recursive = TRUE)
-
-#Csv file for output useless but needed for linter
-write.csv(data_source, "Mission.csv")
-
-# define raster path
-if (data_source == "SAFE") {
-    path_s2 <- file.path("data_dir", list.files("data_dir", pattern = ".SAFE"))
-    #To define the level and know if a correction is needed (convert not ready yet)
-    level_info <- get_s2_level(path_s2)
-    if (level_info == "L1C") {
-        stop("! This tool works for data of L2A level and NOT for the L1C level which is currently a work in progress !")
-    }
-}else {
-    path_s2 <- file.path("data_dir")
-}
-
-##____________________________________________________________________##
-##                  Extract, resample & stack data                    ##
-##--------------------------------------------------------------------##
-# define resolution
-resolution <- 10
-# define source of data
-s2source <- data_source
-
-s2obj <- extract_from_s2_l2a(path_dir_s2 = path_s2,
-                                        path_vector = NULL,
-                                        s2source = s2source,
-                                        resolution = resolution)
-
-##____________________________________________________________________##
-##                        Write CLOUD MASK                            ##
-##--------------------------------------------------------------------##
-
-# directory for cloud mask
-cloud_path <- file.path(result_path, "CloudMask")
-dir.create(path = cloud_path, showWarnings = FALSE, recursive = TRUE)
-# Filename for cloud mask
-cloudmasks <- save_cloud_s2(s2_stars = s2obj$s2_stack,
-                                       cloud_path = cloud_path,
-                                       s2source = s2source, saveraw = TRUE)
-
-zip_cloud <- file.path("Cloud.zip")
-zip::zip(zip_cloud, cloud_path)
-##____________________________________________________________________##
-##                        Write REFLECTANCE                           ##
-##--------------------------------------------------------------------##
-
-# directory for Reflectance
-refl_dir <- file.path(result_path, "Reflectance")
-dir.create(path = refl_dir, showWarnings = FALSE, recursive = TRUE)
-
-if (data_source == "SAFE") {
-    # filename for Reflectance
-    refl_path <- file.path(refl_dir, paste(basename(s2obj$s2_bands$GRANULE), "_Refl", sep = ""))
-
-    # Save Reflectance file as ENVI image with BIL interleaves
-    tile_s2 <- substring(strsplit(basename(s2obj$s2_bands$GRANULE), "_")[[1]][2], 2)
-    dateacq_s2 <- as.Date(substring(strsplit(basename(s2obj$s2_bands$GRANULE), "_")[[1]][4], 1, 8), format = "%Y%m%d")
-}else {
-    # filename for Reflectance
-    refl_path <- file.path(refl_dir, paste(basename(s2obj$s2_bands$path_tile_s2), "_Refl", sep = ""))
-
-    # Save Reflectance file as ENVI image with BIL interleaves
-    tile_s2 <- substring(strsplit(basename(s2obj$s2_bands$path_tile_s2), "_")[[1]][2], 2)
-    dateacq_s2 <- as.Date(substring(strsplit(basename(s2obj$s2_bands$path_tile_s2), "_")[[1]][4], 1, 8), format = "%Y%m%d")
-}
-
-save_data <- save_reflectance_s2(s2_stars = s2obj$s2_stack, refl_path = refl_path,
-                               s2sat = sat_type, tile_s2 = tile_s2, dateacq_s2 = dateacq_s2,
-                               format = "ENVI", datatype = "Int16", mtd = s2obj$s2_bands$metadata, mtd_msi = s2obj$s2_bands$metadata_MSI)
-
-zip_files <- file.path("Refl.zip")
-zip::zip(zip_files, refl_dir)

--- a/preprocess_S2.xml	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,109 +0,0 @@
-<tool id="srs_preprocess_s2" name="Preprocess sentinel 2 data" version="@VERSION@" profile="20.01">
-    <description>read, crop, resample and write it as a raster stack</description>
-    <macros>
-        <import>macro.xml</import>
-    </macros>
-    <requirements>
-        <requirement type="package" version="4.2.2">r-base</requirement>
-        <requirement type="package" version="1.0_7">r-sf</requirement>
-        <requirement type="package" version="0.6_1">r-rgeos</requirement>
-        <requirement type="package" version="0.5_5">r-stars</requirement>
-        <requirement type="package" version="1.5.0">r-stringr</requirement>
-        <requirement type="package" version="1.2.3">r-jqr</requirement>
-        <requirement type="package" version="2.1.3">r-protolite</requirement>
-        <requirement type="package" version="1.5.3">r-sen2r</requirement>
-        <requirement type="package" version="1.2.2">r-progress</requirement>
-        <requirement type="package" version="2.0.3">r-gdalutils</requirement>
-        <requirement type="package" version="1.0.4">r-fasterize</requirement> 
-        <requirement type="package" version="3.99_0.13">r-xml</requirement> 
-        <requirement type="package" version="1.3.3">r-xml2</requirement>
-        <requirement type="package" version="2.2.2">r-zip</requirement>
-    </requirements>
-    <command detect_errors="exit_code"><![CDATA[
-        Rscript
-            '$__tool_directory__/preprocess_S2.r'
-            '$input'
-            '$__tool_directory__/Lib_preprocess_S2.r'
-            '$source_s2'
-            '$type_s2'
-            '$output'
-            '$output_refl'
-            '$output_cloud'
-        ]]>
-    </command>
-    <inputs>
-        <param name="input" type="data" format="zip" multiple="true" label="Input data"/>
-        <param name="source_s2" type="select" label="Where does your data come from ?">
-            <option value="SAFE">From Scihub or Peps</option>
-            <option value="THEIA">From Theia</option>
-        </param>     
-        <param name="type_s2" type="select" label="From which satellite are the data ?">
-            <option value="2A">Sentinel 2A</option>
-            <option value="2B">Sentinel 2B</option>
-        </param> 
-    </inputs>
-    <outputs>
-        <data name="output" from_work_dir="Mission.csv" format="csv" label="Mission source"/>
-        <data name="output_refl" from_work_dir="Refl.zip" format="zip" label="Reflectance"/>
-        <data name="output_cloud" from_work_dir="Cloud.zip" format="zip" label="Cloud Mask"/>
-    </outputs>
-    <tests>
-        <test>
-            <param name="input" value="S2A_MSIL2A_20200306T015621_N0214_R117_T51JXN_20200306T034744.zip"/>
-            <param name="source_s2" value="SAFE"/>
-            <param name="type_s2" value="2A"/>
-            <output name="output" value="Mission.csv"/>
-            <output name="output_refl">
-                <assert_contents>
-                    <has_size value="12923807" delta="100"/>
-                </assert_contents>
-            </output>
-            <output name="output_cloud">
-                <assert_contents>
-                    <has_size value="235912" delta="100"/>
-                </assert_contents>
-            </output>
-        </test>
-    </tests>
-    <help><![CDATA[
-==========================
-Preprocess Sentinel 2 data
-==========================
-
-
-**What it does**
-
-The goal of this tool (using the package preprocS2) is to provide a common framework for the preprocessing of Level-2A Sentinel-2 images (does not work yet for Level 1C).
-Sentinel-2 L2A images can be produced or obtained from various data hubs or atmospheric correction methods. PreprocS2 provides a unique function to read, crop, resample the original image directory, and write it as a raster stack.
-
-**Input description**
-
-A zip folder with Sentinel 2 data. 
-These data can be dowloaded on 3 different platforms: 
-
-- Copernicus Open Access Hub, Scihub, which provides complete, free and open access to Sentinel-2 data : https://scihub.copernicus.eu/dhus/#/home
-
-- PEPS, la 'Plateforme d'Exploitation de Produits Sentinel' : https://peps.cnes.fr/rocket/#/search
-
-- Theia : https://theia.cnes.fr/atdistrib/rocket/#/search
-
-For each of those 3 platforms you will need to create an account.
-
-You need to select from where you dowloaded your data in order for the tool to know the format of the folder once it unzips it.
-
-+----------+----------+
-|   ZIP    |  Source  |
-+==========+==========+
-|folder.zip|Character |
-+----------+----------+
-|   ...    |    ...   |
-+----------+----------+
-
-**Output**
-
-A zip folder containing the refelctance data with a raster stack pile in BIL format (no extension), a header with the metadata (.hdr extension) and 2 other files in .xml.
-
-A zip folder containing the cloud mask with the 2 files for the raw data (.RAW) and 2 files containing the metadata (.hdr).
-    ]]>    </help>
-        <expand macro="SRS_S2ref"/>
-</tool>

--- a/prosail-master/R/Lib_PROSAIL.R	Sun Jan 08 23:03:35 2023 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,1233 +0,0 @@
-# ============================================================================= =
-# prosail
-# Lib_PROSAIL.R
-# ============================================================================= =
-# PROGRAMMERS:
-# Jean-Baptiste FERET <jb.feret@teledetection.fr >
-# Florian de BOISSIEU <fdeboiss@gmail.com >
-# copyright 2019 / 11 Jean-Baptiste FERET
-# ============================================================================= =
-# This Library includes functions dedicated to PROSAIL simulation
-# SAIL versions available are 4SAIL and 4SAIL2
-# ============================================================================= =
-
-#" computes bidirectional reflectance factor based on outputs from PROSAIL and sun position
-#"
-#" The direct and diffuse light are taken into account as proposed by:
-#" Francois et al. (2002) conversion of 400-1100 nm vegetation albedo
-#" measurements into total shortwave broadband albedo using a canopy
-#" radiative transfer model,  Agronomie
-#" Es = direct
-#" Ed = diffuse
-#"
-#" @param rdot numeric. Hemispherical-directional reflectance factor in viewing direction
-#" @param rsot numeric. Bi-directional reflectance factor
-#" @param tts numeric. Solar zenith angle
-#" @param specatm_sensor list. direct and diffuse radiation for clear conditions
-#" @return brf numeric. Bidirectional reflectance factor
-#" @export
-compute_brf  <- function(rdot, rsot, tts, specatm_sensor) {
-
-  ############################## #
-  ##	direct  /  diffuse light	##
-  ############################## #
-  es <- specatm_sensor$Direct_Light
-  ed <- specatm_sensor$Diffuse_Light
-  rd <- pi / 180
-  skyl <- 0.847 - 1.61 * sin((90 - tts) * rd) +  1.04 * sin((90 - tts) * rd) * sin((90 - tts) * rd) # diffuse radiation (Francois et al.,  2002)
-  pardiro <- (1 - skyl) * es
-  pardifo <- skyl * ed
-  brf <- (rdot * pardifo + rsot * pardiro) / (pardiro + pardifo)
-  return(brf)
-}
-
-#" Performs PROSAIL simulation based on a set of combinations of input parameters
-#" @param spec_sensor list. Includes optical constants required for PROSPECT
-#" refractive index,  specific absorption coefficients and corresponding spectral bands
-#" @param input_prospect  list. PROSPECT input variables
-#" @param n numeric. Leaf structure parameter
-#" @param chl numeric. chlorophyll content (microg.cm-2)
-#" @param car numeric. carotenoid content (microg.cm-2)
-#" @param ant numeric. anthocyain content (microg.cm-2)
-#" @param brown numeric. brown pigment content (Arbitrary units)
-#" @param ewt numeric. Equivalent Water Thickness (g.cm-2)
-#" @param lma numeric. Leaf Mass per Area (g.cm-2)
-#" @param prot numeric. protein content  (g.cm-2)
-#" @param cbc numeric. nonprotcarbon-based constituent content (g.cm-2)
-#" @param alpha numeric. Solid angle for incident light at surface of leaf (simulation of roughness)
-#" @param typelidf numeric. Type of leaf inclination distribution function
-#" @param lidfa numeric.
-#" if typelidf  == 1,  controls the average leaf slope
-#" if typelidf  == 2,  corresponds to average leaf angle
-#" @param lidfb numeric.
-#" if typelidf  == 1,  unused
-#" if typelidf  == 2,  controls the distribution"s bimodality
-#" @param lai numeric. Leaf Area Index
-#" @param q numeric. Hot Spot parameter
-#" @param tts numeric. Sun zeith angle
-#" @param tto numeric. Observer zeith angle
-#" @param psi numeric. Azimuth Sun  /  Observer
-#" @param rsoil numeric. Soil reflectance
-#" @param fraction_brown numeric. Fraction of brown leaf area
-#" @param diss numeric. Layer dissociation factor
-#" @param cv numeric. vertical crown cover percentage
-#" = % ground area covered with crowns as seen from nadir direction
-#" @param zeta numeric. Tree shape factor
-#" = ratio of crown diameter to crown height
-#" @param sailversion character. choose between 4SAIL and 4SAIL2
-#" @param brownvegetation list. Defines optical properties for brown vegetation,  if not nULL
-#" - WVL,  reflectance,  Transmittance
-#" - Set to nULL if use PROSPECT to generate it
-#"
-#" @return list. rdot, rsot, rddt, rsdt
-#" rdot: hemispherical-directional reflectance factor in viewing direction
-#" rsot: bi-directional reflectance factor
-#" rsdt: directional-hemispherical reflectance factor for solar incident flux
-#" rddt: bi-hemispherical reflectance factor
-#" @import prospect
-#" @export
-pro4sail  <- function(spec_sensor, input_prospect = nULL, n = 1.5, chl = 40.0,
-                     car = 8.0, ant = 0.0, brown = 0.0, ewt = 0.01,
-                     lma = 0.008, prot = 0.0, cbc = 0.0, alpha = 40.0,
-                     typelidf = 2, lidfa = nULL, lidfb = nULL, lai = nULL,
-                     q = nULL, tts = nULL, tto = nULL, psi = nULL, rsoil = nULL,
-                     fraction_brown = 0.0,  diss = 0.0,  cv = 1, zeta = 1,
-                     sailversion = "4SAIL", brownvegetation = nULL) {
-
-  ############################ #
-  #	LEAF OPTICAL PROPERTIES	##
-  ############################ #
-  if (is.null(input_prospect)) {
-    input_prospect <- data.frame("chl" = chl, "car" = car, "ant" = ant, "brown" = brown, "ewt" = ewt,
-                                "lma" = lma, "prot" = prot, "cbc" = cbc, "n" = n, "alpha" = alpha)
-  }
-  greenvegetation <- prospect::PROSPECT(SpecPROSPECT = spec_sensor,
-                                        n = input_prospect$n[1],
-                                        chl = input_prospect$chl[1],
-                                        car = input_prospect$car[1],
-                                        ant = input_prospect$ant[1],
-                                        brown = input_prospect$brown[1],
-                                        ewt = input_prospect$ewt[1],
-                                        lma = input_prospect$lma[1],
-                                        prot = input_prospect$prot[1],
-                                        cbc = input_prospect$cbc[1],
-                                        alpha = input_prospect$alpha[1])
-
-  if (sailversion  ==  "4SAIL2") {
-    # 4SAIL2 requires one of the following combination of input parameters
-    # Case #1: valid optical properties for brown vegetation
-    if (!is.null(brownvegetation)) {
-      # need to define reflectance and Transmittance for brownvegetation
-      if (length(grep("reflectance", names(brownvegetation))) == 0 || length(grep("Transmittance", names(brownvegetation))) == 0) {
-        message("Please define brownvegetation as a list including reflectance and Transmittance")
-        stop()
-      }
-      # check if spectral domain for optical properties of brown vegetation match
-      # with spectral domain for optical properties of green vegetation
-      if (length(setdiff(spec_sensor$lambda, brownvegetation$wvl)) > 0) {
-        message("Please define same spectral domain for brownvegetation and SpecPROSPECT")
-        stop()
-      }
-      if (length(unique(lengths(input_prospect))) == 1) {
-        if (!unique(lengths(input_prospect)) == 1) {
-          message("brownvegetation defined along with multiple leaf chemical properties")
-          message("Only first set of leaf chemical properties will be used to simulate green vegetation")
-        }
-      }
-    # if no leaf optical properties brown vegetation defined
-    } else if (is.null(brownvegetation)) {
-      # if all PROSPECT input parameters have the same length
-      if (length(unique(lengths(input_prospect))) == 1) {
-        # if all PROSPECT input parameters are unique (no possibility to simulate 2 types of leaf optics)
-        if (unique(lengths(input_prospect)) == 1) {
-          # if fraction_brown set to 0,  then assign green vegetation optics to brown vegetation optics
-          if (fraction_brown == 0) {
-            brownvegetation <- greenvegetation
-          # else run 4SAIL
-          } else {
-            message("4SAIL2 needs two sets of optical properties for green and brown vegetation")
-            message("Currently one set is defined. will run 4SAIL instead of 4SAIL2")
-            sailversion <- "4SAIL"
-          }
-        # if all PROSPECT parameters have at least 2 elements
-        } else if (unique(lengths(input_prospect)) >= 2) {
-          # compute leaf optical properties
-          brownvegetation <- prospect::PROSPECT(SpecPROSPECT = spec_sensor,
-                                                n = input_prospect$n[2],
-                                                chl = input_prospect$chl[2],
-                                                car = input_prospect$car[2],
-                                                ant = input_prospect$ant[2],
-                                                brown = input_prospect$brown[2],
-                                                ewt = input_prospect$ewt[2],
-                                                lma = input_prospect$lma[2],
-                                                prot = input_prospect$prot[2],
-                                                cbc = input_prospect$cbc[2],
-                                                alpha = input_prospect$alpha[2])
-          if (unique(lengths(input_prospect)) > 2) {
-            message("4SAIL2 needs two sets of optical properties for green and brown vegetation")
-            message("Currently more than 2 sets are defined. will only use the first 2")
-          }
-        }
-      }
-    }
-  }
-  if (sailversion  ==  "4SAIL") {
-    if (length(unique(lengths(input_prospect))) == 1) {
-      if (unique(lengths(input_prospect)) > 1) {
-        message("4SAIL needs only one set of optical properties")
-        message("Currently more than one set of leaf chemical constituents is defined.")
-        message("Will run 4SAIL with the first set of leaf chemical constituents")
-      }
-    }
-  }
-
-  if (sailversion  ==  "4SAIL") {
-    # run 4SAIL
-    ref <- foursail(leafoptics = greenvegetation,
-                    typelidf,  lidfa,  lidfb,  lai,  q,  tts,  tto,  psi,  rsoil)
-  } else if (sailversion  ==  "4SAIL2") {
-    # run 4SAIL2
-    ref <- foursail2(leafgreen = greenvegetation,  leafbrown = brownvegetation,
-                     typelidf,  lidfa,  lidfb,  lai,  q,  tts,  tto,  psi,  rsoil,
-                     fraction_brown,  diss,  cv,  zeta)
-  }
-  return(ref)
-}
-
-#" Performs PROSAIL simulation based on a set of combinations of input parameters
-#" @param leafoptics list. Includes leaf optical properties (reflectance and transmittance)
-#" and corresponding spectral bands
-#" @param typelidf numeric. Type of leaf inclination distribution function
-#" @param lidfa numeric.
-#" if typelidf  == 1,  controls the average leaf slope
-#" if typelidf  == 2,  corresponds to average leaf angle
-#" @param lidfb numeric.
-#" if typelidf  == 1,  unused
-#" if typelidf  == 2,  controls the distribution"s bimodality
-#" @param lai numeric. Leaf Area Index
-#" @param q numeric. Hot Spot parameter
-#" @param tts numeric. Sun zeith angle
-#" @param tto numeric. Observer zeith angle
-#" @param psi numeric. Azimuth Sun  /  Observer
-#" @param rsoil numeric. Soil reflectance
-#"
-#" @return list. rdot, rsot, rddt, rsdt
-#" rdot: hemispherical-directional reflectance factor in viewing direction
-#" rsot: bi-directional reflectance factor
-#" rsdt: directional-hemispherical reflectance factor for solar incident flux
-#" rddt: bi-hemispherical reflectance factor
-#" @export
-
-foursail  <- function(leafoptics,  typelidf = 2,  lidfa = nULL,  lidfb = nULL,  lai = nULL,
-                      q = nULL,  tts = nULL,  tto = nULL,  psi = nULL,  rsoil = nULL) {
-
-  ############################## #
-  #	LEAF OPTICAL PROPERTIES	##
-  ############################## #
-  rho <- leafoptics$Reflectance
-  tau <- leafoptics$Transmittance
-
-  #	Geometric quantities
-  rd <- pi / 180
-  cts <- cos(rd * tts)
-  cto <- cos(rd * tto)
-  ctscto <- cts * cto
-  tants <- tan(rd * tts)
-  tanto <- tan(rd * tto)
-  cospsi <- cos(rd * psi)
-  dso <- sqrt(tants * tants + tanto * tanto - 2. * tants * tanto * cospsi)
-
-  #	Generate leaf angle distribution from average leaf angle (ellipsoidal) or (a, b) parameters
-  if (typelidf == 1) {
-    foliar_distrib <- dladgen(lidfa, lidfb)
-    lidf <- foliar_distrib$lidf
-    litab <- foliar_distrib$litab
-
-  } else if (typelidf == 2) {
-    foliar_distrib <- campbell(lidfa)
-    lidf <- foliar_distrib$lidf
-    litab <- foliar_distrib$litab
-  }
-
-  # angular distance,  compensation of shadow length
-  #	Calculate geometric factors associated with extinction and scattering
-  #	Initialise sums
-  ks <- 0
-  ko <- 0
-  bf <- 0
-  sob <- 0
-  sof <- 0
-
-  #	Weighted sums over LIDF
-  na <- length(litab)
-  for (i in 1:na) {
-    ttl <- litab[i]	    # leaf inclination discrete values
-    ctl <- cos(rd * ttl)
-    #	SAIL volume scattering phase function gives interception and portions to be
-    #	multiplied by rho and tau
-    resvolscatt <- volscatt(tts, tto, psi, ttl)
-    chi_s <- resvolscatt$chi_s
-    chi_o <- resvolscatt$chi_o
-    frho <- resvolscatt$frho
-    ftau <- resvolscatt$ftau
-
-    #********************************************************************************
-    #*                   SUITS SYSTEM coEFFICIEnTS
-    #*
-    #*	ks  : Extinction coefficient for direct solar flux
-    #*	ko  : Extinction coefficient for direct observed flux
-    #*	att : Attenuation coefficient for diffuse flux
-    #*	sigb : Backscattering coefficient of the diffuse downward flux
-    #*	sigf : Forwardscattering coefficient of the diffuse upward flux
-    #*	sf  : Scattering coefficient of the direct solar flux for downward diffuse flux
-    #*	sb  : Scattering coefficient of the direct solar flux for upward diffuse flux
-    #*	vf   : Scattering coefficient of upward diffuse flux in the observed direction
-    #*	vb   : Scattering coefficient of downward diffuse flux in the observed direction
-    #*	w   : Bidirectional scattering coefficient
-    #********************************************************************************
-
-    #	Extinction coefficients
-    ksli <- chi_s / cts
-    koli <- chi_o / cto
-
-    #	Area scattering coefficient fractions
-    sobli <- frho * pi / ctscto
-    sofli <- ftau * pi / ctscto
-    bfli <- ctl * ctl
-    ks <- ks + ksli * lidf[i]
-    ko <- ko + koli * lidf[i]
-    bf <- bf + bfli * lidf[i]
-    sob <- sob + sobli * lidf[i]
-    sof <- sof + sofli * lidf[i]
-  }
-
-  #	Geometric factors to be used later with rho and tau
-  sdb <- 0.5 * (ks + bf)
-  sdf <- 0.5 * (ks - bf)
-  dob <- 0.5 * (ko + bf)
-  dof <- 0.5 * (ko - bf)
-  ddb <- 0.5 * (1. + bf)
-  ddf <- 0.5 * (1. - bf)
-
-  #	Here rho and tau come in
-  sigb <- ddb * rho + ddf * tau
-  sigf <- ddf * rho + ddb * tau
-  att <- 1 - sigf
-  m2 <- (att + sigb) * (att - sigb)
-  m2[which(m2 <= 0)] <- 0
-  m <- sqrt(m2)
-
-  sb <- sdb * rho + sdf * tau
-  sf <- sdf * rho + sdb * tau
-  vb <- dob * rho + dof * tau
-  vf <- dof * rho + dob * tau
-  w <- sob * rho + sof * tau
-
-  #	Here the LAI comes in
-  #   Outputs for the case LAI = 0
-  if (lai < 0) {
-    tss <- 1
-    too <- 1
-    tsstoo <- 1
-    rdd <- 0
-    tdd <- 1
-    rsd <- 0
-    tsd <- 0
-    rdo <- 0
-    tdo <- 0
-    rso <- 0
-    rsos <- 0
-    rsod <- 0
-
-    rddt <- rsoil
-    rsdt <- rsoil
-    rdot <- rsoil
-    rsodt <- 0 * rsoil
-    rsost <- rsoil
-    rsot <- rsoil
-  } else {
-    #	Other cases (LAI  >  0)
-    e1 <- exp(-m * lai)
-    e2 <- e1 * e1
-    rinf <- (att - m) / sigb
-    rinf2 <- rinf * rinf
-    re <- rinf * e1
-    denom <- 1. - rinf2 * e2
-
-    j1ks <- jfunc1(ks, m, lai)
-    j2ks <- jfunc2(ks, m, lai)
-    j1ko <- jfunc1(ko, m, lai)
-    j2ko <- jfunc2(ko, m, lai)
-
-    ps <- (sf + sb * rinf) * j1ks
-    qs <- (sf * rinf + sb) * j2ks
-    pv <- (vf + vb * rinf) * j1ko
-    qv <- (vf * rinf + vb) * j2ko
-
-    rdd <- rinf * (1. - e2) / denom
-    tdd <- (1. - rinf2) * e1 / denom
-    tsd <- (ps - re * qs) / denom
-    rsd <- (qs - re * ps) / denom
-    tdo <- (pv - re * qv) / denom
-    rdo <- (qv - re * pv) / denom
-
-    tss <- exp(-ks * lai)
-    too <- exp(-ko * lai)
-    z <- jfunc3(ks, ko, lai)
-    g1 <- (z - j1ks * too) / (ko + m)
-    g2 <- (z - j1ko * tss) / (ks + m)
-
-    tv1 <- (vf * rinf + vb) * g1
-    tv2 <- (vf + vb * rinf) * g2
-    t1 <- tv1 * (sf + sb * rinf)
-    t2 <- tv2 * (sf * rinf + sb)
-    t3 <- (rdo * qs + tdo * ps) * rinf
-
-    #	Multiple scattering contribution to bidirectional canopy reflectance
-    rsod <- (t1 + t2 - t3) / (1. - rinf2)
-
-    #	Treatment of the hotspot-effect
-    alf <- 1e6
-    #	Apply correction 2 / (K + k) suggested by F.-M. Breon
-    if (q > 0) {
-      alf <- (dso / q) * 2. / (ks + ko)
-    }
-    if (alf > 200) {
-      # inserted H. Bach 1 / 3 / 04
-      alf <- 200
-    }
-    if (alf == 0) {
-      #	The pure hotspot - no shadow
-      tsstoo <- tss
-      sumint <- (1 - tss) / (ks * lai)
-    } else {
-      #	Outside the hotspot
-      fhot <- lai * sqrt(ko * ks)
-      #	Integrate by exponential Simpson method in 20 steps
-      #	the steps are arranged according to equal partitioning
-      #	of the slope of the joint probability function
-      x1 <- 0
-      y1 <- 0
-      f1 <- 1
-      fint <- (1. - exp(-alf)) * 0.05
-      sumint <- 0
-      for (i in 1:20) {
-        if (i < 20) {
-          x2 <- -log(1. - i * fint) / alf
-        } else {
-          x2 <- 1
-        }
-        y2 <- -(ko + ks) * lai * x2 + fhot * (1. - exp(-alf * x2)) / alf
-        f2 <- exp(y2)
-        sumint <- sumint + (f2 - f1) * (x2 - x1) / (y2 - y1)
-        x1 <- x2
-        y1 <- y2
-        f1 <- f2
-      }
-      tsstoo <- f1
-    }
-    #	Bidirectional reflectance
-    #	Single scattering contribution
-    rsos <- w * lai * sumint
-    #	Total canopy contribution
-    rso <- rsos + rsod
-    #	Interaction with the soil
-    dn <- 1. - rsoil * rdd
-    # rddt: bi-hemispherical reflectance factor
-    rddt <- rdd + tdd * rsoil * tdd / dn
-    # rsdt: directional-hemispherical reflectance factor for solar incident flux
-    rsdt <- rsd + (tsd + tss) * rsoil * tdd / dn
-    # rdot: hemispherical-directional reflectance factor in viewing direction
-    rdot <- rdo + tdd * rsoil * (tdo + too) / dn
-    # rsot: bi-directional reflectance factor
-    rsodt <- rsod + ((tss + tsd) * tdo + (tsd + tss * rsoil * rdd) * too) * rsoil / dn
-    rsost <- rsos + tsstoo * rsoil
-    rsot <- rsost + rsodt
-  }
-  my_list <- list("rdot" = rdot, "rsot" = rsot, "rddt" = rddt, "rsdt" = rsdt)
-  return(my_list)
-}
-
-#" Performs pro4sail2 simulation based on a set of combinations of input parameters
-#" @param leafgreen list. includes relfectance and transmittance for vegetation #1 (e.g. green vegetation)
-#" @param leafbrown list. includes relfectance and transmittance for vegetation #2 (e.g. brown vegetation)
-#" @param typelidf numeric. Type of leaf inclination distribution function
-#" @param lidfa numeric.
-#" if typelidf  == 1,  controls the average leaf slope
-#" if typelidf  == 2,  corresponds to average leaf angle
-#" @param lidfb numeric.
-#" if typelidf  == 1,  unused
-#" if typelidf  == 2,  controls the distribution"s bimodality
-#" @param lai numeric. Leaf Area Index
-#" @param hot numeric. Hot Spot parameter = ratio of the correlation length of leaf projections in the horizontal plane and the canopy height (doi:10.1016 / j.rse.2006.12.013)
-#" @param tts numeric. Sun zeith angle
-#" @param tto numeric. Observer zeith angle
-#" @param psi numeric. Azimuth Sun  /  Observer
-#" @param rsoil numeric. Soil reflectance
-#" @param fraction_brown numeric. Fraction of brown leaf area
-#" @param diss numeric. Layer dissociation factor
-#" @param cv numeric. vertical crown cover percentage
-#" = % ground area covered with crowns as seen from nadir direction
-#" @param zeta numeric. Tree shape factor
-#" = ratio of crown diameter to crown height
-#"
-#" @return list. rdot, rsot, rddt, rsdt
-#" rdot: hemispherical-directional reflectance factor in viewing direction
-#" rsot: bi-directional reflectance factor
-#" rsdt: directional-hemispherical reflectance factor for solar incident flux
-#" rddt: bi-hemispherical reflectance factor
-#" alfast: canopy absorptance for direct solar incident flux
-#" alfadt: canopy absorptance for hemispherical diffuse incident flux
-#" @export
-
-foursail2  <- function(leafgreen,  leafbrown,
-                       typelidf = 2, lidfa = nULL, lidfb = nULL,
-                       lai = nULL,  hot = nULL, tts = nULL, tto = nULL, psi = nULL, rsoil = nULL,
-                       fraction_brown = 0.5,  diss = 0.5,  cv = 1, zeta = 1) {
-
-  #	This version does not include non-Lambertian soil properties.
-  #	original codes do,  and only need to add the following variables as input
-  rddsoil <- rdosoil <- rsdsoil <- rsosoil <- rsoil
-
-  #	Geometric quantities
-  rd <- pi / 180
-
-  #	Generate leaf angle distribution from average leaf angle (ellipsoidal) or (a, b) parameters
-  if (typelidf == 1) {
-    foliar_distrib <- dladgen(lidfa, lidfb)
-    lidf <- foliar_distrib$lidf
-    litab <- foliar_distrib$litab
-
-  } else if (typelidf == 2) {
-    foliar_distrib <- campbell(lidfa)
-    lidf <- foliar_distrib$lidf
-    litab <- foliar_distrib$litab
-  }
-
-  if (lai < 0) {
-    message("Please define positive LAI value")
-    rddt <- rsdt <- rdot <- rsost <- rsot <- rsoil
-    alfast <- alfadt <- 0 * rsoil
-  } else if (lai == 0) {
-    tss <- too <- tsstoo <- tdd <- 1.0
-    rdd <- rsd <- tsd <- rdo <- tdo <- 0.0
-    rso <- rsos <- rsod <- rsodt <- 0.0
-    rddt <- rsdt <- rdot <- rsost <- rsot <- rsoil
-    alfast <- alfadt <- 0 * rsoil
-  } else if (lai > 0) {
-    cts <- cos(rd * tts)
-    cto <- cos(rd * tto)
-    ctscto <- cts * cto
-    tants <- tan(rd * tts)
-    tanto <- tan(rd * tto)
-    cospsi <- cos(rd * psi)
-    dso <- sqrt(tants * tants + tanto * tanto - 2.0 * tants * tanto * cospsi)
-
-    # Clumping effects
-    cs <- co <- 1.0
-    if (cv <= 1.0) {
-      cs <- 1.0 - (1.0 - cv)^(1.0 / cts)
-      co <- 1.0 - (1.0 - cv)^(1.0 / cto)
-    }
-    overlap <- 0.0
-    if (zeta > 0.0) {
-      overlap <- min(cs * (1.0 - co), co * (1.0 - cs)) * exp(-dso / zeta)
-    }
-    fcd <- cs * co + overlap
-    fcs <- (1.0 - cs) * co - overlap
-    fod <- cs * (1.0 - co) - overlap
-    fos <- (1.0 - cs) * (1.0 - co) + overlap
-    fcdc <- 1.0 - (1.0 - fcd)^(0.5 / cts + 0.5 / cto)
-
-    #	Part depending on diss,  fraction_brown,  and leaf optical properties
-    #	First save the input fraction_brown as the old fraction_brown,  as the following change is only artificial
-    # Better define an fraction_brown that is actually used: fb,  so that the input is not modified!
-
-    fb <- fraction_brown
-    # if only green leaves
-    if (fraction_brown == 0.0) {
-      fb <- 0.5
-      leafbrown$Reflectance <- leafgreen$Reflectance
-      leafbrown$Transmittance <- leafgreen$Transmittance
-    }
-    if (fraction_brown == 1.0) {
-      fb <- 0.5
-      leafgreen$Reflectance <- leafbrown$Reflectance
-      leafgreen$Transmittance <- leafbrown$Transmittance
-    }
-    s <- (1.0 - diss) * fb * (1.0 - fb)
-    # rho1 && tau1 : green foliage
-    # rho2 && tau2 : brown foliage (bottom layer)
-    rho1 <- ((1 - fb - s) * leafgreen$Reflectance + s * leafbrown$Reflectance) / (1 - fb)
-    tau1 <- ((1 - fb - s) * leafgreen$Transmittance + s * leafbrown$Transmittance) / (1 - fb)
-    rho2 <- (s * leafgreen$Reflectance + (fb - s) * leafbrown$Reflectance) / fb
-    tau2 <- (s * leafgreen$Transmittance + (fb - s) * leafbrown$Transmittance) / fb
-
-    # angular distance,  compensation of shadow length
-    #	Calculate geometric factors associated with extinction and scattering
-    #	Initialise sums
-    ks <- ko <- bf <- sob <- sof <- 0
-
-    # Weighted sums over LIDF
-
-    for (i in 1:seq_along(litab)) {
-      ttl <- litab[i]
-      ctl <- cos(rd * ttl)
-      # SAIL volscatt function gives interception coefficients
-      # and two portions of the volume scattering phase function to be
-      # multiplied by rho and tau,  respectively
-      resvolscatt <- volscatt(tts, tto, psi, ttl)
-      chi_s <- resvolscatt$chi_s
-      chi_o <- resvolscatt$chi_o
-      frho <- resvolscatt$frho
-      ftau <- resvolscatt$ftau
-      # Extinction coefficients
-      ksli <- chi_s / cts
-      koli <- chi_o / cto
-      # Area scattering coefficient fractions
-      sobli <- frho * pi / ctscto
-      sofli <- ftau * pi / ctscto
-      bfli <- ctl * ctl
-      ks <- ks + ksli * lidf[i]
-      ko <- ko + koli * lidf[i]
-      bf <- bf + bfli * lidf[i]
-      sob <- sob + sobli * lidf[i]
-      sof <- sof + sofli * lidf[i]
-    }
-    # Geometric factors to be used later in combination with rho and tau
-    sdb <- 0.5 * (ks + bf)
-    sdf <- 0.5 * (ks - bf)
-    dob <- 0.5 * (ko + bf)
-    dof <- 0.5 * (ko - bf)
-    ddb <- 0.5 * (1. + bf)
-    ddf <- 0.5 * (1. - bf)
-
-    # LAIs in two layers
-    lai1 <- (1 - fb) * lai
-    lai2 <- fb * lai
-
-    tss <- exp(-ks * lai)
-    ck <- exp(-ks * lai1)
-    alf <- 1e6
-    if (hot > 0.0) {
-      alf <- (dso / hot) * 2.0 / (ks + ko)
-    }
-    if (alf > 200.0) {
-      alf <- 200.0     # inserted H. Bach 1 / 3 / 04
-    }
-    if (alf == 0.0) {
-      # The pure hotspot
-      tsstoo <- tss
-      s1 <- (1 - ck) / (ks * lai)
-      s2 <- (ck - tss) / (ks * lai)
-    } else {
-      # Outside the hotspot
-      fhot <- lai * sqrt(ko * ks)
-      # Integrate 2 layers by exponential simpson method in 20 steps
-      # the steps are arranged according to equal partitioning
-      # of the derivative of the joint probability function
-      x1 <- y1 <- 0.0
-      f1 <- 1.0
-      ca <- exp(alf * (fb - 1.0))
-      fint <- (1.0 - ca) * .05
-      s1 <- 0.0
-      for (istep in 1:20) {
-        if (istep < 20) {
-          x2 <- -log(1. - istep * fint) / alf
-        } else {
-          x2 <- 1. - fb
-        }
-        y2 <- -(ko + ks) * lai * x2 + fhot * (1.0 - exp(-alf * x2)) / alf
-        f2 <- exp(y2)
-        s1 <- s1 + (f2 - f1) * (x2 - x1) / (y2 - y1)
-        x1 <- x2
-        y1 <- y2
-        f1 <- f2
-      }
-      fint <- (ca - exp(-alf)) * .05
-      s2 <- 0.0
-      for (istep in 1:20) {
-        if (istep < 20) {
-          x2 <- -log(ca - istep * fint) / alf
-        } else {
-          x2 <- 1.0
-        }
-        y2 <- -(ko + ks) * lai * x2 + fhot * (1.0 - exp(-alf * x2)) / alf
-        f2 <- exp(y2)
-        s2 <- s2 + (f2 - f1) * (x2 - x1) / (y2 - y1)
-        x1 <- x2
-        y1 <- y2
-        f1 <- f2
-      }
-      tsstoo <- f1
-    }
-
-    # Calculate reflectances and transmittances
-    # Bottom layer
-    tss <- exp(-ks * lai2)
-    too <- exp(-ko * lai2)
-    sb <- sdb * rho2 + sdf * tau2
-    sf <- sdf * rho2 + sdb * tau2
-
-    vb <- dob * rho2 + dof * tau2
-    vf <- dof * rho2 + dob * tau2
-
-    w2 <- sob * rho2 + sof * tau2
-
-    sigb <- ddb * rho2 + ddf * tau2
-    sigf <- ddf * rho2 + ddb * tau2
-    att <- 1.0 - sigf
-    m2 <- (att + sigb) * (att - sigb)
-    m2[m2 < 0] <- 0
-    m <- sqrt(m2)
-    which_ncs <- which(m > 0.01)
-    which_cs <- which(m <= 0.01)
-
-    tdd <- rdd <- tsd <- rsd <- tdo <- rdo <- 0 * m
-    rsod <- 0 * m
-    if (length(which_ncs) > 0) {
-      resncs <- nonconservativescattering(m[which_ncs], lai2, att[which_ncs], sigb[which_ncs],
-                                          ks, ko, sf[which_ncs], sb[which_ncs], vf[which_ncs], vb[which_ncs], tss, too)
-      tdd[which_ncs] <- resncs$tdd
-      rdd[which_ncs] <- resncs$rdd
-      tsd[which_ncs] <- resncs$tsd
-      rsd[which_ncs] <- resncs$rsd
-      tdo[which_ncs] <- resncs$tdo
-      rdo[which_ncs] <- resncs$rdo
-      rsod[which_ncs] <- resncs$rsod
-    }
-    if (length(which_cs) > 0) {
-      rescs <- conservativescattering(m[which_cs], lai2, att[which_cs], sigb[which_cs],
-                                      ks, ko, sf[which_cs], sb[which_cs], vf[which_cs], vb[which_cs], tss, too)
-      tdd[which_cs] <- rescs$tdd
-      rdd[which_cs] <- rescs$rdd
-      tsd[which_cs] <- rescs$tsd
-      rsd[which_cs] <- rescs$rsd
-      tdo[which_cs] <- rescs$tdo
-      rdo[which_cs] <- rescs$rdo
-      rsod[which_cs] <- rescs$rsod
-    }
-
-    # Set background properties equal to those of the bottom layer on a black soil
-    rddb <- rdd
-    rsdb <- rsd
-    rdob <- rdo
-    rsodb <- rsod
-    tddb <- tdd
-    tsdb <- tsd
-    tdob <- tdo
-    toob <- too
-    tssb <- tss
-    # Top layer
-    tss <- exp(-ks * lai1)
-    too <- exp(-ko * lai1)
-
-    sb <- sdb * rho1 + sdf * tau1
-    sf <- sdf * rho1 + sdb * tau1
-
-    vb <- dob * rho1 + dof * tau1
-    vf <- dof * rho1 + dob * tau1
-
-    w1 <- sob * rho1 + sof * tau1
-
-    sigb <- ddb * rho1 + ddf * tau1
-    sigf <- ddf * rho1 + ddb * tau1
-    att <- 1.0 - sigf
-
-    m2 <- (att + sigb) * (att - sigb)
-    m2[m2 < 0] <- 0
-    m <- sqrt(m2)
-    which_ncs <- which(m > 0.01)
-    which_cs <- which(m <= 0.01)
-
-    tdd <- rdd <- tsd <- rsd <- tdo <- rdo <- 0 * m
-    rsod <- 0 * m
-    if (length(which_ncs) > 0) {
-      resncs <- nonconservativescattering(m[which_ncs], lai1, att[which_ncs], sigb[which_ncs],
-                                          ks, ko, sf[which_ncs], sb[which_ncs], vf[which_ncs], vb[which_ncs], tss, too)
-      tdd[which_ncs] <- resncs$tdd
-      rdd[which_ncs] <- resncs$rdd
-      tsd[which_ncs] <- resncs$tsd
-      rsd[which_ncs] <- resncs$rsd
-      tdo[which_ncs] <- resncs$tdo
-      rdo[which_ncs] <- resncs$rdo
-      rsod[which_ncs] <- resncs$rsod
-    }
-    if (length(which_cs) > 0) {
-      rescs <- conservativescattering(m[which_cs], lai1, att[which_cs], sigb[which_cs],
-                                      ks, ko, sf[which_cs], sb[which_cs], vf[which_cs], vb[which_cs], tss, too)
-      tdd[which_cs] <- rescs$tdd
-      rdd[which_cs] <- rescs$rdd
-      tsd[which_cs] <- rescs$tsd
-      rsd[which_cs] <- rescs$rsd
-      tdo[which_cs] <- rescs$tdo
-      rdo[which_cs] <- rescs$rdo
-      rsod[which_cs] <- rescs$rsod
-    }
-
-    # combine with bottom layer reflectances and transmittances (adding method)
-    rn <- 1.0 - rdd * rddb
-    tup <- (tss * rsdb + tsd * rddb) / rn
-    tdn <- (tsd + tss * rsdb * rdd) / rn
-    rsdt <- rsd + tup * tdd
-    rdot <- rdo + tdd * (rddb * tdo + rdob * too) / rn
-    rsodt <- rsod + (tss * rsodb + tdn * rdob) * too + tup * tdo
-
-    rsost <- (w1 * s1 + w2 * s2) * lai
-
-    rsot <- rsost + rsodt
-
-    # Diffuse reflectances at the top and the bottom are now different
-    rddt_t <- rdd + tdd * rddb * tdd / rn
-    rddt_b <- rddb + tddb * rdd * tddb / rn
-
-    # Transmittances of the combined canopy layers
-    tsst <- tss * tssb
-    toot <- too * toob
-    tsdt <- tss * tsdb + tdn * tddb
-    tdot <- tdob * too + tddb * (tdo + rdd * rdob * too) / rn
-    tddt <- tdd * tddb / rn
-
-    # Apply clumping effects to vegetation layer
-    rddcb <- cv * rddt_b
-    rddct <- cv * rddt_t
-    tddc <- 1 - cv + cv * tddt
-    rsdc <- cs * rsdt
-    tsdc <- cs * tsdt
-    rdoc <- co * rdot
-    tdoc <- co * tdot
-    tssc <- 1 - cs + cs * tsst
-    tooc <- 1 - co + co * toot
-
-    # new weight function fcdc for crown contribution (W. Verhoef,  22-05-08)
-    rsoc <- fcdc * rsot
-    tssooc <- fcd * tsstoo + fcs * toot + fod * tsst + fos
-    # Canopy absorptance for black background (W. Verhoef,  02-03-04)
-    alfas <- 1. - tssc - tsdc - rsdc
-    alfad <- 1. - tddc - rddct
-    # Add the soil background
-    rn <- 1 - rddcb * rddsoil
-    tup <- (tssc * rsdsoil + tsdc * rddsoil) / rn
-    tdn <- (tsdc + tssc * rsdsoil * rddcb) / rn
-
-    rddt <- rddct + tddc * rddsoil * tddc / rn
-    rsdt <- rsdc + tup * tddc
-    rdot <- rdoc + tddc * (rddsoil * tdoc + rdosoil * tooc) / rn
-    rsot <- rsoc + tssooc * rsosoil + tdn * rdosoil * tooc + tup * tdoc
-
-    # Effect of soil background on canopy absorptances (W. Verhoef,  02-03-04)
-    alfast <- alfas + tup * alfad
-    alfadt <- alfad * (1. + tddc * rddsoil / rn)
-  }
-  my_list <- list("rdot" = rdot, "rsot" = rsot, "rddt" = rddt, "rsdt" = rsdt,
-                  "alfast" = alfast,  "alfadt" = alfadt)
-  return(my_list)
-}
-
-
-
-#" computes non conservative scattering conditions
-#" @param m numeric.
-#" @param lai numeric. Leaf Area Index
-#" @param att numeric.
-#" @param sigb numeric.
-#" @param ks numeric.
-#" @param ko numeric.
-#" @param sf numeric.
-#" @param sb numeric.
-#" @param vf numeric.
-#" @param vb numeric.
-#" @param tss numeric.
-#" @param too numeric.
-#"
-#" @return list. tdd,  rdd,  tsd,  rsd,  tdo,  rdo,  rsod
-#"
-#" @export
-nonconservativescattering <- function(m, lai, att, sigb, ks, ko, sf, sb, vf, vb, tss, too) {
-
-  e1 <- exp(-m * lai)
-  e2 <- e1 * e1
-  rinf <- (att - m) / sigb
-  rinf2 <- rinf * rinf
-  re <- rinf * e1
-  denom <- 1. - rinf2 * e2
-
-  j1ks <- jfunc1(ks, m, lai)
-  j2ks <- jfunc2(ks, m, lai)
-  j1ko <- jfunc1(ko, m, lai)
-  j2ko <- jfunc2(ko, m, lai)
-
-  ps <- (sf + sb * rinf) * j1ks
-  qs <- (sf * rinf + sb) * j2ks
-  pv <- (vf + vb * rinf) * j1ko
-  qv <- (vf * rinf + vb) * j2ko
-
-  tdd <- (1. - rinf2) * e1 / denom
-  rdd <- rinf * (1. - e2) / denom
-  tsd <- (ps - re * qs) / denom
-  rsd <- (qs - re * ps) / denom
-  tdo <- (pv - re * qv) / denom
-  rdo <- (qv - re * pv) / denom
-
-  z <- jfunc2(ks, ko, lai)
-  g1 <- (z - j1ks * too) / (ko + m)
-  g2 <- (z - j1ko * tss) / (ks + m)
-
-  tv1 <- (vf * rinf + vb) * g1
-  tv2 <- (vf + vb * rinf) * g2
-
-  t1 <- tv1 * (sf + sb * rinf)
-  t2 <- tv2 * (sf * rinf + sb)
-  t3 <- (rdo * qs + tdo * ps) * rinf
-
-  # Multiple scattering contribution to bidirectional canopy reflectance
-  rsod <- (t1 + t2 - t3) / (1. - rinf2)
-  my_list <- list("tdd" = tdd,  "rdd" = rdd,  "tsd" = tsd,
-                  "rsd" = rsd,  "tdo" = tdo,  "rdo" = rdo,
-                  "rsod" = rsod)
-  return(my_list)
-}
-
-#" computes conservative scattering conditions
-#" @param m numeric.
-#" @param lai numeric. Leaf Area Index
-#" @param att numeric.
-#" @param sigb numeric.
-#" @param ks numeric.
-#" @param ko numeric.
-#" @param sf numeric.
-#" @param sb numeric.
-#" @param vf numeric.
-#" @param vb numeric.
-#" @param tss numeric.
-#" @param too numeric.
-#"
-#" @return list. tdd,  rdd,  tsd,  rsd,  tdo,  rdo,  rsod
-#"
-#" @export
-conservativescattering <- function(m, lai, att, sigb, ks, ko, sf, sb, vf, vb, tss, too) {
-
-  # near or complete conservative scattering
-  j4 <- jfunc4(m, lai)
-  amsig <- att - sigb
-  apsig <- att + sigb
-  rtp <- (1 - amsig * j4) / (1 + amsig * j4)
-  rtm <- (-1 + apsig * j4) / (1 + apsig * j4)
-  rdd <- 0.5 * (rtp + rtm)
-  tdd <- 0.5 * (rtp - rtm)
-
-  dns <- ks * ks - m * m
-  dno <- ko * ko - m * m
-  cks <- (sb * (ks - att) - sf * sigb) / dns
-  cko <- (vb * (ko - att) - vf * sigb) / dno
-  dks <- (-sf * (ks + att) - sb * sigb) / dns
-  dko <- (-vf * (ko + att) - vb * sigb) / dno
-  ho <- (sf * cko + sb * dko) / (ko + ks)
-
-  rsd <- cks * (1 - tss * tdd) - dks * rdd
-  rdo <- cko * (1 - too * tdd) - dko * rdd
-  tsd <- dks * (tss - tdd) - cks * tss * rdd
-  tdo <- dko * (too - tdd) - cko * too * rdd
-  # Multiple scattering contribution to bidirectional canopy reflectance
-  rsod <- ho * (1 - tss * too) - cko * tsd * too - dko * rsd
-
-  my_list <- list("tdd" = tdd,  "rdd" = rdd,  "tsd" = tsd,
-                  "rsd" = rsd,  "tdo" = tdo,  "rdo" = rdo,
-                  "rsod" = rsod)
-  return(my_list)
-}
-
-
-
-
-
-
-#" computes the leaf angle distribution function value (freq)
-#"
-#" Ellipsoidal distribution function characterised by the average leaf
-#" inclination angle in degree (ala)
-#" Campbell 1986
-#" @param ala average leaf angle
-#" @return foliar_distrib list. lidf and litab
-#" @export
-campbell  <- function(ala) {
-
-  tx1 <- c(10., 20., 30., 40., 50., 60., 70., 80., 82., 84., 86., 88., 90.)
-  tx2 <- c(0., 10., 20., 30., 40., 50., 60., 70., 80., 82., 84., 86., 88.)
-  litab <- (tx2 + tx1) / 2
-  n <- length(litab)
-  tl1 <- tx1 * (pi / 180)
-  tl2 <- tx2 * (pi / 180)
-  excent <- exp(-1.6184e-5 * ala**3 + 2.1145e-3 * ala**2 - 1.2390e-1 * ala + 3.2491)
-  sum0 <- 0
-
-  freq <- c()
-  for (i in 1:n) {
-    x1 <- excent / (sqrt(1. + excent**2. * tan(tl1[i])**2))
-    x2 <- excent / (sqrt(1. + excent**2. * tan(tl2[i])**2))
-    if (excent == 1) {
-      freq[i] <- abs(cos(tl1[i]) - cos(tl2[i]))
-    } else {
-      alpha <- excent / sqrt(abs(1 - excent**2))
-      alpha2 <- alpha**2
-      x12 <- x1**2
-      x22 <- x2**2
-      alpx1 <- 0 * alpha2
-      alpx2 <- 0 * alpha2
-      almx1 <- 0 * alpha2
-      almx2 <- 0 * alpha2
-      if (excent > 1) {
-        alpx1 <- sqrt(alpha2[excent > 1] + x12[excent > 1])
-        alpx2[excent > 1] <- sqrt(alpha2[excent > 1] + x22[excent > 1])
-        dum <- x1 * alpx1 + alpha2 * log(x1 + alpx1)
-        freq[i] <- abs(dum - (x2 * alpx2 + alpha2 * log(x2 + alpx2)))
-      } else {
-        almx1 <- sqrt(alpha2 - x12)
-        almx2 <- sqrt(alpha2 - x22)
-        dum <- x1 * almx1 + alpha2 * asin(x1 / alpha)
-        freq[i] <- abs(dum - (x2 * almx2 + alpha2 * asin(x2 / alpha)))
-      }
-    }
-  }
-  sum0 <- sum(freq)
-  freq0 <- freq / sum0
-  foliar_distrib <- list("lidf" = freq0, "litab" = litab)
-  return(foliar_distrib)
-}
-
-#" computes the leaf angle distribution function value (freq)
-#"
-#" Using the original bimodal distribution function initially proposed in SAIL
-#"  references
-#"  ----------
-#"  (Verhoef1998) Verhoef,  Wout. Theory of radiative transfer models applied
-#"  in optical remote sensing of vegetation canopies.
-#"  nationaal Lucht en Ruimtevaartlaboratorium,  1998.
-#"  http: /  / library.wur.nl / WebQuery / clc / 945481.
-#" @param a controls the average leaf slope
-#" @param b controls the distribution"s bimodality
-#" LIDF type 		  a 		b
-#" Planophile 	  1		  0
-#" Erectophile    -1	 	0
-#" Plagiophile 	  0		  -1
-#" Extremophile 	0		  1
-#" Spherical 	    -0.35 -0.15
-#" Uniform        0     0
-#" requirement: ||lidfa||  +  ||lidfb|| < 1
-#"
-#" @return foliar_distrib list. lidf and litab
-#" @export
-dladgen  <- function(a, b) {
-  litab <- c(5., 15., 25., 35., 45., 55., 65., 75., 81., 83., 85., 87., 89.)
-  freq <- c()
-  for (i1 in 1:8) {
-    t <- i1 * 10
-    freq[i1] <- dcum(a, b, t)
-  }
-  for (i2 in 9:12) {
-    t <- 80. + (i2 - 8) * 2.
-    freq[i2] <- dcum(a, b, t)
-  }
-  freq[13] <- 1
-  for (i in 13:2) {
-    freq[i] <- freq[i] - freq[i - 1]
-  }
-  foliar_distrib <- list("lidf" = freq, "litab" = litab)
-  return(foliar_distrib)
-}
-
-#" dcum function
-#" @param a numeric. controls the average leaf slope
-#" @param b numeric. controls the distribution"s bimodality
-#" @param t numeric. angle
-#" @return f
-#" @export
-dcum <- function(a, b, t) {
-  rd <- pi / 180
-  if (a >= 1) {
-    f <- 1 - cos(rd * t)
-  } else {
-    eps <- 1e-8
-    delx <- 1
-    x <- 2 * rd * t
-    p <- x
-    while (delx  >= eps) {
-      y <- a * sin(x) + .5 * b * sin(2. * x)
-      dx <- .5 * (y - x + p)
-      x <- x + dx
-      delx <- abs(dx)
-    }
-    f <- (2. * y + p) / pi
-  }
-  return(f)
-}
-
-#" J1 function with avoidance of singularity problem
-#"
-#" @param k numeric. Extinction coefficient for direct (solar or observer) flux
-#" @param l numeric.
-#" @param t numeric. Leaf Area Index
-#" @return jout numeric.
-#" @export
-jfunc1 <- function(k, l, t) {
-  # J1 function with avoidance of singularity problem
-  del <- (k - l) * t
-  jout <- 0 * l
-  jout[which(abs(del) > 1e-3)] <- (exp(-l[which(abs(del) > 1e-3)] * t) - exp(-k * t)) / (k - l[which(abs(del) > 1e-3)])
-  jout[which(abs(del) <= 1e-3)] <- 0.5 * t * (exp(-k * t) + exp(-l[which(abs(del) <= 1e-3)] * t)) * (1 - del[which(abs(del) <= 1e-3)] * del[which(abs(del) <= 1e-3)] / 12)
-  return(jout)
-}
-
-#" J2 function with avoidance of singularity problem
-#"
-#" @param k numeric. Extinction coefficient for direct (solar or observer) flux
-#" @param l numeric.
-#" @param t numeric. Leaf Area Index
-#" @return jout numeric.
-#" @export
-jfunc2 <- function(k, l, t) {
-  #	J2 function
-  jout <- (1. - exp(-(k + l) * t)) / (k + l)
-  return(jout)
-}
-
-#" J3 function with avoidance of singularity problem
-#"
-#" @param k numeric. Extinction coefficient for direct (solar or observer) flux
-#" @param l numeric.
-#" @param t numeric. Leaf Area Index
-#" @return jout numeric.
-#" @export
-jfunc3 <- function(k, l, t) {
-  out <- (1. - exp(-(k + l) * t)) / (k + l)
-  return(out)
-}
-
-
-#" j4 function for treating (near) conservative scattering
-#"
-#" @param m numeric. Extinction coefficient for direct (solar or observer) flux
-#" @param t numeric. Leaf Area Index
-#" @return jout numeric.
-#" @export
-jfunc4 <- function(m, t) {
-
-  del <- m * t
-  out <- 0 * del
-  out[del > 1e-3] <- (1 - exp(-del)) / (m * (1 + exp(-del)))
-  out[del <= 1e-3] <- 0.5 * t * (1. - del * del / 12.)
-  return(out)
-}
-
-
-#" compute volume scattering functions and interception coefficients
-#" for given solar zenith,  viewing zenith,  azimuth and leaf inclination angle.
-#"
-#" @param tts numeric. solar zenith
-#" @param tto numeric. viewing zenith
-#" @param psi numeric. azimuth
-#" @param ttl numeric. leaf inclination angle
-#" @return res list. includes chi_s,  chi_o,  frho,  ftau
-#" @export
-volscatt  <- function(tts, tto, psi, ttl) {
-  #********************************************************************************
-  #*	chi_s	= interception functions
-  #*	chi_o	= interception functions
-  #*	frho	= function to be multiplied by leaf reflectance rho
-  #*	ftau	= functions to be multiplied by leaf transmittance tau
-  #********************************************************************************
-  #	Wout Verhoef,  april 2001,  for CROMA
-
-  rd <- pi / 180
-  costs <- cos(rd * tts)
-  costo <- cos(rd * tto)
-  sints <- sin(rd * tts)
-  sinto <- sin(rd * tto)
-  cospsi <- cos(rd * psi)
-  psir <- rd * psi
-  costl <- cos(rd * ttl)
-  sintl <- sin(rd * ttl)
-  cs <- costl * costs
-  co <- costl * costo
-  ss <- sintl * sints
-  so <- sintl * sinto
-
-  #c ..............................................................................
-  #c     betas -bts- and betao -bto- computation
-  #c     Transition angles (beta) for solar (betas) and view (betao) directions
-  #c     if thetav + thetal > pi / 2,  bottom side of the leaves is observed for leaf azimut
-  #c     interval betao + phi<leaf azimut<2pi-betao + phi.
-  #c     if thetav + thetal<pi / 2,  top side of the leaves is always observed,  betao=pi
-  #c     same consideration for solar direction to compute betas
-  #c ..............................................................................
-
-  cosbts <- 5
-  if (abs(ss) > 1e-6) {
-    cosbts <- -cs / ss
-  }
-  cosbto <- 5
-  if (abs(so) > 1e-6) {
-    cosbto <- -co / so
-  }
-
-  if (abs(cosbts) < 1) {
-    bts <- acos(cosbts)
-    ds <- ss
-  } else {
-    bts <- pi
-    ds <- cs
-  }
-  chi_s <- 2. / pi * ((bts - pi * .5) * cs + sin(bts) * ss)
-  if (abs(cosbto) < 1) {
-    bto <- acos(cosbto)
-    doo <- so
-  } else if (tto < 90) {
-    bto <- pi
-    doo <- co
-  } else {
-    bto <- 0
-    doo <- -co
-  }
-  chi_o <- 2. / pi * ((bto - pi * .5) * co + sin(bto) * so)
-
-  #c ..............................................................................
-  #c   computation of auxiliary azimut angles bt1,  bt2,  bt3 used
-  #c   for the computation of the bidirectional scattering coefficient w
-  #c .............................................................................
-
-  btran1 <- abs(bts - bto)
-  btran2 <- pi - abs(bts + bto - pi)
-
-  if (psir <= btran1) {
-    bt1 <- psir
-    bt2 <- btran1
-    bt3 <- btran2
-  } else {
-    bt1 <- btran1
-    if (psir <= btran2) {
-      bt2 <- psir
-      bt3 <- btran2
-    } else {
-      bt2 <- btran2
-      bt3 <- psir
-    }
-  }
-  t1 <- 2. * cs * co + ss * so * cospsi
-  t2 <- 0
-  if (bt2 > 0) {
-    t2 <- sin(bt2) * (2. * ds * doo + ss * so * cos(bt1) * cos(bt3))
-  }
-
-  denom <- 2. * pi * pi
-  frho <- ((pi - bt2) * t1 + t2) / denom
-  ftau <- (-bt2 * t1 + t2) / denom
-
-  if (frho < 0) {
-    frho <- 0
-  }
-  if (ftau < 0) {
-    ftau <- 0
-  }
-  res <- list("chi_s" = chi_s, "chi_o" = chi_o, "frho" = frho, "ftau" = ftau)
-  return(res)
-}