In jbferet/prospect: PROSPECT leaf radiative transfer model and inversion routines

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval=FALSE
)

`prospect` inversion: generalities

PROSPECT can be inverted using iterative optimization, by calling the function Invert_PROSPECT. This iterative optimization is based on the function fmincon included in the package pracma.

By default, the merit function used for the inversion named Merit_PROSPECT_RMSE minimizes the RMSE between the simulated and the measured leaf optical properties. However, users can define their own merit function with associated criterion to minimize and define it as input variable for Invert_PROSPECT: MeritFunction = "MyOwnMeritFunction".

As optimal spectral domains were identified when inverting PROSPECT with Merit_PROSPECT_RMSE, the function Invert_PROSPECT_OPT does not accept other merit functions. However, users can easily develop their own script to identify optimal spectral domains with the merit function of their choice.

The results of inversion will also depend on the parameters to be assessed during the prospect model. The following table shows the correspondence between model versions available in prospect and input parameters accepted. Each parameter can be considered as a parameter to estimate or a parameter with a value set to a fixed value during inversion. The choice depends on different factors, such as prior knowledge, or spectral domain considered (e.g.: as pigments do not absorb in the SWIR domain, they cannot be assessed when inverting PROSPECT with SWIR data. Hence, user can set their value with no consequence for inversion).

Input variables

