The PEcAn RTM package
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 = gfortranOverview 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).
##
## Attaching package: 'PEcAnRTM'## The following object is masked from 'package:graphics':
##
## matplotwl <- 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.2788870Notice 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.
## [1] "spectra" "matrix"This class allows traditional subsetting via index with single
brackets ([]), as well as subsetting by wavelength with
double brackets.
## 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## 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.07136411The 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.)
## 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.000p4s <- 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.
## 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=1Inversion
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.
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.
## function (params, seed = NULL)
## prospect(params, 5)[, 1]
## <bytecode: 0x634e6a95fa70>
## <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$nchains <- 2 # Number of MCMC chains
invert.options$ngibbs.max <- 5e4 # Maximum number of iterations per chain (fails if no convergence by then)
invert.options$do.lsq <- TRUE # Initialize with results from a fast least-squares optimization algorithm
invert.options$threshold <- 1.3 # Maximum value for Gelman diagnosticThe full list of inversion options is below. See the documentation
for invert.auto for a full description of each of
these.
## [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" "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")
}## Loading required namespace: minpack.lm## [1] "Running chain 1 of 2"## Loading required namespace: MASS## [1] "Running chain 2 of 2"
## [1] "The following parameters did not converge: Cw (1.419), Cm (1.752)"
## [1] "The following parameters did not converge: Car (1.335), Cm (1.392), residual (1.399)"
## [1] "The following parameters did not converge: Cm (1.418), residual (1.560)"
## [1] "The following parameters did not converge: residual (1.595)"
## [1] "Converged with all Gelman diag <= 1.017"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.080740 1.169298 1.157836 1.050614 1.223484 1.536570 FALSE
## 47 1.044447 1.125787 1.159229 1.023841 1.122472 1.405202 FALSE
## 48 1.027697 1.145893 1.172236 1.019656 1.097492 1.357509 FALSE
## 49 1.016517 1.159094 1.183522 1.012208 1.067158 1.256458 FALSE
## 50 1.029338 1.235793 1.160377 1.025443 1.096639 1.268767 FALSE
## 51 1.026910 1.229342 1.164557 1.016737 1.091276 1.256598 FALSE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 TRUE TRUE FALSE TRUE
## 47 TRUE TRUE FALSE TRUE
## 48 TRUE TRUE FALSE FALSE
## 49 TRUE TRUE FALSE FALSE
## 50 TRUE TRUE FALSE FALSE
## 51 TRUE TRUE FALSE FALSE
## PSRF residual > 1.10
## 46 TRUE
## 47 TRUE
## 48 TRUE
## 49 TRUE
## 50 TRUE
## 51 TRUE
## [1] "Converged with all Gelman diag <= 1.004"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.028729 1.230528 1.163625 1.016446 1.092404 1.258088 FALSE
## 47 1.030283 1.251963 1.174008 1.014015 1.066117 1.183010 FALSE
## 48 1.021176 1.273824 1.175712 1.005429 1.036672 1.156235 FALSE
## 49 1.015635 1.285901 1.184449 1.005218 1.030687 1.109799 FALSE
## 50 1.007945 1.294600 1.178792 1.010350 1.032911 1.072760 FALSE
## 51 1.003906 1.291644 1.168939 1.010051 1.018304 1.049877 FALSE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 TRUE TRUE FALSE FALSE
## 47 TRUE TRUE FALSE FALSE
## 48 TRUE TRUE FALSE FALSE
## 49 TRUE TRUE FALSE FALSE
## 50 TRUE TRUE FALSE FALSE
## 51 TRUE TRUE FALSE FALSE
## PSRF residual > 1.10
## 46 TRUE
## 47 TRUE
## 48 TRUE
## 49 TRUE
## 50 FALSE
## 51 FALSE
## [1] "The following parameters did not converge: Cab (1.374)"
## [1] "Converged with all Gelman diag <= 1.005"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.003888 1.375278 1.121231 1.017124 1.012614 1.012503 FALSE
## 47 1.016331 1.365004 1.087012 1.015157 1.013974 1.005970 FALSE
## 48 1.031465 1.342348 1.076987 1.015742 1.036220 1.002120 FALSE
## 49 1.037575 1.360289 1.058824 1.016697 1.034113 1.001785 FALSE
## 50 1.036279 1.282889 1.047067 1.022808 1.031465 1.006187 FALSE
## 51 1.007447 1.277659 1.045203 1.028888 1.024956 1.005488 FALSE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 TRUE TRUE FALSE FALSE
## 47 TRUE FALSE FALSE FALSE
## 48 TRUE FALSE FALSE FALSE
## 49 TRUE FALSE FALSE FALSE
## 50 TRUE FALSE FALSE FALSE
## 51 TRUE FALSE FALSE FALSE
## PSRF residual > 1.10
## 46 FALSE
## 47 FALSE
## 48 FALSE
## 49 FALSE
## 50 FALSE
## 51 FALSE
## [1] "Converged with all Gelman diag <= 1.013"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.033448 1.301569 1.045310 1.021594 1.030594 1.003962 FALSE
## 47 1.010808 1.280539 1.044255 1.027388 1.025542 1.006440 FALSE
## 48 1.006783 1.247267 1.038093 1.039421 1.032377 1.015311 FALSE
## 49 1.006808 1.257933 1.043904 1.026112 1.038902 1.017236 FALSE
## 50 1.009356 1.218080 1.035193 1.021091 1.049528 1.026883 FALSE
## 51 1.012987 1.191433 1.026135 1.024215 1.048546 1.040513 FALSE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 TRUE FALSE FALSE FALSE
## 47 TRUE FALSE FALSE FALSE
## 48 TRUE FALSE FALSE FALSE
## 49 TRUE FALSE FALSE FALSE
## 50 TRUE FALSE FALSE FALSE
## 51 TRUE FALSE FALSE FALSE
## PSRF residual > 1.10
## 46 FALSE
## 47 FALSE
## 48 FALSE
## 49 FALSE
## 50 FALSE
## 51 FALSE
## [1] "Converged with all Gelman diag <= 1.014"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.011753 1.263190 1.040519 1.022794 1.045857 1.021908 FALSE
## 47 1.012837 1.171759 1.030761 1.022459 1.045332 1.031823 FALSE
## 48 1.009181 1.188760 1.022057 1.019581 1.051695 1.041481 FALSE
## 49 1.006556 1.223649 1.016482 1.020675 1.055084 1.073490 FALSE
## 50 1.013542 1.183314 1.015976 1.018780 1.065703 1.065061 FALSE
## 51 1.029480 1.107558 1.014456 1.019019 1.081059 1.038339 FALSE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 TRUE FALSE FALSE FALSE
## 47 TRUE FALSE FALSE FALSE
## 48 TRUE FALSE FALSE FALSE
## 49 TRUE FALSE FALSE FALSE
## 50 TRUE FALSE FALSE FALSE
## 51 TRUE FALSE FALSE FALSE
## PSRF residual > 1.10
## 46 FALSE
## 47 FALSE
## 48 FALSE
## 49 FALSE
## 50 FALSE
## 51 FALSE
## [1] "Converged with all Gelman diag <= 1.004"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.007250 1.194474 1.016018 1.020884 1.058372 1.074005 FALSE
## 47 1.011772 1.161997 1.015045 1.022007 1.068776 1.049461 FALSE
## 48 1.029366 1.105246 1.014017 1.017895 1.081555 1.036914 FALSE
## 49 1.084736 1.076918 1.011814 1.008438 1.129401 1.033577 FALSE
## 50 1.122661 1.058782 1.011729 1.009151 1.164941 1.029861 TRUE
## 51 1.102884 1.042066 1.012917 1.015952 1.184660 1.004136 TRUE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 TRUE FALSE FALSE FALSE
## 47 TRUE FALSE FALSE FALSE
## 48 TRUE FALSE FALSE FALSE
## 49 FALSE FALSE FALSE TRUE
## 50 FALSE FALSE FALSE TRUE
## 51 FALSE FALSE FALSE TRUE
## PSRF residual > 1.10
## 46 FALSE
## 47 FALSE
## 48 FALSE
## 49 FALSE
## 50 FALSE
## 51 FALSE
## [1] "Converged with all Gelman diag <= 1.002"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.085195 1.076856 1.011815 1.008319 1.129644 1.033451 FALSE
## 47 1.110924 1.055592 1.011984 1.009568 1.160053 1.025880 TRUE
## 48 1.106289 1.040648 1.012061 1.015059 1.186917 1.005648 TRUE
## 49 1.107562 1.043891 1.009146 1.003551 1.158366 1.010693 TRUE
## 50 1.111014 1.045049 1.015131 1.001583 1.131815 1.007262 TRUE
## 51 1.133145 1.024904 1.008148 1.002403 1.132599 1.004778 TRUE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 FALSE FALSE FALSE TRUE
## 47 FALSE FALSE FALSE TRUE
## 48 FALSE FALSE FALSE TRUE
## 49 FALSE FALSE FALSE TRUE
## 50 FALSE FALSE FALSE TRUE
## 51 FALSE FALSE FALSE TRUE
## PSRF residual > 1.10
## 46 FALSE
## 47 FALSE
## 48 FALSE
## 49 FALSE
## 50 FALSE
## 51 FALSE
## [1] "Converged with all Gelman diag <= 1.000"## N Cab Car Cw Cm residual PSRF N > 1.10
## 46 1.108690 1.043450 1.008282 1.005284 1.169422 1.0110327 TRUE
## 47 1.111353 1.043898 1.014220 1.001502 1.140563 1.0066172 TRUE
## 48 1.133271 1.025959 1.008141 1.002512 1.134123 1.0050223 TRUE
## 49 1.137303 1.016542 1.009298 1.002279 1.112392 1.0016647 TRUE
## 50 1.130122 1.023832 1.010357 1.004909 1.093035 0.9999536 TRUE
## 51 1.128091 1.027770 1.005208 1.007650 1.080391 1.0002770 TRUE
## PSRF Cab > 1.10 PSRF Car > 1.10 PSRF Cw > 1.10 PSRF Cm > 1.10
## 46 FALSE FALSE FALSE TRUE
## 47 FALSE FALSE FALSE TRUE
## 48 FALSE FALSE FALSE TRUE
## 49 FALSE FALSE FALSE TRUE
## 50 FALSE FALSE FALSE FALSE
## 51 FALSE FALSE FALSE FALSE
## PSRF residual > 1.10
## 46 FALSE
## 47 FALSE
## 48 FALSE
## 49 FALSE
## 50 FALSE
## 51 FALSE
## [1] "Converged with all Gelman diag <= 1.009"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(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)