← Back to RTM-Suite

This is the SCOPEinR-specific companion to How-in-R.Rmd (which covers ToolsRTM’s leaf/canopy models alone). SCOPE (Soil Canopy Observation, Photochemistry and Energy fluxes) couples soil reflectance, canopy optical BRDF, leaf photosynthesis, and a full energy-balance closure loop into one simulation – given one LUT input row (leaf/canopy/meteo/soil parameters), get.SCOPE() returns TOC reflectance, converged leaf/soil temperatures, energy fluxes, and (optionally) canopy-level fluorescence.

1. Load the packages

library(ToolsRTM)
library(SCOPEinR)

2. Run one full SCOPE simulation

Load SCOPE’s own bundled example LUT row and model options, and run the full leaf optics -> soil -> optical BRDF -> energy balance -> fluorescence pipeline in one get.SCOPE() call.

table.with.opts <- read.table(system.file("input", "setoptions.csv", package = "SCOPEinR"),
                               header = TRUE, sep = ",")
Table.LUT <- read.table(system.file("input", "LUT_input.csv", package = "SCOPEinR"),
                         header = TRUE, sep = ",")

invisible(capture.output(
  db.sim <- SCOPEinR::get.SCOPE(
    LUT           = Table.LUT[1, ],
    options.SCOPE = table.with.opts,
    optipar       = SCOPEinR::optipar2021.Pro.CX,
    leaf.model    = "fluspect-CX",
    canopy.model  = "fourSAIL",
    get.outputs   = "ALL",
    get.plots     = FALSE
  )
))
res <- db.sim[[1]]  # get.SCOPE() always returns a list-of-results, even for one LUT row

cat("Canopy layers:", res$data.canopy$nlayers, "\n")
## Canopy layers: 30
cat("TOC reflectance at 550nm:", round(res$data.rad$refl[550 - 400 + 1], 4), "\n")
## TOC reflectance at 550nm: 0.0438
cat("TOC reflectance at 800nm:", round(res$data.rad$refl[800 - 400 + 1], 4), "\n")
## TOC reflectance at 800nm: 0.3588
cat("Canopy-average leaf temperature (Tcave):", round(res$data.fluxes$Tcave, 2), "degC\n")
## Canopy-average leaf temperature (Tcave): 22.04 degC
cat("Soil temperature, sunlit/shaded:", round(res$data.thermal$Tsu, 2), "/",
    round(res$data.thermal$Tsh, 2), "degC\n")
## Soil temperature, sunlit/shaded: 39.84 / 23.79 degC
cat("Net radiation, total (Rntot):", round(res$data.fluxes$Rntot, 1), "W/m2\n")
## Net radiation, total (Rntot): 496.7 W/m2
cat("Latent heat flux, total (lEtot):", round(res$data.fluxes$lEtot, 1), "W/m2\n")
## Latent heat flux, total (lEtot): 181.7 W/m2
cat("Sensible heat flux, total (Htot):", round(res$data.fluxes$Htot, 1), "W/m2\n")
## Sensible heat flux, total (Htot): 275 W/m2
cat("Canopy photosynthesis (Actot):", round(res$data.fluxes$Actot, 2), "umol m-2 s-1\n")
## Canopy photosynthesis (Actot): 19.89 umol m-2 s-1
cat("Canopy fluorescence flux (EoutF):", round(res$data.rad$EoutF, 4), "W/m2\n")
## Canopy fluorescence flux (EoutF): 0.3889 W/m2
cat("F685/F740:", round(res$data.rad$F685, 3), "/", round(res$data.rad$F740, 3), "\n")
## F685/F740: 0.311 / 1.808

A note on get.SCOPE()’s return structure: unlike toolsrtm’s Python port (a clean ScopeResult dataclass, see the Python version of this tutorial), the R version returns SCOPE’s own internal nested-list structure directly (data.rad, data.thermal, data.fluxes, data.canopy, … – 19 top-level fields). This is deliberate: it’s exactly what the original SCOPE model itself produces, so it stays directly comparable to SCOPE documentation/literature, at the cost of needing to know these field names (see below for a lookup table for the fields used in this tutorial).

