R/Lib_PROSPECT.R

Defines functions colour_to_ansi Complete_Input_PROSPECT download_LeafDB PROSPECT_LUT FitSpectralData define_Input_PROSPECT check_version_prospect calctav PROSPECT

Documented in calctav check_version_prospect colour_to_ansi Complete_Input_PROSPECT define_Input_PROSPECT download_LeafDB FitSpectralData PROSPECT_LUT

#' #############################################################################
#' prospect
#' Lib_PROSPECT.R
#' #############################################################################
#' PROGRAMMERS:
#' Jean-Baptiste FERET jb.feret@teledetection.fr
#' Florian de Boissieu florian.deboissieu@inrae.fr
#' Copyright 2020/04 Jean-Baptiste FERET
#' =============================================================================
#' This Library includes functions dedicated to PROSPECT simulation
#' =============================================================================

#' core function running PROSPECT
#' This function allows simulations using PROSPECT-D or PROSPECT-PRO depending
#' on the parameterization.
#  This code includes numerical optimizations proosed in the FLUSPECT code
#  Authors: Wout Verhoef, Christiaan van der Tol (tol@itc.nl), Joris Timmermans,
#  Date: 2007
#  Update from PROSPECT to FLUSPECT: January 2011 (CvdT)
#'
#' @param SpecPROSPECT list. Includes spectral constants derived from
#' SpecPROSPECT_FullRange: refractive index, specific absorption coefficients
#' and corresponding spectral bands
#' @param Input_PROSPECT list. Includes all prospect input parameters
#' @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. Anthocyanin 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. NonProt Carbon-based constituent content (g.cm-2)
#' @param alpha numeric. Solid angle for incident light at surface of leaf
#' @param check boolean. set to TRUE to check input data format
#'
#' @return leaf directional-hemispherical reflectance and transmittance
#' @importFrom expint expint
#' @export
#'
PROSPECT <- function(SpecPROSPECT = NULL, Input_PROSPECT = NULL,
                     N = 1.5, CHL = 40.0, CAR = 8.0, ANT = 0.0, BROWN = 0.0,
                     EWT = 0.01, LMA = NULL, PROT = 0, CBC = 0, alpha = 40.0,
                     check = TRUE) {

  # define PROSPECT input in a dataframe
  Input_PROSPECT <- define_Input_PROSPECT(Input_PROSPECT, CHL, CAR,
                                          ANT, BROWN, EWT, LMA, PROT,
                                          CBC, N, alpha)
  # if (check) Input_PROSPECT <- define_Input_PROSPECT(Input_PROSPECT, CHL, CAR,
  #                                                    ANT, BROWN, EWT, LMA, PROT,
  #                                                    CBC, N, alpha)
  # default: simulates leaf optics using full spectral range
  if (is.null(SpecPROSPECT)) SpecPROSPECT <- prospect::SpecPROSPECT_FullRange
  # compute total absorption corresponding to each homogeneous layer
  Kall <- (Input_PROSPECT$CHL * SpecPROSPECT$SAC_CHL +
             Input_PROSPECT$CAR * SpecPROSPECT$SAC_CAR +
             Input_PROSPECT$ANT * SpecPROSPECT$SAC_ANT +
             Input_PROSPECT$BROWN * SpecPROSPECT$SAC_BROWN +
             Input_PROSPECT$EWT * SpecPROSPECT$SAC_EWT +
             Input_PROSPECT$LMA * SpecPROSPECT$SAC_LMA +
             Input_PROSPECT$PROT * SpecPROSPECT$SAC_PROT +
             Input_PROSPECT$CBC * SpecPROSPECT$SAC_CBC) / Input_PROSPECT$N

  # Non-conservative scattering (normal case) when Kall > 0
  j <- which(Kall <= 0)
  t1 <- (1 - Kall) * exp(-Kall)
  t2 <- (Kall * Kall) * expint(Kall)
  tau <- t1 + t2
  if (length(j)>0) tau[j] <- 1
  # tau <- matrix(1, ncol = 1, nrow = length(t1))
  # tau[j] <- t1[j] + t2[j]

  # ***********************************************************************
  # reflectance and transmittance of one layer
  # ***********************************************************************
  # Allen W.A., Gausman H.W., Richardson A.J., Thomas J.R. (1969),
  # Interaction of isotropic ligth with a compact plant leaf, J. Opt.
  # Soc. Am., 59(10):1376-1379.
  # ***********************************************************************
  # reflectivity and transmissivity at the interface
  # ***********************************************************************
  if (Input_PROSPECT$alpha == 40) {
    talf <- SpecPROSPECT$calctav_40
  } else {
    talf <- calctav(Input_PROSPECT$alpha, SpecPROSPECT$nrefrac)
  }
  ralf <- 1 - talf
  t12 <- SpecPROSPECT$calctav_90
  r12 <- 1 - t12
  t21 <- t12 / (SpecPROSPECT$nrefrac**2)
  r21 <- 1 - t21

  # top surface side
  denom <- 1 - (r21 * r21 * (tau**2))
  Ta <- (talf * tau * t21) / denom
  Ra <- ralf + (r21 * tau * Ta)
  # bottom surface side
  t <- t12 * tau * t21 / denom
  r <- r12 + (r21 * tau * t)

  # ***********************************************************************
  # reflectance and transmittance of N layers
  # Stokes equations to compute properties of next N-1 layers (N real)
  # Normal case
  # ***********************************************************************
  # Stokes G.G. (1862), On the intensity of the light reflected from
  # or transmitted through a pile of plates, Proc. Roy. Soc. Lond.,
  # 11:545-556.
  # ***********************************************************************
  D <- sqrt((1 + r + t) * (1 + r - t) * (1 - r + t) * (1 - r - t))
  rq <- r**2
  tq <- t**2
  a <- (1 + rq - tq + D) / (2 * r)
  b <- (1 - rq + tq + D) / (2 * t)

  bNm1 <- b**(Input_PROSPECT$N - 1)
  bN2 <- bNm1**2
  a2 <- a**2
  denom <- a2 * bN2 - 1
  Rsub <- a * (bN2 - 1) / denom
  Tsub <- bNm1 * (a2 - 1) / denom

  # 	Case of zero absorption
  j <- which(r + t >= 1)
  Tsub[j] <- t[j] / (t[j] + (1 - t[j]) * (Input_PROSPECT$N - 1))
  Rsub[j] <- 1 - Tsub[j]

  # leaf reflectance and transmittance : combine top layer with next N-1 layers
  denom <- 1 - Rsub * r
  tran <- Ta * Tsub / denom
  refl <- Ra + (Ta * Rsub * t) / denom
  return(data.frame('wvl' = SpecPROSPECT$lambda,
                    'Reflectance' = refl,
                    'Transmittance' = tran))
}

