GitHub

4SAIL Bidirectional Canopy Reflectance

foursail is the package's implementation of Verhoef's 4SAIL canopy bidirectional reflectance model. It returns a FourSAILResult carrying both the full BRF and the explicit decomposition into direct/direct and diffuse contributions, which makes it easy to diagnose what each pathway contributes at a given wavelength.

using CanopyOptics
using CairoMakie

Geometry, leaf, soil

Geometry and leaf-angle integration are wavelength-independent; build them once with FourSAILGeometry and reuse across wavelengths or moisture/structure sweeps.

LD   = spherical_leaves()
geom = FourSAILGeometry(LD; sza_deg = 30.0, vza_deg = 20.0, raa_deg = 40.0)
FourSAILGeometry{Float64}(0.8660254037844386, 0.9396926207859084, 0.6981317007977318, 0.3792849736778069, 0.5764375927354312, 0.5310341999699336, 0.3877030319495611, 0.0016122646253205672, 0.45439553219196605, 0.12204206054346511, 0.4316938358092173, 0.09934036416071634, 0.6661767358242505, 0.33382326417574953)

Use a small synthetic spectrum (red, NIR, SWIR1) so the table fits in one screen. In practice you would feed in a full PROSPECT-PRO leaf.

λ      = ["red 670 nm", "NIR 800 nm", "SWIR 1.6 µm"]
leaf_R = [0.06, 0.45, 0.30]      # leaf reflectance ρ
leaf_T = [0.02, 0.45, 0.25]      # leaf transmittance τ
soil_R = [0.10, 0.18, 0.30]
lai    = 3.0

result = foursail(leaf_R, leaf_T, soil_R, geom, lai; hotspot = 0.05)
FourSAILResult{Vector{Float64}}([0.024374017638279313, 0.4772551932668812, 0.202720598887185], [0.020333686336998856, 0.3863433353220289, 0.16184893701656813], [0.019804949372136785, 0.37199404469972364, 0.15601155814642953], [0.0270877730285333, 0.37081375916765136, 0.17445637283238835], [0.18275953097652572, 0.42245984858565505, 0.2578972036277135], [0.02231101214773805, 0.16779585112276874, 0.11178669074605693], [0.026382125305708812, 0.17512385480711612, 0.12400003021996922], [0.0007056477228244871, 0.1956899043605352, 0.05045634261241914])

result.rsot is the full bidirectional reflectance factor — what a real instrument would observe at the chosen sun/view geometry.

println("Full BRF (rsot):")
for (name, v) in zip(λ, result.rsot)
    println("  ", rpad(name, 12), round(v, digits = 4))
end
Full BRF (rsot):
  red 670 nm  0.0271
  NIR 800 nm  0.3708
  SWIR 1.6 µm 0.1745

The rsot = rsost + rsodt decomposition

4SAIL splits the bidirectional reflectance into two physically distinct branches:

  • rsost — the direct/direct branch: single scattering from leaves plus directly transmitted soil reflectance reaching the sensor without any diffuse scattering. Treat it as a lower bound on what the sensor sees if you were to ignore multiple scattering.
  • rsodt — the diffuse / multiple-scattering remainder that has to be added to rsost to recover the full BRF.

By construction rsot ≈ rsost + rsodt. In the NIR band, where leaves are highly scattering (ρ + τ ≈ 0.9), most of the signal arrives through diffuse multiple scattering — rsost alone severely underestimates the canopy reflectance.

println("\nrsost (direct/direct) vs rsodt (diffuse) vs rsot (full):")
for (name, rost, rodt, rot) in zip(λ, result.rsost, result.rsodt, result.rsot)
    println("  ", rpad(name, 12),
            "  rsost=", round(rost, digits = 4),
            "  rsodt=", round(rodt, digits = 4),
            "  rsot=",  round(rot,  digits = 4))
end

rsost (direct/direct) vs rsodt (diffuse) vs rsot (full):
  red 670 nm    rsost=0.0264  rsodt=0.0007  rsot=0.0271
  NIR 800 nm    rsost=0.1751  rsodt=0.1957  rsot=0.3708
  SWIR 1.6 µm   rsost=0.124  rsodt=0.0505  rsot=0.1745

Confirm the decomposition holds to machine precision.

maximum(abs, result.rsot .- (result.rsost .+ result.rsodt))
0.0

Wavelength sweep — visualizing the decomposition