What you want Where it is
TOC reflectance (full spectrum) res$data.rad$refl
Bi-hemispherical / directional components res$data.rad$rdd / rsd / rdo / rso
Wavelength grid for the above res$data.spectral$wlS (400-2400nm @ 1nm, then coarser out to 50000nm – NOT uniform 1nm throughout)
Fluorescence spectrum res$data.rad$LoF_, on res$data.spectral$wlF (640-850nm @ 1nm)
Sunlit/shaded leaf temperature per layer res$data.thermal$Tcu / Tch (length = res$data.canopy$nlayers)
Soil temperature res$data.thermal$Tsu / Tsh
Energy-balance totals res$data.fluxes$Rntot / lEtot / Htot / Gtot

3. Explore the outputs: reflectance, fluorescence, temperature, fluxes

wl_optical <- 400:2400  # the uniform-1nm part of data.spectral$wlS
n <- length(wl_optical)

op <- par(mfrow = c(2, 2))

plot(wl_optical, res$data.rad$refl[1:n], type = "l", col = "black", lwd = 1.5,
     xlab = "Wavelength (nm)", ylab = "Reflectance", main = "TOC reflectance components")
lines(wl_optical, res$data.rad$rdd[1:n], col = "steelblue")
lines(wl_optical, res$data.rad$rso[1:n], col = "firebrick")
legend("topright", c("refl (apparent)", "rdd (bi-hemispherical)", "rso (bidirectional)"),
       col = c("black", "steelblue", "firebrick"), lty = 1, cex = 0.7)

wlF <- res$data.spectral$wlF
plot(wlF, res$data.rad$LoF_, type = "l", col = "red3", lwd = 1.5,
     xlab = "Wavelength (nm)", ylab = "Fluorescence radiance (mW m-2 nm-1 sr-1)",
     main = sprintf("TOC fluorescence (F685=%.2f, F740=%.2f)", res$data.rad$F685, res$data.rad$F740))

layers <- seq_len(res$data.canopy$nlayers)
plot(layers, res$data.thermal$Tcu, type = "o", pch = 16, col = "darkorange",
     ylim = range(c(res$data.thermal$Tcu, res$data.thermal$Tch, res$data.thermal$Tsu, res$data.thermal$Tsh)),
     xlab = "Canopy layer (1 = top)", ylab = "Temperature (degC)",
     main = sprintf("Converged temperature profile (%d iterations)", res$iter.ebal$counter))
lines(layers, res$data.thermal$Tch, type = "o", pch = 16, col = "steelblue")
abline(h = res$data.thermal$Tsu, col = "brown", lty = 2)
abline(h = res$data.thermal$Tsh, col = "tan", lty = 2)
legend("topright", c("Tcu (sunlit)", "Tch (shaded)", "Tsu (soil, sunlit)", "Tsh (soil, shaded)"),
       col = c("darkorange", "steelblue", "brown", "tan"), lty = c(1, 1, 2, 2), cex = 0.7)

fluxes <- c(Rntot = res$data.fluxes$Rntot, lEtot = res$data.fluxes$lEtot,
            Htot = res$data.fluxes$Htot, Gtot = res$data.fluxes$Gtot)
barplot(fluxes, col = c("darkorange", "steelblue", "firebrick", "grey60"),
        ylab = "W m-2", main = "Canopy energy-balance totals")

par(op)

Four different views of the same one simulation: the reflectance components (Section 2’s refl line among them), the fluorescence spectrum, the converged sunlit/shaded canopy and soil temperatures, and where the absorbed energy actually goes.

4. Many simulations (100-row LUT)