#' computation of transmissivity of a dielectric plane surface,
#' averaged over all directions of incidence and over all polarizations.
#'
#' @param alpha numeric. max incidence angle of solid angle of incident light
#' @param nr numeric. refractive index
#'
#' @return numeric. Transmissivity of a dielectric plane surface
#' @export
calctav <- function(alpha, nr) {
  # Stern F. (1964), Transmission of isotropic radiation across an
  # interface between two dielectrics, Appl. Opt., 3(1):111-113.
  # Allen W.A. (1973), Transmission of isotropic light across a
  # dielectric surface in two and three dimensions, J. Opt. Soc. Am.,
  # 63(6):664-666.
  # ***********************************************************************

  rd <- pi / 180
  n2 <- nr**2
  np <- n2 + 1
  nm <- n2 - 1
  a <- (nr + 1) * (nr + 1) / 2
  k <- -(n2 - 1) * (n2 - 1) / 4
  sa <- sin(alpha * rd)

  b2 <- (sa**2) - (np / 2)
  if (alpha == 90) {
    b1 <- 0 * b2
  } else {
    b1 <- sqrt((b2**2) + k)
  }
  b <- b1 - b2
  b3 <- b**3
  a3 <- a**3
  ts <- ((k**2) / (6 * b3) + (k / b) - b / 2) - ((k**2) / (6 * a3) + (k / a) - (a / 2))

  tp1 <- -2 * n2 * (b - a) / (np**2)
  tp2 <- -2 * n2 * np * log(b / a) / (nm**2)
  tp3 <- n2 * ((1 / b) - (1 / a)) / 2
  tp4 <- 16 * n2**2 * ((n2**2) + 1) * log(((2 * np * b) - (nm**2)) / ((2 * np * a) - (nm**2))) / ((np**3) * (nm**2))
  tp5 <- 16 * (n2**3) * (1 / ((2 * np * b) - (nm**2)) - (1 / (2 * np * a - (nm**2)))) / (np**3)
  tp <- tp1 + tp2 + tp3 + tp4 + tp5
  tav <- (ts + tp) / (2 * (sa**2))
  return(tav)
}

