Simulating reflectance

Available radiative transfer models are called by passing a vector of their parameters. Similar models with different versions (e.g. PROSPECT 4, 5, 5B) also have a “version” argument. These models return a matrix of reflectance, transmittance, and/or absorption spectra (depending on the model) at 1 nm resolution over the wavelength range 400 to 2500 nm.

The PROSPECT family of models returns the reflectance (column 1) and transmittance (column 2) for an individual leaf as a function of 4 to 6 parameters (depending on the version).

library(PEcAnRTM)

## 
## Attaching package: 'PEcAnRTM'

## The following object is masked from 'package:graphics':
## 
##     matplot

wl <- 400:2500
params <- c(
  "N" = 1.4, "Cab" = 40, "Car" = 15,
  "Cbrown" = 0.5, "Cw" = 0.002, "Cm" = 0.004
)

p4 <- prospect(params[c(-3, -4)], version = 4)
p4

##     reflectance transmittance
## 400  0.04099898  0.0001814716
## 401  0.04103723  0.0001699630
## 402  0.04106248  0.0001575602
## 403  0.04111370  0.0001452140
## 404  0.04120401  0.0001348522
## ...
##      reflectance transmittance
## 2496   0.1422667     0.2801837
## 2497   0.1417965     0.2792691
## 2498   0.1420623     0.2795135
## 2499   0.1413868     0.2782765
## 2500   0.1418849     0.2788870

Notice only the first and last few lines of the spectra are printed, and the wavelength is annotated on the side. This is is because all RTM simulations in this package are of a special matrix class spectra.

class(p4)

## [1] "spectra" "matrix"

This class allows traditional subsetting via index with single brackets ([]), as well as subsetting by wavelength with double brackets.

p4[1:50, 1]

##     reflectance
## 400  0.04099898
## 401  0.04103723
## 402  0.04106248
## 403  0.04111370
## 404  0.04120401
## ...
##     reflectance
## 445  0.04545556
## 446  0.04542921
## 447  0.04540287
## 448  0.04537665
## 449  0.04533726

p4[[500:520, 2]]

##     transmittance
## 500   0.008937186
## 501   0.010279518
## 502   0.011772740
## 503   0.013438540
## 504   0.015276985
## ...
##     transmittance
## 516    0.05336301
## 517    0.05776873
## 518    0.06224375
## 519    0.06679217
## 520    0.07136411

The package also provides special plotting (plot, matplot) and combining (cbind) methods. Note that the cbind method automatically matches wavelengths, which facilitates working with spectra with different wavelengths.

p5 <- prospect(params[-4], version = 5)
p5b <- prospect(params, version = "5B")
p_multi <- cbind(p4, p5, p5b)
matplot(p_multi[, c(1, 3, 5)], lty = 1:3, col = 2, ylim = c(0, 1))
matplot(1 - p_multi[, c(2, 4, 6)], lty = 1:3, col = 3,  add = TRUE)
legend("topright", c("Reflectance", "Transmittance"), col=c(2, 3), lty = 1)
legend("top", c("4", "5", "5B"), lty = 1:3)

The SAIL family of models returns the bidirectional (1), hemispherical directional (2), directional hemispherical (3), and bidirectional hemispherical (4) reflectance factors for a canopy with a given set of approximately 20 parameters. It is often coupled to the PROSPECT model as PRO4SAIL. (Again, note that the return type is a spectra, which leads matplot to automatically put wavelengths on the x axis.)

sail.params <- defparam("pro4sail")
print(sail.params)

##        N      Cab      Car   Cbrown       Cw       Cm    LIDFa    LIDFb 
##    1.500   40.000    8.000    0.000    0.010    0.009   -0.350   -0.150 
## TypeLIDF      LAI        q      tts      tto      psi    psoil 
##    1.000    3.000    0.010   30.000   10.000    0.000    1.000

p4s <- pro4sail(sail.params)
matplot(p4s, xlab="Wavelength (nm)", ylab="Reflectance")
legend("topright", as.character(1:4), col=1:4, lty=1:4)

The above example illustrates the use of defparam to get the default parameters for a particular model. Similarly, model.list is a data.frame containing all currently available models.

print(model.list)

##          modname modcode         fullname
## 1     prospect_d    1153       PROSPECT D
## 2    prospect_5b    1152      PROSPECT 5B
## 3     prospect_5    1151       PROSPECT 5
## 4     prospect_4    1141       PROSPECT 4
## 5       pro4sail    2111         PRO4SAIL
## 6      pro4saild    2113        PRO4SAILD
##                                                                          par.names
## 1                                                     N Cab Car Canth Cbrown Cw Cm
## 2                                                           N Cab Car Cbrown Cw Cm
## 3                                                                  N Cab Car Cw Cm
## 4                                                                      N Cab Cw Cm
## 5            N Cab Car Cbrown Cw Cm LIDFa LIDFb TypeLIDF_r LAI q tts tto psi psoil
## 6      N Cab Car Canth Cbrown Cw Cm LIDFa LIDFb TypeLIDF_r LAI q tts tto psi psoil
##                                                                                                                         par.default
## 1                                                                              N=1.5 Cab=40 Car=8 Canth=8 Cbrown=0 Cw=0.01 Cm=0.009
## 2                                                                                      N=1.5 Cab=40 Car=8 Cbrown=0 Cw=0.01 Cm=0.009
## 3                                                                                               N=1.5 Cab=40 Car=8 Cw=0.01 Cm=0.009
## 4                                                                                                     N=1.5 Cab=40 Cw=0.01 Cm=0.009
## 5          N=1.5 Cab=40 Car=8 Cbrown=0 Cw=0.01 Cm=0.009 LIDFa=-0.35 LIDFb=-0.15 TypeLIDF=1 LAI=3 q=0.01 tts=30 tto=10 psi=0 psoil=1
## 6  N=1.5 Cab=40 Car=8 Canth=1 Cbrown=0 Cw=0.01 Cm=0.009 LIDFa=-0.35 LIDFb=-0.15 TypeLIDF=1 LAI=3 q=0.01 tts=30 tto=10 psi=0 psoil=1

