initial commit

.gitignore

+.Rproj.user
+.Rbuildignore
+.Rhistory

DESCRIPTION

+Package: DiversityMappR
+Title: DiversityMappR to map biodiversity from remotely sensed data
+Version: 0.0.0.9000
+Authors@R: 
+    person(given = "Jean-Baptiste",
+           family = "Feret",
+           role = c("aut", "cre"),
+           email = "jb.feret@teledetection.fr")
+Description: this packages allows processing image data based on the method described in the following publication:
+    Féret, J.-B., Asner, G.P., 2014. Mapping tropical forest canopy diversity using high-fidelity imaging spectroscopy. Ecol. Appl. 24, 1289–1296. https://doi.org/10.1890/13-1824.1
+    It expects an image file as input, with a specific data format.
+    ENVI HDR image with BIL interleave required
+License: What license it uses
+Encoding: UTF-8
+LazyData: true
+Imports: 
+    R.utils,
+    matlab,
+    snow,
+    future.apply,
+    matrixStats,
+    fields,
+    ecodist,
+    zip,
+    vegan,
+    dissUtils,
+    labdsv,
+    raster,
+    rgdal
+RoxygenNote: 6.1.1

DiversityMappR.Rproj

+Version: 1.0
+
+RestoreWorkspace: No
+SaveWorkspace: No
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX
+
+AutoAppendNewline: Yes
+StripTrailingWhitespace: Yes
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source
+PackageRoxygenize: rd,collate,namespace

NAMESPACE

+# Generated by roxygen2: do not edit by hand
+
+S3method(split,line)
+export(Check.Data.Format)
+export(Convert.Raster2BIL)
+export(Map.Alpha.Diversity)
+export(Map.Beta.Diversity)
+export(Map.Spectral.Species)
+export(Perform.Radiometric.Filtering)

R/Lib_CheckConvertData.R