#' This function checks if the input parameters are defined as expected
#' to run either PROSPECT-D or PROSPECT-PRO
#' @param LMA numeric. content corresponding to LMA
#' @param PROT numeric. content corresponding to protein content
#' @param CBC numeric. content corresponding to carbon based constituents
#'
#' @return list. updated LMA, PROT and CBC
#' @export

check_version_prospect <- function(LMA, PROT, CBC){
  # PROSPECT-D as default value
  if (is.null(LMA) & PROT == 0 & CBC == 0) LMA <- 0.008
  # PROSPECT-PRO if PROT or CBC are not NULL
  if (is.null(LMA) & (PROT > 0 | CBC > 0)) LMA <- 0
  # if calling PROSPECT-PRO (protein content or CBC defined by user)
  # then set LMA to 0 in any case
  if (!LMA==0 & (PROT > 0 | CBC > 0)) {
    print_msg('version_PROSPECT')
    LMA <- 0
  }
  return(list('LMA' = LMA, 'PROT' = PROT, 'CBC' = CBC))
}


#' This function produces a data frame from all prospect input variables if not
#' defined already
#'
#'

#' @param Input_PROSPECT numeric. content corresponding to LMA
#' @param CHL numeric. Chlorophyll content (microg.cm-2)
#' @param CAR numeric. Carotenoid content (microg.cm-2)
#' @param ANT numeric. Anthocyanin 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. NonProt Carbon-based constituent content (g.cm-2)
#' @param N numeric. Leaf structure parameter
#' @param alpha numeric. Solid angle for incident light at surface of leaf
#'
#' @return list. updated LMA, PROT and CBC
#' @export

define_Input_PROSPECT <- function(Input_PROSPECT, CHL = NULL, CAR = NULL,
                                  ANT = NULL, BROWN = NULL, EWT = NULL,
                                  LMA = NULL, PROT = NULL, CBC = NULL,
                                  N = NULL, alpha = NULL){

  default_PROSPECT <- data.frame('CHL' = 40.0, 'CAR' = 8.0, 'ANT' = 0.0,
                                 'BROWN' = 0.0, 'EWT' = 0.01, 'LMA' = 0.0,
                                 'PROT'= 0.0, 'CBC' = 0.0, 'N' = 1.5,
                                 'alpha' = 40.0)
  if (is.null(Input_PROSPECT)){
    dm_val <- prospect::check_version_prospect(LMA, PROT, CBC)
    Input_PROSPECT <- data.frame('CHL' = CHL, 'CAR' = CAR, 'ANT' = ANT,
                                 'BROWN' = BROWN, 'EWT' = EWT, 'LMA' = dm_val$LMA,
                                 'PROT'= dm_val$PROT, 'CBC' = dm_val$CBC,
                                 'N' = N, 'alpha' = alpha)
  } else if (!is.null(Input_PROSPECT)){
    missing <- which(!names(default_PROSPECT)%in%names(Input_PROSPECT))
    if (length(missing)>0) Input_PROSPECT <- cbind(Input_PROSPECT, default_PROSPECT[missing])
    dm_val <- prospect::check_version_prospect(Input_PROSPECT$LMA,
                                               Input_PROSPECT$PROT,
                                               Input_PROSPECT$CBC)
    Input_PROSPECT$LMA <- dm_val$LMA
    Input_PROSPECT$PROT <- dm_val$PROT
    Input_PROSPECT$CBC <- dm_val$CBC
  }
  return(Input_PROSPECT)
}


#' This function adapts SpecPROSPECT accordingly to experimental data
#' or to a spectral domain defined by UserDomain
#' @param lambda numeric. Spectral bands corresponding to experimental data
#' @param SpecPROSPECT list. Includes optical constants: refractive index,
#' specific absorption coefficients and corresponding spectral bands
#' @param Refl numeric. Measured reflectance data
#' @param Tran numeric. Measured Transmittance data
#' @param UserDomain numeric. either Lower and upper bounds for domain of
#' interest (optional) or list of spectral bands of interest
#' @param UL_Bounds boolean. set to TRUE if UserDomain only includes lower and
#' upper band, set to FALSE if UserDomain is a list of spectral bands (in nm)
#'
#' @return list including spectral properties at the new resolution
#' @import dplyr
#' @export