Sweep a synthetic ρ(λ), τ(λ) that interpolates from a red-band absorption regime through the NIR plateau. The gap between rsost (the "if everything were single-scattering" curve) and rsot is the diffuse contribution rsodt.

λs    = range(0.4, 2.5, length = 80)
ρ_λ   = @. 0.04 + 0.42 / (1 + exp(-(λs - 0.72) * 25))  # red dip → NIR plateau
τ_λ   = 0.9 .* ρ_λ                                     # τ tracks ρ in this toy example
soil_λ = fill(0.20, length(λs))

sweep = foursail(ρ_λ, τ_λ, soil_λ, geom, lai; hotspot = 0.05)

fig = Figure(size = (640, 360))
ax = Axis(fig[1, 1]; xlabel = "wavelength (µm)", ylabel = "reflectance",
          title = "4SAIL BRF and its decomposition  (LAI = 3, sza/vza/raa = 30/20/40°)")
lines!(ax, λs, sweep.rsost; label = "rsost (direct/direct)", linestyle = :dash)
lines!(ax, λs, sweep.rsodt; label = "rsodt (diffuse)",       linestyle = :dot)
lines!(ax, λs, sweep.rsot;  label = "rsot  (full BRF)",      linewidth = 2)
axislegend(ax; position = :rt)
fig

Batched geometries — single-pass directional response

For BRDF studies pass a FourSAILGeometrySet instead of a single geometry: 4SAIL evaluates all (wavelength × angle) combinations in one batched pass (CPU or GPU via KernelAbstractions, depending on the array backend).

vzas  = collect(0.0:5.0:60.0)
raas  = fill(40.0, length(vzas))
geoms = FourSAILGeometrySet(LD; sza_deg = 30.0, vza_deg = vzas, raa_deg = raas)
brdf  = foursail(leaf_R, leaf_T, soil_R, geoms, lai; hotspot = 0.05)
FourSAILResult{Matrix{Float64}}([0.024374017638279313 0.024374017638279313 … 0.024374017638279313 0.024374017638279313; 0.4772551932668812 0.4772551932668812 … 0.4772551932668812 0.4772551932668812; 0.202720598887185 0.202720598887185 … 0.202720598887185 0.202720598887185], [0.020333686336998856 0.020333686336998856 … 0.020333686336998856 0.020333686336998856; 0.3863433353220289 0.3863433353220289 … 0.3863433353220289 0.3863433353220289; 0.16184893701656813 0.16184893701656813 … 0.16184893701656813 0.16184893701656813], [0.019428441720794717 0.01944894130448475 … 0.023306738847595158 0.024379716587449217; 0.3613804184831989 0.361966974061223 … 0.4561523039346512 0.4773629128017519; 0.1517761361662266 0.15200846654989766 … 0.19250592194276608 0.20277414579266625], [0.024374842645832003 0.02510740688459012 … 0.02676299052380305 0.026620872945358445; 0.34697679748300897 0.3520448592885843 … 0.4089608864430321 0.4154833347602742; 0.16059198492652932 0.16398656321165428 … 0.1828566334387175 0.18432291691840796], [0.18275953097652572 0.18275953097652572 … 0.18275953097652572 0.18275953097652572; 0.42245984858565505 0.42245984858565505 … 0.42245984858565505 0.42245984858565505; 0.2578972036277135 0.2578972036277135 … 0.2578972036277135 0.2578972036277135], [0.019407145346435364 0.020119688150113668 … 0.024606696164362736 0.024904437042259742; 0.14721641497624208 0.15204578146867767 … 0.18820775320038763 0.19211202827965423; 0.09786713917116524 0.10117250092198092 … 0.12486224680564721 0.1271865604426518], [0.023677001142877102 0.024409130244416037 … 0.026001187306798016 0.02584679767729551; 0.1549021554098372 0.15976677723842192 … 0.19071783725677113 0.1938082774227186; 0.11067670656049045 0.11404082720488802 … 0.12904572023295305 0.1300136423477591], [0.0006978415029549012 0.0006982766401740805 … 0.000761803217005035 0.0007740752680629343; 0.19207464207317176 0.19227808205016236 … 0.21824304918626092 0.2216750573375556; 0.04991527836603889 0.049945736006766245 … 0.05381091320576445 0.054309274570648866])

brdf.rsot is now a (nwavelength × nangle) matrix.

size(brdf.rsot)
(3, 13)

This page was generated using Literate.jl.