Examples

Every snippet on this page has actually been run against the real toolsrtm/scopeinpython packages – not hand-written pseudocode. These mirror the R tutorial series (ToolsRTM, SCOPEinR) topic-for-topic where a Python port exists; see each package’s own README.md for the full R-tutorial-to-Python-module bridge table, including the gaps called out at the bottom of this page.

Leaf + canopy (toolsrtm)

Mirrors R Tutorials 01-02 (prospect_d = leaf optics, foursail = canopy BRDF).

import numpy as np
from toolsrtm import prospect_d, foursail

leaf = prospect_d(N=1.5, Cab=40, Car=8, Anth=1, Cbrown=0, EWT=0.01, LMA=0.009, alpha=40)
print(leaf.lambda_[:3], leaf.refl[:3], leaf.tran[:3])

inputLUT = dict(
    N=1.5, Cab=40, Car=8, Anth=1, Cbrown=0, EWT=0.01, LMA=0.009, alpha=40,
    Prot=0.002, CBC=0.007,          # only used if leaf_model='PROSPECT-PRO'
    LIDFa=-0.35, LIDFb=-0.15, TypeLidf=1,
    LAI=3, hspot=0.01, tts=30, tto=0, psi=0,
)
rsoil = np.full(2101, 0.15)
sail = foursail(inputLUT, rsoil, leaf_model="PROSPECT-D", spectrum_all=True)
print("TOC bidirectional reflectance factor (rsot) at 550 nm:", sail.rsot[550 - 400])

Soil (BSM) + optical canopy BRDF (scopeinpython)

Mirrors SCOPEinR Tutorials 01-02 (soil model + canopy optics via rtmo).

import numpy as np
from toolsrtm import prospect_d
from toolsrtm.canopy import dladgen
from scopeinpython import SoilParams, WettingParams, get_bsm, CanopyStructure, get_spectra_scope, run_rtmo

spectral = get_spectra_scope()
rsoil = get_bsm(SoilParams(BSMBrightness=0.5, BSMlat=25, BSMlon=45),
                 WettingParams(SMp=15, SMC=25, film=0.015))

leaf = prospect_d(N=1.5, Cab=40, Car=8, Anth=1, Cbrown=0, EWT=0.01, LMA=0.009, alpha=40)
refl_leaf, tran_leaf = leaf.refl[:2001], leaf.tran[:2001]

lidf = dladgen(-0.35, -0.15).lidf
canopy = CanopyStructure(LAI=3, lidf=lidf, hot=0.1 / 2.0)

# Esun_/Esky_ must be supplied by the caller -- see python/README.md
result = run_rtmo(
    spectral=spectral, leaf_refl=refl_leaf, leaf_tran=tran_leaf,
    rho_thermal=0.01, tau_thermal=0.01, rsoil=rsoil, canopy=canopy,
    tts=30, tto=0, psi=0, Esun_=Esun_, Esky_=Esky_,
)
print("TOC reflectance at 550 nm:", result.refl[550 - 400])

Sensor convolution + vegetation indices (toolsrtm)

Mirrors R Tutorials 07-09: convolve a simulated hyperspectral canopy spectrum onto real Sentinel-2A band spectral response functions, then compute vegetation indices from the convolved bands.

import numpy as np
from toolsrtm import prospect_d, foursail
from toolsrtm.srf import srf_sentinel2a, spectral_convolution_srf
from toolsrtm.indices import get_indices

inputLUT = dict(
    N=1.5, Cab=40, Car=8, Anth=1, Cbrown=0, EWT=0.01, LMA=0.009, alpha=40,
    LIDFa=-0.35, LIDFb=-0.15, TypeLidf=1,
    LAI=3, hspot=0.01, tts=30, tto=0, psi=0,
)
rsoil = np.full(2101, 0.15)
sail = foursail(inputLUT, rsoil, leaf_model="PROSPECT-D", spectrum_all=True)
wave = np.arange(400, 2501)

s2a = srf_sentinel2a()
conv = spectral_convolution_srf(wave, sail.rsot, s2a)
print("Sentinel-2A bands:", s2a.band_names)
print("Convolved TOC reflectance:", np.round(conv.rfl, 4))

indices = get_indices(conv.wl, conv.rfl, spectral_domain="VNIR")
for name in ("NDVI", "MSAVI", "REP"):
    print(name, "=", round(float(indices[name][0]), 4))

Machine-learning trait inversion (toolsrtm)

Mirrors R Tutorials 11-12: build a small LUT, extract Sentinel-2-like bands, and invert Cab with a PLSR model (get_inversion dispatches to 12 algorithms in total – see toolsrtm.inversion.ALGORITHMS).

import numpy as np
import pandas as pd
from toolsrtm import prospect_d, foursail
from toolsrtm.inversion import get_inversion