FitSpectralData <- function(lambda, SpecPROSPECT = NULL,
                            Refl = NULL, Tran = NULL,
                            UserDomain = NULL, UL_Bounds = FALSE) {
  # default: simulates leaf optics using full spectral range
  if (is.null(SpecPROSPECT)) SpecPROSPECT <- prospect::SpecPROSPECT_FullRange
  # convert Refl and Tran into dataframe if needed
  if (inherits(Refl, what = c('numeric', 'matrix'))) Refl <- data.frame(Refl)
  if (inherits(Tran, what = c('numeric', 'matrix'))) Tran <- data.frame(Tran)
  # if (class(Refl)[1]%in%c('numeric', 'matrix')) Refl <- data.frame(Refl)
  # if (class(Tran)[1]%in%c('numeric', 'matrix')) Tran <- data.frame(Tran)
  # Adjust LOP: check common spectral domain between PROSPECT and leaf optics
  if (!is.null(Refl)) Refl <- Refl %>% filter(lambda%in%SpecPROSPECT$lambda)
  if (!is.null(Tran)) Tran <- Tran %>% filter(lambda%in%SpecPROSPECT$lambda)
  lambda <- lambda[lambda%in%SpecPROSPECT$lambda]
  # Adjust PROSPECT
  lb <- lambda
  SpecPROSPECT  <- SpecPROSPECT %>% filter(SpecPROSPECT$lambda%in%lb)
  # if UserDomain is defined
  if (is.null(UserDomain)) UserDomain <- lambda
  if (UL_Bounds==TRUE) UserDomain <- seq(min(UserDomain), max(UserDomain))
  if (!is.null(Refl)) Refl <- Refl %>% filter(lambda%in%UserDomain)
  if (!is.null(Tran)) Tran <- Tran %>% filter(lambda%in%UserDomain)
  lambda <- lambda[lambda%in%UserDomain]
  # Adjust PROSPECT
  SpecPROSPECT  <- SpecPROSPECT %>% filter(SpecPROSPECT$lambda%in%UserDomain)
  if (any(!UserDomain%in%lambda)){
    message('leaf optics out of range defined by UserDomain')
  }
  RT <- reshape_lop4inversion(Refl = Refl,
                              Tran = Tran,
                              SpecPROSPECT = SpecPROSPECT)
  return(list("SpecPROSPECT" = SpecPROSPECT, "lambda" = lambda,
              "Refl" = RT$Refl, "Tran" = RT$Tran, "nbSamples" = RT$nbSamples))
}

#' computation of a LUT of leaf optical properties using a set of
#' leaf chemical & structural parameters
#'
#' @param Input_PROSPECT dataframe. list of PROSPECT input parameters.
#' @param SpecPROSPECT list. spectral constants
#' refractive index, specific absorption coefficients & spectral bands
#'
#' @return list. LUT including leaf reflectance and transmittance
#' @export
PROSPECT_LUT <- function(Input_PROSPECT, SpecPROSPECT = NULL) {

  if (is.null(SpecPROSPECT)) SpecPROSPECT <- prospect::SpecPROSPECT_FullRange
  # expected PROSPECT input parameters
  ExpectedParms <- data.frame('CHL' = 0, 'CAR' = 0, 'ANT' = 0, 'BROWN' = 0,
                              'EWT' = 0, 'LMA' = 0, 'PROT' = 0, 'CBC' = 0,
                              'N' = 1.5, 'alpha' = 40)
  # parameters provided
  inOK <- names(Input_PROSPECT)
  # identify missing elements
  Parm2Add <- which(!names(ExpectedParms) %in% inOK)
  # check if all parameters are included.
  if (length(Parm2Add) > 0) print_msg(cause = 'Missing_Input',
                                      args = list('Input' = inOK,
                                                  'Expected' = ExpectedParms))

  # re-order missing elements the end of the list using default value
  Input_PROSPECT <- Complete_Input_PROSPECT(Input_PROSPECT = Input_PROSPECT,
                                            Parm2Add = Parm2Add,
                                            ExpectedParms = ExpectedParms)
  # print number of samples to be simulated
  nbSamples <- nrow(Input_PROSPECT)
  messageLUT <- paste('A LUT with', nbSamples, 'samples will be produced')
  cat(colour_to_ansi('green'), messageLUT, "\033[0m\n")

  # run PROSPECT for nbSamples
  run_list_PROSPECT <- function(Input_PROSPECT, SpecPROSPECT){
    LUT_tmp <- do.call(PROSPECT,
                       c(list(SpecPROSPECT = SpecPROSPECT),
                         Input_PROSPECT))
    return(LUT_tmp)
  }
  indiv_leaves <- split(Input_PROSPECT,
                        factor(seq(1:nbSamples)))
  LUT_tmp <- lapply(X = indiv_leaves,
                    FUN = run_list_PROSPECT,
                    SpecPROSPECT = SpecPROSPECT)
  LUT_Refl <- as.data.frame(lapply(LUT_tmp,'[[', 'Reflectance'))
  LUT_Tran <- as.data.frame(lapply(LUT_tmp,'[[', 'Transmittance'))
  names(LUT_Refl) <- names(LUT_Tran) <- paste0('sample_',seq(1,nbSamples))
  return(list('Reflectance' = LUT_Refl,
              'Transmittance' = LUT_Tran,
              'Input_PROSPECT' = Input_PROSPECT))
}