+# ==============================================================================
+# DiversityMappR
+# Lib_CheckConvertData.R
+# ==============================================================================
+# PROGRAMMERS:
+# Jean-Baptiste FERET <jb.feret@irstea.fr>
+# Copyright 2018/07 Jean-Baptiste FERET
+# ==============================================================================
+# This library checks if a raster in format expected for diversity mapping
+# if not, can convert the raster into BIL format expected for diversity mapping
+# ==============================================================================
+
+#' converts a raster into BIL format as expected by DivMapping codes
+#'
+#' @param Raster.Path full path for the raster to be converted
+#' @param Sensor name of the sensor
+#' @param Convert.Integer should data be converted into integer ?
+#' @param Multiplying.Factor multiplying factor (eg convert real reflectance values between 0 and 1 into integer between 0 and 10000)
+#' @param Output.Directory
+#' @param Multiplying.Factor.Last multiplying factor for last band
+#'
+#' @return Output.Path path for the image converted into ENVI BIL format
+#' @export
+Convert.Raster2BIL = function(Raster.Path,Sensor = 'unknown',Output.Directory = FALSE,Convert.Integer = TRUE,Multiplying.Factor = 1,Multiplying.Factor.Last = 1){
+
+  # get directory and file name of original image
+  Input.File  = basename(Raster.Path)
+  Input.Dir   = dirname(Raster.Path)
+  # define path where data will be stored
+  if (Output.Directory == FALSE){
+    Output.Path = paste(Input.Dir,'/Converted_DivMAP/',file_path_sans_ext(Input.File),sep='')
+  } else {
+    dir.create(Output.Directory, showWarnings = FALSE,recursive = TRUE)
+    Output.Path = paste(Output.Directory,'/',file_path_sans_ext(Input.File),sep='')
+  }
+  message('The converted file will be written in the following location:')
+  print(Output.Path)
+  Output.Dir   = dirname(Output.Path)
+  dir.create(Output.Dir, showWarnings = FALSE,recursive = TRUE)
+
+  # apply multiplying factors
+  message('reading initial file')
+  Output.Img          = Multiplying.Factor*brick(Raster.Path)
+  Last.Band.Name      = Output.Img@data@names[length(Output.Img@data@names)]
+  Output.Img[[Last.Band.Name]] = Multiplying.Factor.Last*Output.Img[[Last.Band.Name]]
+
+  # convert into integer
+  if (Convert.Integer==TRUE){
+    Output.Img = round(Output.Img)
+  }
+
+  # write raster
+  message('writing converted file')
+  if (Convert.Integer ==TRUE){
+    r = writeRaster(Output.Img, filename = Output.Path, format="EHdr", overwrite=TRUE,datatype='INT2S')
+  } else {
+    r = writeRaster(Output.Img, filename = Output.Path, format="EHdr", overwrite=TRUE)
+  }
+  hdr(r, format="ENVI")
+
+  # remove unnecessary files
+  File2Remove   = paste(Output.Path,'.aux.xml',sep='')
+  File2Remove2  = paste(file_path_sans_ext(Output.Path),'.aux.xml',sep='')
+  file.remove(File2Remove)
+
+  File2Remove   = paste(Output.Path,'.prj',sep='')
+  File2Remove2  = paste(file_path_sans_ext(Output.Path),'.prj',sep='')
+  if (file.exists(File2Remove)){
+    file.remove(File2Remove)
+  } else if  (file.exists(File2Remove2)){
+    file.remove(File2Remove2)
+  }
+
+  File2Remove = paste(Output.Path,'.sta',sep='')
+  File2Remove = paste(file_path_sans_ext(Output.Path),'.sta',sep='')
+  if (file.exists(File2Remove)){
+    file.remove(File2Remove)
+  } else if (file.exists(File2Remove2)){
+    file.remove(File2Remove2)
+  }
+
+  File2Remove = paste(Output.Path,'.stx',sep='')
+  File2Remove2 = paste(file_path_sans_ext(Output.Path),'.stx',sep='')
+  if (file.exists(File2Remove)){
+    file.remove(File2Remove)
+  } else if (file.exists(File2Remove2)){
+    file.remove(File2Remove2)
+  }
+
+  File2Rename = paste(file_path_sans_ext(Output.Path),'.hdr',sep='')
+  File2Rename2= paste(Output.Path,'.hdr',sep='')
+  if (file.exists(File2Rename)){
+    file.rename(from = File2Rename,to = File2Rename2)
+  }
+
+  # change dot into underscore
+  Output.Path.US = gsub(Output.Path,pattern = '[.]',replacement = '_')
+  if (!Output.Path.US==Output.Path){
+    file.rename(from = Output.Path,to = Output.Path.US)
+  }
+
+  Output.Path.US.HDR = paste(gsub(Output.Path,pattern = '[.]',replacement = '_'),'.hdr',sep='')
+  if (!Output.Path.US.HDR==paste(Output.Path,'.hdr',sep='')){
+    file.rename(from = paste(Output.Path,'.hdr',sep=''),to = Output.Path.US.HDR)
+    file.rename(from = Output.Path.US.HDR,to = Output.Path.US.HDR)
+  }
+
+
+  if (!Sensor=='unknown'){
+    # ! HARD-CODED
+    # HDR.Temp.Path = paste('DiversityMappR/data-raw/HDR/',Sensor,'.hdr',sep='')
+    HDR.Temp.Path = paste('D:/PUBLICATIONS/2019d_DiversityMappR/01_CODES/DiversityMappR/data-raw/HDR/',Sensor,'.hdr',sep='')
+    if (file.exists(HDR.Temp.Path)){
+      message('reading header template corresponding to the sensor located here:')
+      print(HDR.Temp.Path)
+      # get raster template corresponding to the sensor
+      HDR.Template  = read.ENVI.header(HDR.Temp.Path)
+      # get info to update hdr file
+      # read hdr
+      HDR.input     = read.ENVI.header(Get.HDR.Name(Output.Path))
+      if (!is.null(HDR.Template$wavelength)){
+        HDR.input$wavelength    = HDR.Template$wavelength
+      }
+      if (!is.null(HDR.Template$`sensor type`)){
+        HDR.input$`sensor type` = HDR.Template$`sensor type`
+      }
+      if (!is.null(HDR.Template$`band names`)){
+        HDR.input$`band names`  = HDR.Template$`band names`
+      }
+      if (!is.null(HDR.Template$`wavelength units`)){
+        HDR.input$`wavelength units`  = HDR.Template$`wavelength units`
+      }
+      # define visual stretch in the VIS domain
+      HDR.input$`default stretch` = '0 1000 linear'
+      # write corresponding hdr file
+      write.ENVI.header(HDR.input, Get.HDR.Name(Output.Path))
+
+    } else if (!file.exists(HDR.Temp.Path)){
+      message('Header template corresponding to the sensor expected to be found here')
+      print(HDR.Temp.Path)
+      message('please provide this header template in order to write info in HDR file')
+      print(Get.HDR.Name(Output.Path))
+      message('or manually add wavelength location in HDR file, if relevant')
+    }
+  } else if (Sensor=='unknown'){
+    message('please make sure that the follozing header file contains information required')
+    print(Get.HDR.Name(Output.Path))
+    message('or manually add wavelength location in HDR file, if relevant')
+  }
+  return(Output.Path)
+}
+
+#' Checks if the data to be processed has the format type expected
+#'
+#' @param Raster.Path full path for the raster to be converted
+#' @param Mask is the raster a mask?
+#'
+#' @return Output.Path path for the image converted into ENVI BIL format
+#' @export
+Check.Data.Format = function(Raster.Path,Mask = FALSE){
+
+  HDR.Path = Get.HDR.Name(Raster.Path)
+  # check if the hdr file exists
+  if (file.exists(HDR.Path)){
+    HDR     = read.ENVI.header(HDR.Path)
+    if (Mask == FALSE & (!HDR$interleave=='bil') & (!HDR$interleave=='BIL')){
+      message('')
+      message('*********************************************************')
+      message('The image format may not compatible with the processing chain')
+      message('Image format expected:')
+      message('ENVI hdr file with band interleaved by line (BIL) file format')
+      message('')
+      message('Current Image format')
+      print(HDR$interleave)
+      message('Please run the function named ')
+      print('Convert.Raster2BIL')
+      message('in order to convert your raster data')
+      message('or use appropriate software')
+      message('*********************************************************')
+      message('')
+      stop()
+    } else if (Mask == FALSE & ((HDR$interleave=='bil') | (HDR$interleave=='BIL'))){
+      if (HDR$`wavelength units`=='Unknown'){
+        message('')
+        message('*********************************************************')
+        message('Please make sure the wavelengths are in nanometers')
+        message('if not, stop processing and convert wavelengths in nanometers in HDR file')
+        message('*********************************************************')
+        message('')
+      }
+      if ((!HDR$`wavelength units`=='Nanometers') & (!HDR$`wavelength units`=='nanometers')){
+        message('')
+        message('*********************************************************')
+        message('Please make sure the wavelengths are in nanometers')
+        message('if not, stop processing and convert wavelengths in nanometers in HDR file')
+        message('*********************************************************')
+        message('')
+      }
+      if (HDR$`wavelength units`=='micrometers'){
+        message('')
+        message('*********************************************************')
+        message('Please convert wavelengths in nanometers in HDR file')
+        message('*********************************************************')
+        message('')
+        stop()
+      }
+      if ((HDR$`wavelength units`=='nanometers') | (HDR$`wavelength units`=='Nanometers')){
+        message('')
+        message('*********************************************************')
+        message('       All information seem OK for image processing      ')
+        message('*********************************************************')
+        message('')
+      }
+    }
+  } else {
+    message('')
+    message('*********************************************************')
+    message('The following HDR file was expected, but could not be found:')
+    print(HDR.Path)
+    message('The image format may not compatible with the processing chain')
+    message('Image format expected:')
+    message('ENVI hdr file with band interleaved by line (BIL) file format')
+    message('')
+    message('Please run the function named ')
+    print('Convert.Raster2BIL')
+    message('in order to convert your raster data')
+    message('*********************************************************')
+    message('')
+    stop()
+  }
+  return()
+}