Inversion

A novel feature of this package is the ability to perform a Bayesian inversion of a Radiative Transfer Model. Here are several advantages of the Bayesian approach:

• Parameter uncertainty: The output of a Bayesian analysis is a full joint probability distribution of the model parameters, which includes a robust estimate of their uncertainy and covariance between parameters.
• Prior knowledge: If previous, independent estimates of parameters are available, these can be used to inform the model.
• Partitioning variability: Random effects models provide a powerful framework for understanding the sources of variability and uncertainty in a data set.

An inversion can be performed using the invert.auto function, which uses the Metropolis Hastings MCMC algorithm to invert an arbitrary model. There are a lot of configuration options to invert.auto, so the recommended way to perform an inversion is to start with a default settings list, provided by the package itself. In the following sample, we demonstrate the default inversion of the PROSPECT model.

invert.options <- default.settings.prospect

The model to be inverted is, in this case, just a one-line call to the PROSPECT 5 model with the params vector. It returns a vector of reflectance values. The requirement is that the “model” output be as long as the length (or number of rows) of the observation vector (or matrix). In this case, the PROSPECT model returns 2101 reflectance values, so our observation also has to have this many points.

invert.options$model

## function (params, seed = NULL) 
## prospect(params, 5)[, 1]
## <bytecode: 0x556bb29d6cc0>
## <environment: namespace:PEcAnRTM>

The recommended way to set settings is to start with a default object and modify it. For example, in the following block, we reduce the length of the run to make it go a little faster. threshold = ... sets the maximum value of the multivariate Gelman Diagnostic used to assess convergence – it defaults to 1.1, but we set it higher here for demonstrative purposes).

invert.options$n.tries <- 1      # Number of attempts
invert.options$nchains <- 2      # Number of MCMC chains
invert.options$ngibbs <- 5000    # Number of iterations per chain
invert.options$burnin <- 1000    # Length of burnin period
invert.options$do.lsq.first <- TRUE # Initialize with results from a fast least-squares optimization algorithm
invert.options$threshold <- 1.3  # Maximum value for Gelman diagnostic

The full list of inversion options is below. See the documentation for invert.auto for a full description of each of these.

names(invert.options)

##  [1] "model"          "inits.function" "prior.function" "param.mins"    
##  [5] "param.maxs"     "ngibbs"         "nchains"        "burnin"        
##  [9] "ngibbs.max"     "ngibbs.min"     "ngibbs.step"    "return.samples"
## [13] "target"         "do.lsq"         "save.samples"   "quiet"         
## [17] "adapt"          "adj_min"        "n.tries"        "do.lsq.first"  
## [21] "threshold"

Now, we load some test data (for Acer rubrum leaves (testspec_ACRU, accessible through data(testspec)).

data(testspec)
observed <- testspec_ACRU[,1]
plot(wl, observed, xlab="Wavelength", ylab="Reflectance", type='l')

To perform an inversion, just pass the observation matrix and the inversion settings into invert.auto. quiet = TRUE suppresses the progress bar, which is done here to clean up the knitted document. We also set parallel = FALSE here because this vignette is rebuilt on shared hardware where only one processor might be available; for runs on your own machine, you will probably want to set it TRUE to use all the cores on your machine.

if(file.exists("inversion.output.rds")){
    inversion.output <- readRDS("inversion.output.rds")
} else {
    inversion.output <- invert.auto(observed = observed,
                                    invert.options = invert.options,
                                    quiet = TRUE,
                                    parallel = FALSE)
    saveRDS(inversion.output, "inversion.output.rds")
}

## [1] "Running chain 1 of 2"

## Loading required namespace: MASS

## [1] "Running chain 2 of 2"
## [1] "The following parameters did not converge: N (4.474), Cab (46.096), Car (12.711), Cw (2.077), Cm (2.083), residual (34.387), mpsrf (14.114)"
## [1] "The following parameters did not converge: N (3.794), Cab (9.790), Car (6.178), Cw (1.905), Cm (1.441), residual (17.537), mpsrf (2.584)"
## [1] "The following parameters did not converge: N (2.550), Cab (7.033), Car (2.153), Cw (1.407), residual (10.135), mpsrf (1.588)"
## [1] "The following parameters did not converge: N (1.802), Cab (4.856), Cw (1.441), residual (7.295), mpsrf (1.411)"
## [1] "The following parameters did not converge: N (1.345), Cab (3.041), residual (5.119)"
## [1] "The following parameters did not converge: Cab (1.709), residual (3.268)"
## [1] "Converged with all Gelman diag <= 1.000"

The output of invert.auto is a list of two objects: A list of summary statistics for each parameter and an mcmc.list object of the samples for diagnostic purposes or calculation of other summary statistics.

par(mfrow=c(2,2))
plot(inversion.output$samples, auto.layout=FALSE)

par(mfrow=c(1,1))
samples.mat <- as.matrix(inversion.output$samples)[-(2000:0),1:5]
colnames(samples.mat) <- params.prospect5
pairs(samples.mat, pch=".")

means <- unlist(inversion.output$results[grep("mu", names(inversion.output$results))])[1:5]
prospect.sim <- prospect(means, 5)[,1]  # reflectance

plot(wl, observed, type='l', col=1, xlab="wavelength (nm)", ylab="reflectance")
lines(wl, prospect.sim, type='l', col=2)
legend("topright", c("observed", "predicted"), lty=1, col=1:2)