#' Complete the list of PROSPECT parameters with default values
#'
#' @param urldb character. URL for online repository where to download data
#' @param dbName character. name of the database available online
#'
#' @return list. Includes leaf chemistry, refl, tran & number of samples
#' @export

download_LeafDB <- function(urldb = NULL,
                            dbName = 'ANGERS'){
  # repository where data are stored
  if (is.null(urldb)) urldb <- 'https://gitlab.com/jbferet/myshareddata/raw/master/LOP/'
  # download leaf chemistry and optical properties
  DataBioch <- data.table::fread(file.path(urldb,dbName,'DataBioch.txt'))
  Refl <- data.table::fread(file.path(urldb,dbName,'ReflectanceData.txt'))
  Tran <- data.table::fread(file.path(urldb,dbName,'TransmittanceData.txt'))
  # Get wavelengths corresponding to the reflectance & transmittance measurements
  lambda <- Refl$wavelength
  Refl$wavelength <- Tran$wavelength <- NULL
  # Get the number of samples
  nbSamples <- ncol(Refl)
  return(list('DataBioch' = DataBioch,
              'lambda' = lambda,
              'Refl' = Refl,
              'Tran' = Tran,
              'nbSamples' = nbSamples))
}


#' Complete the list of PROSPECT parameters with default values
#'
#' @param Input_PROSPECT input parameters sent to PROSPECT by user
#' @param Parm2Add Parameters to be added to input parameters
#' @param ExpectedParms full set of parameters expected to run PROSPECT
#'
#' @return Input_PROSPECT
#' @export

Complete_Input_PROSPECT <- function(Input_PROSPECT, Parm2Add, ExpectedParms) {
  ii <- 0
  nbSamples <- length(Input_PROSPECT[[1]])
  nbInputs <- length(Input_PROSPECT)
  for (i in Parm2Add) {
    ii <- ii + 1
    nbInputs <- nbInputs + 1
    Input_PROSPECT[[nbInputs]] <- matrix(ExpectedParms[[i]],
                                         ncol = 1, nrow = nbSamples)
    names(Input_PROSPECT)[[nbInputs]] <- names(ExpectedParms)[[i]]
  }
  return(data.frame('CHL' = matrix(Input_PROSPECT$CHL, ncol = 1),
                    'CAR' = matrix(Input_PROSPECT$CAR, ncol = 1),
                    'ANT' = matrix(Input_PROSPECT$ANT, ncol = 1),
                    'BROWN' = matrix(Input_PROSPECT$BROWN, ncol = 1),
                    'EWT' = matrix(Input_PROSPECT$EWT, ncol = 1),
                    'LMA' = matrix(Input_PROSPECT$LMA, ncol = 1),
                    'PROT' = matrix(Input_PROSPECT$PROT, ncol = 1),
                    'CBC' = matrix(Input_PROSPECT$CBC, ncol = 1),
                    'N' = matrix(Input_PROSPECT$N, ncol = 1),
                    'alpha' = matrix(Input_PROSPECT$alpha, ncol = 1)))
}


#' Convert plain text colour to ANSI code
#'
#' @param colour colour in plain text ("red", "green", etc.) to convert to ANSI
#'
#' @return string representing provided colour as ANSI encoding
#'
#' @examples
#' colour_to_ansi("red") # gives: "\033[31m"
#' @export

colour_to_ansi <- function(colour) {
  # Note ANSI colour codes
  colour_codes <- list("black" = 30,
                       "red" = 31,
                       "green" = 32,
                       "yellow" = 33,
                       "blue" = 34,
                       "magenta" = 35,
                       "cyan" = 36,
                       "white" = 37)

  # Create ANSI version of colour
  ansi_colour <- paste0("\033[", colour_codes[[colour]], "m")
  return(ansi_colour)
}
jbferet/prospect documentation built on June 12, 2025, 10:05 a.m.