R/Lib_ContinuumRemoval.R

+# ===============================================================================
+# DiversityMappR
+# Lib_ContinuumRemoval.R
+# ===============================================================================
+# PROGRAMMERS:
+# Jean-Baptiste FERET <jb.feret@irstea.fr>
+# Copyright 2018/07 Jean-Baptiste FERET
+# ===============================================================================
+# This Library is dedicated to the computation of the continuum removal
+# ===============================================================================
+
+#' prepares data to run multithreaded continuum removal
+#'
+#' @param Spectral.Data initial data matrix (nb samples x nb bands)
+#' @param nbCPU
+#' @param Spectral information about spectral bands
+#'
+#' @return samples from image and updated number of pixels to sampel if necessary
+Apply.Continuum.Removal=function(Spectral.Data,Spectral,nbCPU=1){
+
+  if (!length(Spectral$WaterVapor)==0){
+    Spectral.Data=Spectral.Data[,-Spectral$WaterVapor]
+  }
+
+  # split data to perform continuum removal on into reasonable amount of data
+  nb.Values     = size(Spectral.Data)[1]*size(Spectral.Data)[2]
+  if (nb.Values>0){
+    # corresponds to ~ 40 Mb data, but CR tends to requires ~ 10 times memory
+    # avoids memory crash
+    Max.nb.Values   = 2e6
+    nb.CR           = ceiling(nb.Values/Max.nb.Values)
+    Spectral.Data   = splitRows(Spectral.Data,nb.CR)
+    # perform multithread continuum removal
+    plan(multiprocess, workers = nbCPU) ## Parallelize using four cores
+    Schedule.Per.Thread = ceiling(nb.CR/nbCPU)
+    Spectral.Data.tmp=future_lapply(Spectral.Data,FUN = ContinuumRemoval,Spectral.Bands=Spectral$Wavelength,future.scheduling = Schedule.Per.Thread)
+    plan(sequential)
+    Spectral.Data = do.call("rbind",Spectral.Data.tmp)
+    rm(Spectral.Data.tmp)
+  } else {
+    # edit 31-jan-2018
+    # resize to delete first and last band as in continuum removal
+    Spectral.Data=Spectral.Data[,-c(1,2)]
+  }
+  gc()
+  return(Spectral.Data)
+}
+
+#' Computes continuum removal for matrix shaped data: more efficient than
+#' processing individual spectra
+#' the convex hull is based on the computation of the derivative between R at a
+#' given spectral band and R at the following bands
+#'
+#' @param Minit initial data matrix (nb samples x nb bands)
+#' @param Spectral.Bands information about spectral bands
+#'
+#' @return samples from image and updated number of pixels to sampel if necessary
+ContinuumRemoval=function(Minit,Spectral.Bands){
+
+  # Filter and prepare data prior to continuum removal
+  CR.data         = Filter.Prior.CR(Minit,Spectral.Bands)
+  Minit           = CR.data$Minit
+  nb.Bands        = size(Minit)[2]
+  CR.data$Minit   = c()
+  Spectral.Bands  = CR.data$Spectral.Bands
+  nb.Samples      = CR.data$nb.Samples
+  # if samples to be considered
+  if (nb.Samples>0){
+    # initialization:
+    # spectral band corresponding to each element of the data matrix
+    Lambda          = repmat(matrix(Spectral.Bands,nrow=1),nb.Samples,1)
+    # prepare matrices used to check evolution of the CR process
+    # - elements still not processed through continuum removal: initialization to 1
+    Still.Need.CR   = matrix(1,nrow=nb.Samples,ncol=nb.Bands)
+    # - value of the convex hull: initially set to 0
+    Convex.Hull     = matrix(0,nrow=nb.Samples,ncol=nb.Bands)
+    # - reflectance value for latest interception with convex hull:
+    # initialization to value of the first reflectance measurement
+    Intercept.Hull  = repmat(matrix(Minit[,1],ncol=1),1,nb.Bands)
+    # - spectral band of latest interception
+    Latest.Intercept= repmat(matrix(Spectral.Bands[1],ncol=1),nb.Samples,nb.Bands)
+    # number of spectral bands found as intercept
+    nb.Intercept    = 0
+    # continues until arbitrary stopping criterion:
+    # stops when reach last spectral band (all values before last = 0)
+    while (max(Still.Need.CR[,1:(nb.Bands-2)])==1 & (nb.Intercept <=(nb.Bands/2))){
+      nb.Intercept      = nb.Intercept+1
+      # Mstep give the position of the values to be updated
+      Update.Data       = matrix(1,nrow=nb.Samples,ncol=nb.Bands)
+      Update.Data[,nb.Bands]=0
+      # initial step: first column set to 0; following steps: all bands below
+      # max of the convex hull are set to 0
+      Update.Data[which((Lambda-Latest.Intercept)<0)]=0
+      # compute slope for each coordinate
+      Slope             = (Minit-Intercept.Hull)/(Lambda-Latest.Intercept)*Still.Need.CR
+      # set current spectral band and previous bands to -9999
+      if (!length(which(Still.Need.CR==0))==0){
+        Slope[which(Still.Need.CR==0)]=-9999
+      }
+      if (!length(which(is.na(Slope)))==0){
+        Slope[which(is.na(Slope))]=-9999
+      }
+      # get max index for each row and convert into linear index
+      Index.Max.Slope   = RowToLinear(max.col(Slope,ties.method="last"),nb.Samples,nb.Bands)
+      # !!!! OPTIM: replace repmat with column operation
+      # update coordinates of latest intercept
+      Latest.Intercept  = repmat(matrix(Lambda[Index.Max.Slope],ncol=1),1,nb.Bands)
+      # update latest intercept
+      Intercept.Hull    = repmat(matrix(Minit[Index.Max.Slope],ncol=1),1,nb.Bands)
+      # values corresponding to the domain between the two continuum maxima
+      Update.Data[which((Lambda-Latest.Intercept)>=0 | Latest.Intercept==Spectral.Bands[nb.Bands])]=0
+      # values to eliminate for the next analysis: all spectral bands before latest intercept
+      Still.Need.CR[which((Lambda-Latest.Intercept)<0)] = 0
+      # the max slope is known, as well as the coordinates of the beginning and ending
+      # a matrix now has to be built
+      Convex.Hull=Convex.Hull+Update.Data*(Intercept.Hull+sweep((Lambda-Latest.Intercept),1,Slope[Index.Max.Slope],"*"))
+    }
+    CR_Results0               = Minit[,2:(nb.Bands-2)]/Convex.Hull[,2:(nb.Bands-2)]
+    CR_Results                = matrix(0,ncol=(nb.Bands-3),nrow=nb.Samples)
+    CR_Results[CR.data$Samples.To.Keep,] = CR_Results0
+  } else {
+    CR_Results        = matrix(0,ncol=(nb.Bands-3),nrow=nb.Samples)
+  }
+  list          = ls()
+  rm(list=list[-which(list=="CR_Results")])
+  rm(list)
+  gc()
+  return(CR_Results)
+}
+
+#' Filter data prior to continuum removal:
+#' - values are expected to be real reflectance values between 0 and 10000
+#' - negative values may occur, so a +100 value is applied to avoid negative
+#' - possibly remaining negative values are set to 0
+#' - constant spectra are eliminated
+#'
+#' @param Spectral.Bands
+#' @param Minit initial data matrix, n rows = n samples, p cols = p spectral bands
+#'
+#' @return updated Minit
+Filter.Prior.CR = function(Minit,Spectral.Bands){
+
+  # number of samples to be processed
+  nb.Samples  = nrow(Minit)
+  # make sure there is no negative values
+  Minit       = Minit+100.0
+  Minit[which(Minit<0)]=0
+  # eliminate invariant spectra
+  SD          = rowSds(Minit)
+  elim        = which(SD==0 | is.na(SD))
+  keep        = which(!SD==0 & !is.na(SD))
+  nbi0        = nrow(Minit)
+  nbj         = ncol(Minit)
+  if (nbj==1 & nbi0>1){
+    Minit     = matrix(Minit,nrow=1)
+    nbi0      = 1
+    nbj       = ncol(Minit)
+  }
+  Minit       = Minit[keep,]
+  if (length(keep)==1){
+    Minit     = matrix(Minit,nrow=1)
+  }
+  # add negative values to the last column and update spectral bands
+  Minit       = cbind(Minit, matrix(-9999,ncol=1,nrow=nb.Samples))
+  nb.Bands    = ncol(Minit)
+  Spectral.Bands[nb.Bands]=Spectral.Bands[nb.Bands-1]+10
+  my_list <- list("Minit"=Minit,"Spectral.Bands"=Spectral.Bands,"nb.Samples"=nb.Samples,"Samples.To.Keep"=keep)
+  return(my_list)
+}
+
+#'
+#' @param MM
+#' @param nbi
+#' @param nbj
+#'
+#' @return
+RowToLinear=function(MM,nbi,nbj){
+
+  adj=seq(1:nbi)
+  MaxCont=((MM-1)*(nbi))+adj
+  return(MaxCont)
+}
+
+#' R equivalent of repmat (matlab)
+#'
+#' @param X initial matrix
+#' @param m nb of replications in row dimension
+#' @param n nb of replications in column dimension
+#'
+#' @return
+repmat = function(X,m,n){
+  mx = dim(X)[1]
+  nx = dim(X)[2]
+  matrix(t(matrix(X,mx,nx*n)),mx*m,nx*n,byrow=T)
+}

