Introduction

The PEcAnRTM package provides tools for analyses involving common radiative transfer models. The highlights of this package are its ability to efficiently run a suite of related leaf and canopy radiative transfer models, as well as to perform maximum likelihood and, particularly, Bayesian inversions of these models.

Installation

The easiest way to install this package is via install_github from the devtools package.

install.packages("devtools")
library(devtools)
install_github("ashiklom/pecan", subdir = "modules/rtm")
# Defaults to branch 'master'. 
# For custom branches, add `ref = "branchname"`

This package relies on a modern (>= 2003) Fortran compiler, as determined by your R installation. On most Unix systems, this is standard, but R may specify a particular compiler version that you don’t have, resulting in an installation error. To fix this, simply add the following to your system ~/R/Makevars file.

FC = gfortran

Overview of features

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)
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

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)

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

p4[1:50, 1]
p4[[500:520, 2]]

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)
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)

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

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.

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

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)

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 (NOTE: quiet = TRUE suppresses the progress bar, which is done here to clean up the knitted document. 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 demosntrative purposes).

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,
                                    threshold = 1.3)
    saveRDS(inversion.output, "inversion.output.rds")
}

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,1))
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)