getLUT.SCOPE() builds a LUT the same way ToolsRTM::getLUT() does for PROSAIL – sample every parameter’s own default range from SCOPEinR’s bundled inputs_SCOPE.csv. get.SCOPE() itself is much more expensive per call than foursail() (dominated by the energy-balance nonlinear solve, ~0.1-0.3s per call vs ~2ms, so 100 rows takes on the order of a minute) – see Scripts/R/ForSCOPE/1-getSCOPE.R for the same idea at production scale with diagnostic plots and saved outputs.

inputLUT <- read.table(system.file("input", "inputs_SCOPE.csv", package = "SCOPEinR"),
                        header = TRUE, sep = ",")
N_SAMPLES <- 100
Table.LUT.many <- getLUT.SCOPE(inputLUT = inputLUT, nLUT = N_SAMPLES)

start_time <- Sys.time()
db.sim.many <- get.SCOPE(
  LUT = Table.LUT.many, n.LUT = N_SAMPLES, options.SCOPE = table.with.opts,
  optipar = SCOPEinR::optipar2021.Pro.CX, leaf.model = "fluspect-CX",
  canopy.model = "fourSAIL", get.outputs = "ALL", get.plots = FALSE
)
cat("Ran", N_SAMPLES, "simulations in", round(as.numeric(Sys.time() - start_time), 1), "s\n")
## Ran 100 simulations in 1.9 s
Actot_vals <- sapply(db.sim.many, function(r) r$data.fluxes$Actot)
EoutF_vals <- sapply(db.sim.many, function(r) r$data.rad$EoutF)

op <- par(mfrow = c(1, 2))
plot(Table.LUT.many$Vcmax25[seq_len(N_SAMPLES)], Actot_vals, pch = 19, col = "#2166AC",
     xlab = "Vcmax25 (umol m-2 s-1)", ylab = "Canopy photosynthesis, Actot (umol m-2 s-1)",
     main = "Actot vs Vcmax25")
plot(Table.LUT.many$Cab[seq_len(N_SAMPLES)], EoutF_vals, pch = 19, col = "#B2182B",
     xlab = "Cab", ylab = "Canopy fluorescence flux, EoutF (W/m2)",
     main = "SIF vs Cab")

par(op)

Vcmax25 (maximum carboxylation capacity) is the direct physiological driver of Actot in SCOPE’s Farquhar-type photosynthesis model, so across 100 randomly-drawn LUT rows the relationship comes through clearly. The right panel shows the other half of what makes SCOPEinR more than plain PROSAIL: fluorescence (EoutF) is sensitive to leaf pigments, not just to LAI/structure – substitute Table.LUT.many$Anth or Table.LUT.many$Car for Cab above to see the same fluorescence output against anthocyanin or carotenoid content instead.

5. Trait inversion with machine learning

The reverse problem: given ONLY the sensor-band reflectance a satellite would actually measure, retrieve a trait (LAI here – SCOPE’s classic structural retrieval target) without knowing the ground truth. Convolve each of the 100 simulations’ apparent reflectance (reflapp – SCOPE’s radiance-derived, sensor-observed spectrum, Section 3) onto Sentinel-2A’s bands, then train on 70% of the rows and evaluate on the 30% the model never saw:

library(randomForest)

band_refl <- t(sapply(db.sim.many, function(r) {
  rfl_i <- r$data.rad$reflapp[1:n]
  # reflapp is a radiance ratio (Lo_/incident irradiance), and incident
  # irradiance is genuinely near-zero at a handful of water-vapor/O2
  # absorption wavelengths (~849-850, ~1355-1420, ~1800-1950nm) --
  # numerically unstable there in EVERY simulation, including SCOPEinR's
  # own canonical example row, not just unusual random draws. Real
  # sensors avoid placing bands on these features for the same reason;
  # linearly interpolate over just those narrow gaps before convolving
  # rather than let a handful of bad wavelengths poison an otherwise
  # good spectrum.
  bad <- !is.finite(rfl_i)
  if (any(bad)) rfl_i[bad] <- approx(wl_optical[!bad], rfl_i[!bad], xout = wl_optical[bad])$y
  df_i <- data.frame(wave = wl_optical, rfl = rfl_i)
  get.spectral.convolution.srf(df_i, ToolsRTM::srf.sentinel2a)$RFL
}))
colnames(band_refl) <- paste0("B", seq_len(ncol(band_refl)))
ml_data <- data.frame(band_refl, LAI = Table.LUT.many$LAI)

