← 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 2.2 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.

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