rng = np.random.default_rng(1)
rows = []
for _ in range(200):
    Cab, LAI = rng.uniform(10, 80), rng.uniform(0.5, 6)
    inputLUT = dict(
        N=1.5, Cab=Cab, Car=8, Anth=1, Cbrown=0, EWT=0.01, LMA=0.009, alpha=40,
        LIDFa=-0.35, LIDFb=-0.15, TypeLidf=1,
        LAI=LAI, hspot=0.01, tts=30, tto=0, psi=0,
    )
    sail = foursail(inputLUT, np.full(2101, 0.15), leaf_model="PROSPECT-D", spectrum_all=True)
    row = {"Cab": Cab, "LAI": LAI}
    for wl in (490, 560, 665, 705, 740, 783, 842, 865, 1610, 2190):
        row[f"R{wl}"] = sail.rsot[wl - 400]
    rows.append(row)

df = pd.DataFrame(rows)
band_cols = [c for c in df.columns if c.startswith("R")]

result = get_inversion(df, dep_var="Cab", inputs=band_cols, algorithm="PLSR", n_samples=200, seed=1)
print("Test R2:", round(result.statistics["test"]["r2"], 3))
print("Test RMSE:", round(result.statistics["test"]["rmse"], 3))

MARMIT soil moisture model (toolsrtm)

Mirrors R Tutorial 16: build a wetted soil spectrum from a dry reference and couple it into a canopy simulation.

from toolsrtm.marmit import get_marmit_rsoil
from toolsrtm import prospect_d, foursail

soil = get_marmit_rsoil(soil_id=3, L=0.05, eps=0.4, version="marmit1")
print("SMC (soil moisture content):", round(float(soil.smc), 4))
for wl in (550, 850, 1600):
    i = wl - 400
    print(f"{wl}nm: dry={soil.rsoil_dry[i]:.4f}  wet={soil.rsoil_wet[i]:.4f}")

inputLUT = dict(
    N=1.5, Cab=40, Car=8, Anth=1, Cbrown=0, EWT=0.01, LMA=0.009, alpha=40,
    LIDFa=-0.35, LIDFb=-0.15, TypeLidf=1,
    LAI=1.5, hspot=0.01, tts=30, tto=0, psi=0,
)
sail = foursail(inputLUT, soil.rsoil_wet, leaf_model="PROSPECT-D", spectrum_all=True)
print("Canopy TOC reflectance at 850nm with wet MARMIT soil:", round(float(sail.rsot[850 - 400]), 4))

Full SCOPE run: energy balance + fluorescence (scopeinpython)

Mirrors SCOPEinR Tutorials 03-04: one full get_scope() call – optics, energy balance, photosynthesis and fluorescence together – against the same bundled example LUT row SCOPEinR’s own test suite and R vignettes use (SCOPEinR/inst/input/LUT_input.csv).

import csv
from pathlib import Path
from scopeinpython import ScopeOptions, get_scope

with open("SCOPEinR/inst/input/LUT_input.csv", newline="") as f:
    row = next(csv.DictReader(f))

res = get_scope(row, options=ScopeOptions(k_maxit=100, maxEBer=1.0))
print("Canopy layers:", res.nlayers)
print("TOC reflectance at 550/700/850 nm:",
      round(float(res.rtmo.refl[550 - 400]), 4),
      round(float(res.rtmo.refl[700 - 400]), 4),
      round(float(res.rtmo.refl[850 - 400]), 4))
print("Net radiation, total (Rntot):", round(float(res.ebal.Rntot), 2), "W/m2")
print("Total photosynthesis (Actot):", round(float(res.ebal.Actot), 2), "umol CO2/m2/s")
if res.rtmf is not None:
    print("Emitted fluorescence (EoutF):", round(float(res.rtmf.EoutF), 4), "W/m2/sr")

Full simulate -> indices -> ML-invert pipeline

A complete, actually-executed pipeline (100 samples, spectral indices, scikit-learn trait inversion, real R² 0.7-0.9) lives outside this package as plain scripts + a Jupyter notebook, kept separate from the R scripts:

Scripts/Python/ForPROSAIL_fourSAIL/ — 1_simulate_lut.py → 2_spectral_indices.py → 3_inversion_ml.py, plus pipeline.ipynb. See Scripts/Python/README.md for how to run it.

What isn’t ported yet

Two tutorial-level gaps not covered by any example above, on top of the finer-grained implementation gaps in Known limitations:

Sobol/Johnson sensitivity analysis (toolsrtm, R Tutorial 10’s get.sobol.indices()/get.spectral.sensitivity()) and LUT correlation-distribution helpers (get_distributionLUT()/getCor(), R Tutorial 05) have no Python module yet.
The real-Sentinel-2/STAC capstone (SCOPEinR Tutorial 11 – Actot retrieved from a real satellite time series and mapped spatially; ToolsRTM Tutorial 18 – a spatial index map from a real STAC-retrieved image) has no scopeinpython/toolsrtm equivalent yet, even though toolsrtm.satellite already has the STAC retrieval machinery a Python port of it would reuse.

See each package’s own README.md for the full R-tutorial-to-Python-module bridge table.