train_idx <- sample(seq_len(N_SAMPLES), size = round(0.7 * N_SAMPLES))
rf <- randomForest(LAI ~ ., data = ml_data[train_idx, ], ntree = 300)
pred <- predict(rf, ml_data[-train_idx, ])
obs  <- ml_data$LAI[-train_idx]
cat("R2:", cor(pred, obs)^2, " RMSE:", sqrt(mean((pred - obs)^2)), "\n")
## R2: 0.7238131  RMSE: 1.090676
plot(obs, pred, xlab = "Observed LAI", ylab = "Predicted LAI", pch = 19, col = "#2166AC",
     main = "LAI retrieved from Sentinel-2A bands (SCOPE reflapp)")
abline(0, 1, col = "grey40", lty = 2)

Same interface as How-in-R.Rmd’s own step 6 (RandomForest on convolved sensor bands) – what changes here is the input spectrum: reflapp comes out of the full energy-balance/fluorescence-coupled SCOPE run instead of plain PROSAIL, so this trains on what a sensor would see over a canopy whose leaf/soil temperatures were actually solved for, not assumed.

Read the scatter plot, not just R²: the points aren’t scattered evenly around the 1:1 line – low observed LAI (<1.5) is systematically over-predicted, and high observed LAI (>5) is systematically under-predicted, so the point cloud visibly flattens toward the training set’s own mean at both tails. This is real, expected Random Forest behaviour, not a modelling error: an ensemble of trees can only average values it saw in training, so it can never extrapolate past the training range and pulls extreme predictions back toward the bulk of the data – worse here than usual because only 70 rows (0.7 * N_SAMPLES) are available to train on. cor(pred, obs)^2 (a correlation-based R²) is fairly insensitive to exactly this kind of bias, since it only asks whether predicted and observed move together, not whether they agree in absolute terms – which is why the number alone (0.72) undersells how much the retrieval degrades at the LAI extremes visible in the plot above. More/better-distributed training LUT rows (Section 4’s nLUT) or a model less prone to mean-reversion (e.g. a linear/PLSR baseline for comparison) are the two standard fixes, not shown here to keep this tutorial’s scope to “one clean end-to-end example.”

If you compare this R² to the Python version of this section: R’s getLUT.SCOPE() (Section 4) randomizes essentially every free parameter at once – leaf biochemistry, LAI, leaf-angle distribution, and illumination/view geometry (tts/tto/psi) all vary simultaneously across the 100 rows. Python’s version of this tutorial perturbs only Cab/LAI/EWT/Vcmax25 around one fixed baseline, holding geometry and leaf-angle distribution constant. More simultaneously-varying nuisance parameters means more confounding for the model to see through, so R’s LAI retrieval is a genuinely harder task and a lower R² here is expected – not a bug, a real methodological difference between the two tutorials’ sampling designs.

What’s not covered here (and why)

  • Batch LUT runs at course-pipeline scale (N=1000+, chunked/parallel) – see Scripts/R/ForSCOPE/1-getSCOPE.R/2-Explore_outputsSCOPE.R for the full version with diagnostic figures and saved .rds outputs, and SCOPEinR_tutorial.Rmd (Section 7) for get.SCOPE.parallel() and chunked production runs.
  • Convergence diagnostics deep dive – res$iter.ebal already exposes counter/maxEBercu/maxEBerch/maxEBers; see Scripts/R/ForSCOPE/6-validate_ebal_convergence.R for a dedicated analysis of these across many runs.

See CCGCAM/RTM-Suite for the full course materials and the R/Python numerical verification writeup (python/README.md).