R/Lib_FilterData.R

+# ==============================================================================
+# DiversityMappR
+# Lib_FilterData.R
+# ==============================================================================
+# PROGRAMMERS:
+# Jean-Baptiste FERET <jb.feret@irstea.fr>
+# Copyright 2018/07 Jean-Baptiste FERET
+# ==============================================================================
+# This library contains functions to filter raster based on radiometric criteria
+# ==============================================================================
+
+#' Performs radiometric filtering based on three criteria: NDVI, NIR reflectance, Blue reflectance
+#'
+#' @param Image.Path Path of the image to be processed
+#' @param Mask.Path Path of the mask corresponding to the image
+#' @param Output.Dir output directory
+#' @param TypePCA Type of PCA (PCA, SPCA, NLPCA...)
+#' @param NDVI.Thresh NDVI threshold applied to produce a mask (select pixels with NDVI>NDVI.Thresh)
+#' @param Blue.Thresh Blue threshold applied to produce a mask (select pixels with Blue refl < Blue.Thresh --> filter clouds) refl expected between 0 and 10000
+#' @param NIR.Thresh NIR threshold applied to produce a mask (select pixels with NIR refl < NIR.Thresh) refl expected between 0 and 10000
+#' @param Mask.Path
+#'
+#' @return ImPathShade = updated shademask file
+#' @export
+Perform.Radiometric.Filtering = function(Image.Path,Mask.Path,Output.Dir,TypePCA='SPCA',NDVI.Thresh = 0.5, Blue.Thresh = 500, NIR.Thresh = 1500,Blue = 480, Red = 700, NIR = 835){
+  # define full output directory
+  Output.Dir.Full   = Define.Output.Dir(Output.Dir,Image.Path,TypePCA)
+  # define dimensions of the image
+  ImPathHDR         = Get.HDR.Name(Image.Path)
+  HDR               = read.ENVI.header(ImPathHDR)
+  Image.Format      = ENVI.Type2Bytes(HDR)
+  ipix              = as.double(HDR$lines)
+  jpix              = as.double(HDR$samples)
+  Nb.Pixels         = ipix*jpix
+  lenTot            = Nb.Pixels*as.double(HDR$bands)
+  ImSizeGb          = (lenTot*Image.Format$Bytes)/(1024^3)
+
+  # Create / Update shade mask if optical data
+  if (Mask.Path == FALSE | Mask.Path == ''){
+    print('Create mask based on NDVI, NIR and Blue threshold')
+  } else {
+    print('Update mask based on NDVI, NIR and Blue threshold')
+  }
+  Shade.Update  = paste(Output.Dir.Full,'ShadeMask_Update',sep='')
+  Mask.Path   = Create.Mask.Optical(Image.Path,Mask.Path,Shade.Update,NDVI.Thresh,Blue.Thresh,NIR.Thresh,Blue,Red,NIR)
+  return(Mask.Path)
+}
+
+#' create a mask based on NDVI, Green reflectance and NIR reflectance
+#' NDVI (min) threshold eliminates non vegetated pixels
+#' Blue (max) threshold eliminates Clouds
+#' NIR (min) threshold eliminates shadows
+#' ! only valid if Optical data!!
+#'
+#' @param ImPath full path of a raster file
+#' @param ImPathShade full path of the raster mask corresponding to the raster file
+#' @param ImPathShade.Update wavelength (nm) of the spectral bands to be found
+#' @param NDVI.Thresh NDVI threshold applied to produce a mask (select pixels with NDVI>NDVI.Thresh)
+#' @param Blue.Thresh Blue threshold applied to produce a mask (select pixels with Blue refl < Blue.Thresh --> filter clouds) refl expected between 0 and 10000
+#' @param NIR.Thresh NIR threshold applied to produce a mask (select pixels with NIR refl < NIR.Thresh) refl expected between 0 and 10000
+#'
+#' @return ImPathShade path for the updated shademask produced
+Create.Mask.Optical = function(ImPath,ImPathShade,ImPathShade.Update,NDVI.Thresh,Blue.Thresh,NIR.Thresh, Blue = 480, Red = 700, NIR = 835){
+  # define wavelength corresponding to the spectral domains Blue, Red and NIR
+  Spectral.Bands= c(Blue,Red,NIR)
+  ImPathHDR     = Get.HDR.Name(ImPath)
+  Header        = read.ENVI.header(ImPathHDR)
+  # get image bands correponding to spectral bands of interest
+  Image.Bands   = Get.Image.Bands(Spectral.Bands,Header$wavelength)
+  # read band data from image
+  Image.Subset  = Read.Image.Bands(ImPath,Header,Image.Bands$ImBand)
+  # create mask
+  # check if spectral bands required for NDVI exist
+  if (Image.Bands$Distance2WL[2]<25 & Image.Bands$Distance2WL[3]<25){
+    NDVI          = ((Image.Subset[,,3])-(Image.Subset[,,2]))/((Image.Subset[,,3])+(Image.Subset[,,2]))
+  } else {
+    NDVI          = matrix(1,nrow = Header$lines,ncol = Header$samples)
+    message('Could not find the spectral bands required to compute NDVI')
+  }
+  if (Image.Bands$Distance2WL[1]>25){
+    Image.Subset[,,1] = Blue.Thresh+0*Image.Subset[,,1]
+    message('Could not find a spectral band in the blue domain: will not perform filtering based on blue reflectance')
+  }
+  if (Image.Bands$Distance2WL[3]>50){
+    Image.Subset[,,3] = NIR.Thresh+0*Image.Subset[,,3]
+    message('Could not find a spectral band in the NIR domain: will not perform filtering based on NIR reflectance')
+  }
+  Mask          = matrix(0,nrow=Header$lines,ncol=Header$samples)
+  SelPixels     = which(NDVI>NDVI.Thresh & Image.Subset[,,1]<Blue.Thresh & Image.Subset[,,3]>NIR.Thresh)
+  Mask[SelPixels] = 1
+  # update initial shade mask
+  ImPathShade   = Update.Shademask(ImPathShade,Header,Mask,ImPathShade.Update)
+  return(ImPathShade)
+}