The function Invert_PROSPECT requires either reflectance, or transmittance, or both. User can define which input variables of PROSPECT should be assessed during inversion, and which ones should be set to a given value (0 or user's choice). The list of input variables for inversion is :

SpecPROSPECT: data frame including the spectral bands, refractive index and specific absorption coefficients. SpecPROSPECT = NULL uses the spectral domain from 400 nm to 2500 nm. Simulation and inversion on different spectral domains can be performed by adjusting SpecPROSPECT
Refl: numeric: individual leaf reflectance corresponding to the spectral domain defined in SpecPROSPECT. Set to NULL if inversion on transmittance only
Tran: numeric: individual leaf transmittance corresponding to the spectral domain defined in SpecPROSPECT. Set to NULL if inversion on reflectance only
Parms2Estimate list. Parameters to estimate. Set to 'ALL' by default.
InitValues data frame including initial values of PROSPECT input parameters. During optimization, they are used either as initialization values for parameters to estimate, or as fix values for other parameters. Parameters are not taken into account when not compatible with PROSPECT_version.
PROSPECT_version character. Corresponds to the PROSPECT version. Should be one of the following versions: 'D', 'PRO'.
MeritFunction character. name of the function to be used as merit function with given criterion to minimize (default = Merit_PROSPECT_RMSE)
xlub data frame. Boundaries of the parameters to estimate. The data frame must have columns corresponding to \code{Parms2Estimate} first line being the lower boundaries and second line the upper boundaries.
Est_Brown_Pigments boolean. Should BROWN be assessed or not? Brown pigments are found during later stages of senescent leaves, but they are usually not found in juvenile, mature or early senescent leaves.
Est_alpha boolean. Should alpha be assessed or not? Keep in mind that most published results use alpha with its default value.

Output variables

Invert_PROSPECT returns a list of assessed PROSPECT parameters defined in Parms2Estimate.

PROSPECT-D inversion: using full spectral information

Inversion over the full spectral domain

All parameters are assessed, except alpha and BROWN which are set to their default value. Brown pigments are characteristic of late senescent stages. They are absent for all other development stages, hence the assessment of brown pigments is not systematic here.

# simulate leaf optical properties
LRT_D <- PROSPECT(CHL = 45, CAR = 10, ANT = 0.2, 
                  EWT = 0.012, LMA = 0.01, N = 1.3)
# define set of parameters to be assessed
Parms2Estimate  <- 'ALL'
# invert PROSPECT with simulated leaf optical properties
OutPROSPECT1 <- Invert_PROSPECT(Refl = LRT_D$Reflectance, 
                                Tran = LRT_D$Transmittance,
                                Parms2Estimate = Parms2Estimate, 
                                PROSPECT_version = 'D')


# simulate leaf optical properties with brown pigments
LRT_D_brown <- PROSPECT(CHL = 45, CAR = 10, ANT = 0.2, BROWN = 0.1, 
                        EWT = 0.012, LMA = 0.01, N = 1.3)
# define set of parameters to be assessed
Parms2Estimate  <- 'ALL'
# invert PROSPECT with simulated leaf optical properties and specify that brown 
# pigment should be assessed in addition to all other constituents
OutPROSPECT2 <- Invert_PROSPECT(Refl = LRT_D_brown$Reflectance, 
                                Tran = LRT_D_brown$Transmittance,
                                Parms2Estimate = Parms2Estimate, 
                                PROSPECT_version = 'D', 
                                Est_Brown_Pigments = TRUE)

Inversion over the VNIR domain

All parameters are assessed, except alpha and BROWN which are set to their default value. FitSpectralData directly adjusts reflectance, transmittance and optical constants to user defined spectral bands. These spectral bands should be integer values between 400 and 2500, and can be continuous or not.

# define spectral subdomain from 400 nm to 800 nm
SpectralSubDomain <- seq(400,800)
# adjust spectral domain for reflectance, transmittance and SpecPROSPECT
SubData <- FitSpectralData(lambda = LRT_D$wvl, 
                           Refl = LRT_D$Reflectance, 
                           Tran = LRT_D$Transmittance,
                           UserDomain = SpectralSubDomain)

# Note that FitSpectralData can also run with UserDomain defining 
# upper and lower bounds when setting input argument 'UL_Bounds = TRUE'
ULBounds <- c(400,800)
SubData <- FitSpectralData(lambda = LRT_D$wvl, 
                           Refl = LRT_D$Reflectance, 
                           Tran = LRT_D$Transmittance,
                           UserDomain = ULBounds, 
                           UL_Bounds = TRUE)

# invert PROSPECT with simulated leaf optical properties
OutPROSPECT3 <- Invert_PROSPECT(SpecPROSPECT = SubData$SpecPROSPECT, 
                                Refl = SubData$Refl, 
                                Tran = SubData$Tran,
                                Parms2Estimate = Parms2Estimate, 
                                PROSPECT_version = 'D')

Inversion over the VNIR domain with LMA and EWT values set by user

Only pigments and N are assessed. The same optical properties as previous example are used.

# define set of parameters to be assessed
Parms2Estimate <- c('CHL', 'CAR', 'ANT', 'N')
# define set of parameters to be assessed
InitValues <- c('EWT' = 0.01, 'LMA' = 0.01)
# invert PROSPECT with simulated leaf optical properties
OutPROSPECT4 <- Invert_PROSPECT(SpecPROSPECT = SubData$SpecPROSPECT, 
                                Refl = SubData$Refl, 
                                Tran = SubData$Tran,
                                Parms2Estimate = Parms2Estimate, 
                                PROSPECT_version = 'D', 
                                InitValues = InitValues)

Inversion over the SWIR domain between 1700 nm and 2400 nm

EWT, LMA and N are assessed. The same optical properties as previous example are used.

# define set of parameters to be assessed
Parms2Estimate <- c('EWT', 'LMA', 'N')
# define spectral subdomain 
SpectralSubDomain <- seq(1700,2400)
# adjust spectral domain
SubData <- FitSpectralData(lambda = LRT_D$wvl, 
                           Refl = LRT_D$Reflectance, 
                           Tran = LRT_D$Transmittance,
                           UserDomain = SpectralSubDomain)

# invert PROSPECT with simulated leaf optical properties
OutPROSPECT5 <- Invert_PROSPECT(SpecPROSPECT = SubData$SpecPROSPECT, 
                                Refl = SubData$Refl, 
                                Tran = SubData$Tran,
                                Parms2Estimate = Parms2Estimate, 
                                PROSPECT_version = 'D')

PROSPECT-D inversion: using optimal spectral domains

The function Invert_PROSPECT_OPT automatically sets the optimal spectral domains during inversion for all constituents to be assessed.

Optimal spectral domains and configuration are defined in Féret et al. (2019), Féret et al. (2021), and Spafford et al. (2021).

N does not need to be part of Parms2Estimate, as it is automatically assessed when needed.

# define set of parameters to be assessed
Parms2Estimate  <- c('CHL','CAR','ANT','EWT','LMA')
# call Invert_PROSPECT_OPT to get optimal estimation of leaf parameters 
# using optimal spectral domains
ParmEst <- Invert_PROSPECT_OPT(lambda = LRT_D$wvl, 
                               Refl = LRT_D$Reflectance,
                               Tran = LRT_D$Transmittance,
                               PROSPECT_version = 'D',
                               Parms2Estimate = Parms2Estimate)

PROSPECT-PRO inversion

Such definition of optimal spectral domains can also be set manually. For example, here is how to estimate protein content from leaf optical properties using the optimal spectral domain defined in Féret et al. (2021).

Please note that N needs to be added to Parms2Estimate, if user want it to be assessed during the inversion, otherwise it will be set to its default value.

# simulate leaf optical properties
LRT_PRO <- PROSPECT(CHL = 45, CAR = 10, ANT = 0.2, EWT = 0.012,
                    PROT = 0.002, CBC = 0.015, N = 1.3)
# define spectral subdomain 
SpectralSubDomain <- c(2125,2175)
# adjust spectral domain
SubData <- FitSpectralData(lambda = LRT_PRO$wvl, 
                           Refl = LRT_PRO$Reflectance, 
                           Tran = LRT_PRO$Transmittance,
                           UserDomain = SpectralSubDomain, 
                           UL_Bounds = TRUE)

Parms2Estimate <- c('EWT', 'PROT', 'CBC', 'N')
# invert PROSPECT with simulated leaf optical properties
OutPROSPECT6 <- Invert_PROSPECT(SpecPROSPECT = SubData$SpecPROSPECT,
                                Refl = SubData$Refl, 
                                Tran = SubData$Tran, 
                                Parms2Estimate = Parms2Estimate, 
                                PROSPECT_version = 'PRO')

PROSPECT inversion with user-defined merit function

A new merit function can be defined, as long as it uses the same input and output variables as the function Merit_PROSPECT_RMSE.

Here is an example of a function using nMRSE instead of RMSE for the merit function. This means that for each spectral band, the square difference between measured and simulated leaf optics is normalized by the value of the measured leaf optics.

First, the merit function is

Merit_PROSPECT_nRMSE <- function(x, SpecPROSPECT,
                                 Refl, Tran,
                                 Input_PROSPECT,
                                 Parms2Estimate) {
  x[x < 0] <- 0
  Input_PROSPECT[Parms2Estimate] <- x
  RT <- do.call('PROSPECT', c(list(SpecPROSPECT = SpecPROSPECT), Input_PROSPECT))
  fcr <- fct <- 0
  if (!is.null(Refl)) fcr <- sqrt(sum(((Refl - RT$Reflectance)**2)/Refl) / length(RT$Reflectance))
  if (!is.null(Tran)) fct <- sqrt(sum(((Tran - RT$Transmittance)**2)/Tran) / length(RT$Transmittance))
  fc <- fcr + fct
  return(fc)
}

# simulate leaf optical properties
LRT <- PROSPECT(CHL = 45, CAR = 10, ANT = 0.2, 
                EWT = 0.02, LMA = 0.015, N = 1.8)
Parms2Estimate <- c('CHL', 'CAR', 'ANT', 'EWT', 'LMA', 'N')
# invert PROSPECT with simulated leaf optical properties
OutPROSPECT7 <- Invert_PROSPECT(Refl = LRT$Refl, 
                                Tran = LRT$Tran, 
                                Parms2Estimate = Parms2Estimate, 
                                PROSPECT_version = 'D', 
                                MeritFunction = 'Merit_PROSPECT_nRMSE')