From 222d2b32a9e40b89fc4d6b2df67c6e87b61e1c00 Mon Sep 17 00:00:00 2001 From: cfranken Date: Thu, 14 May 2026 15:14:16 -0700 Subject: [PATCH] docs(foursail): add hotspot and polar BRDF tutorials foursail_hotspot.jl walks through the principal-plane BRF at 685 and 800 nm with PROSPECT-PRO leaves, contrasting rsot with and without the Kuusk hotspot and overlaying the embedded rsost branch. foursail_brdf.jl uses FourSAILGeometrySet to evaluate the full (VZA, VAZ) half-sphere in one batched call and visualizes the angular NDVI effect against nadir, including the hotspot in 4SAIL's raa = 0 convention. Co-Authored-By: Claude Opus 4.7 (1M context) --- docs/make.jl | 3 +- docs/src/pages/tutorials/foursail_brdf.jl | 158 +++++++++++++++++ docs/src/pages/tutorials/foursail_hotspot.jl | 169 +++++++++++++++++++ 3 files changed, 329 insertions(+), 1 deletion(-) create mode 100644 docs/src/pages/tutorials/foursail_brdf.jl create mode 100644 docs/src/pages/tutorials/foursail_hotspot.jl diff --git a/docs/make.jl b/docs/make.jl index e24f0c3..058abf5 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -6,7 +6,8 @@ using Literate, UnitfulEquivalences, Distributions @info "Building docs with CanopyOptics source" path=pathof(CanopyOptics) function build() - tutorials = ["bilambertian.jl", "specular.jl", "foursail.jl", "dielectric.jl"] + tutorials = ["bilambertian.jl", "specular.jl", "foursail.jl", + "foursail_hotspot.jl", "foursail_brdf.jl", "dielectric.jl"] tutorials_paths = [joinpath(@__DIR__, "src", "pages", "tutorials", tutorial) for tutorial in tutorials] for tutorial in tutorials_paths diff --git a/docs/src/pages/tutorials/foursail_brdf.jl b/docs/src/pages/tutorials/foursail_brdf.jl new file mode 100644 index 0000000..3f4f683 --- /dev/null +++ b/docs/src/pages/tutorials/foursail_brdf.jl @@ -0,0 +1,158 @@ +# # 4SAIL Polar BRDF and the Angular NDVI Effect + +# Where [the principal-plane tutorial](foursail_hotspot.md) sweeps a +# single line through view-angle space, this one runs 4SAIL over the +# full `(VZA, VAZ)` half-sphere in a single batched call and uses the +# result to look at the *angular* dependence of NDVI. +# +# Two ideas are worth highlighting: +# +# * [`FourSAILGeometrySet`](@ref) lets you evaluate 4SAIL for an entire +# directional grid in one go — no Julia-level loop over view angles. +# * NDVI is not a single number for a canopy: even a vegetation-only +# scene with a fixed leaf spectrum shows several percent of angular +# spread, with a peak on the hotspot side. +using CanopyOptics +using CairoMakie + +# ## Scene +const sza_deg = 30.0 +const lai = 4.0 +const lambda_nm = [685.0, 800.0] +const soil_albedo = 0.0 +const hotspot_h = 0.05 + +# Dense (VZA, VAZ) grid. `FourSAILGeometrySet` consumes flattened +# (VZA, VAZ) vectors; we reshape the results back to a 2D grid after. +vzas_deg = collect(0.0:5.0:80.0) +vazs_deg = collect(0.0:5.0:355.0) +n_vza, n_vaz = length(vzas_deg), length(vazs_deg) + +vza_all = repeat(vzas_deg; inner = n_vaz) +raa_all = repeat(vazs_deg; outer = n_vza) # 4SAIL's raa = azimuth difference + +# ## Leaf optics from PROSPECT-PRO +leaf = CanopyOptics.LeafProspectProProperties{Float64}( + N = 1.5, Ccab = 40.0, Ccar = 8.0, + Canth = 0.0, Cbrown = 0.0, + Cw = 0.012, Cm = 0.009, + Cprot = 0.0, Ccbc = 0.0, +) + +opti = CanopyOptics.createLeafOpticalStruct(400.0:1.0:2500.0) +T_grid, R_grid = CanopyOptics.prospect(leaf, opti) +λ_grid = [Float64(v.val) for v in opti.λ] +grid_idx = [argmin(abs.(λ_grid .- λ)) for λ in lambda_nm] +R_leaf = R_grid[grid_idx] +T_leaf = T_grid[grid_idx] + +# ## One batched 4SAIL call for the entire half-sphere +LD = CanopyOptics.planophile_leaves2(Float64) +geoms = FourSAILGeometrySet(LD; + sza_deg = sza_deg, + vza_deg = vza_all, + raa_deg = raa_all, + quadrature = CanopyOptics.CanopyQuadrature(n_leaf = 64, n_azimuth = 16), +) + +result = foursail(R_leaf, T_leaf, soil_albedo, geoms, lai; + hotspot = hotspot_h, mode = :full) + +# `result.rsot` is (n_wavelength × n_directions). Reshape each +# wavelength row back to (n_vza × n_vaz): +brf = Array{Float64}(undef, length(lambda_nm), n_vza, n_vaz) +for iλ in eachindex(lambda_nm) + brf[iλ, :, :] .= permutedims(reshape(result.rsot[iλ, :], n_vaz, n_vza)) +end + +R_red = @view brf[1, :, :] +R_nir = @view brf[2, :, :] +NDVI = (R_nir .- R_red) ./ (R_nir .+ R_red) + +println("685 nm BRF extrema: ", extrema(R_red)) +println("800 nm BRF extrema: ", extrema(R_nir)) +println("NDVI extrema: ", extrema(NDVI)) + +# Anchor the angular NDVI effect to nadir so we can see *how much* it +# varies across the view hemisphere relative to the most common +# space-based reference geometry. +nadir_ndvi = NDVI[1, 1] +ΔNDVI = NDVI .- nadir_ndvi +println("NDVI at nadir: ", nadir_ndvi, + " ΔNDVI extrema (relative to nadir): ", extrema(ΔNDVI)) + +# ## Polar BRDF plots +# +# Closing the azimuth loop (`hcat(..., grid[:, 1])`) avoids the visible +# seam at `vaz = 0/360°` in the polar surface plots. The red star marks +# the hotspot direction at `(vza = sza, vaz = 0°)`, which in 4SAIL's +# convention is the same azimuth as the sun. +θ_rad = deg2rad.(vcat(vazs_deg, 360.0)) +close_az(grid) = hcat(grid, grid[:, 1]) + +fig = Figure(size = (1500, 700)) + +for (iλ, λ) in pairs(lambda_nm) + grid = brf[iλ, :, :] + ax = PolarAxis(fig[1, 2(iλ-1) + 1]; + title = "$(Int(round(λ))) nm BRF (rsot)", + theta_0 = π / 2, + direction = -1, + rticks = 0:20:80, + rticklabelsize = 10, + thetaticks = (deg2rad.(0:45:315), string.(0:45:315) .* "°")) + hm = surface!(ax, θ_rad, vzas_deg, close_az(grid)'; + colormap = :viridis, shading = NoShading) + scatter!(ax, [deg2rad(0.0)], [sza_deg]; + color = :red, markersize = 12, marker = :star5, + strokecolor = :white, strokewidth = 1.0) + Colorbar(fig[1, 2(iλ-1) + 2], hm; label = "BRF") +end + +Δlim = maximum(abs, ΔNDVI) +ax_ndvi = PolarAxis(fig[2, 1]; + title = "NDVI", + theta_0 = π / 2, direction = -1, + rticks = 0:20:80, rticklabelsize = 10, + thetaticks = (deg2rad.(0:45:315), string.(0:45:315) .* "°")) +hm_ndvi = surface!(ax_ndvi, θ_rad, vzas_deg, close_az(NDVI)'; + colormap = :viridis, shading = NoShading) +Colorbar(fig[2, 2], hm_ndvi; label = "NDVI") + +ax_delta = PolarAxis(fig[2, 3]; + title = "NDVI − NDVI(nadir)", + theta_0 = π / 2, direction = -1, + rticks = 0:20:80, rticklabelsize = 10, + thetaticks = (deg2rad.(0:45:315), string.(0:45:315) .* "°")) +hm_delta = surface!(ax_delta, θ_rad, vzas_deg, close_az(ΔNDVI)'; + colormap = :balance, colorrange = (-Δlim, Δlim), + shading = NoShading) +Colorbar(fig[2, 4], hm_delta; label = "ΔNDVI") + +for ax in (ax_ndvi, ax_delta) + scatter!(ax, [deg2rad(0.0)], [sza_deg]; + color = :red, markersize = 12, marker = :star5, + strokecolor = :white, strokewidth = 1.0) +end + +Label(fig[0, 1:4], + "4SAIL polar BRDF and NDVI angular effect " * + "(SZA = $(Int(sza_deg))°, LAI = $(lai), planophile LAD, soil albedo = $(soil_albedo), h = $(hotspot_h))"; + fontsize = 16) +fig + +# ## Reading the figure +# +# * The 685 nm BRDF is dim and dominated by the hotspot peak: most of +# the red signal *is* the directly back-scattered branch, because +# chlorophyll absorbs almost everything that scatters into the diffuse +# field. +# * The 800 nm BRDF is much brighter and noticeably flatter, with a +# broad bowl shape — multiple scattering smears the angular contrast +# out across the upper hemisphere. +# * The two effects combine in the NDVI panel: NDVI dips on the hotspot +# side (where the red BRF is locally lifted) and rises on the +# forward-scatter side. Even with a fixed leaf spectrum and a +# horizontally homogeneous canopy this is a several-percent angular +# variation — large enough to matter for cross-instrument comparison +# if view geometries differ. diff --git a/docs/src/pages/tutorials/foursail_hotspot.jl b/docs/src/pages/tutorials/foursail_hotspot.jl new file mode 100644 index 0000000..04ab618 --- /dev/null +++ b/docs/src/pages/tutorials/foursail_hotspot.jl @@ -0,0 +1,169 @@ +# # 4SAIL Hotspot in the Solar Principal Plane + +# This tutorial walks through a red/NIR principal-plane BRF diagnostic +# with [`foursail`](@ref), using a realistic PROSPECT-PRO leaf and the +# [`FourSAILGeometrySet`](@ref) batched entry point. It shows three +# things in one figure: +# +# * the full BRF `rsot` with hotspot off (`h = 0`) +# * the full BRF `rsot` with a Kuusk-style hotspot turned on +# * the direct/direct branch `rsost` from the same hotspot solve +# +# The shape `rsot = rsost + rsodt` is what makes 4SAIL convenient as a +# diagnostic tool — the embedded `rsost` is the path that carries the +# hotspot, so you can read off the hotspot increment without rerunning +# the model in a reduced mode. +using CanopyOptics +using CairoMakie + +# ## Geometry — signed VZA in the principal plane +# +# The "principal plane" is the vertical plane containing both the sun +# and the observer. We sweep view zenith from −80° to +80° and split +# the sign by relative azimuth: +# +# * `signed_vza < 0` → `raa = 0°`: backscatter side, same azimuth as +# the sun, the side that carries the hotspot. +# * `signed_vza > 0` → `raa = 180°`: forward-scatter side. +# +# (4SAIL's `raa_deg` argument is the *azimuth difference* between sun +# and view directions, so RAA = 0° puts the observer in the same +# half-plane as the sun.) +const sza_deg = 30.0 +const lai = 4.0 +const hotspot_h = 0.05 +const lambda_nm = [685.0, 800.0] + +signed_vza = collect(-80.0:2.0:80.0) +vza_abs = abs.(signed_vza) +raa_deg = ifelse.(signed_vza .< 0, 0.0, 180.0) + +LD = CanopyOptics.planophile_leaves2(Float64) +geoms = FourSAILGeometrySet(LD; + sza_deg = sza_deg, + vza_deg = vza_abs, + raa_deg = raa_deg, + quadrature = CanopyOptics.CanopyQuadrature(n_leaf = 64, n_azimuth = 16), +) + +# ## Leaf optics from PROSPECT-PRO +# +# Evaluate PROSPECT once on a 1 nm grid in the 400–2500 nm range, then +# pull the two requested wavelengths out by nearest-neighbor index. +leaf = CanopyOptics.LeafProspectProProperties{Float64}( + N = 1.5, Ccab = 40.0, Ccar = 8.0, + Canth = 0.0, Cbrown = 0.0, + Cw = 0.012, Cm = 0.009, + Cprot = 0.0, Ccbc = 0.0, +) + +opti = CanopyOptics.createLeafOpticalStruct(400.0:1.0:2500.0) +T_grid, R_grid = CanopyOptics.prospect(leaf, opti) +λ_grid = [Float64(v.val) for v in opti.λ] +grid_idx = [argmin(abs.(λ_grid .- λ)) for λ in lambda_nm] +R_leaf = R_grid[grid_idx] +T_leaf = T_grid[grid_idx] + +println("PROSPECT leaf optics at the two probe wavelengths:") +for (λ, ρ, τ) in zip(lambda_nm, R_leaf, T_leaf) + println(" λ = ", λ, " nm ρ = ", round(ρ, digits = 4), + " τ = ", round(τ, digits = 4)) +end + +# ## Two 4SAIL solves — hotspot off, hotspot on +# +# `foursail` returns a [`FourSAILResult`](@ref). For the batched +# `FourSAILGeometrySet` path each field is a (n_wavelength × n_angle) +# matrix. +soil_albedo = 0.0 + +sail_nohot = foursail(R_leaf, T_leaf, soil_albedo, geoms, lai; + hotspot = 0.0, mode = :full) +sail_hot = foursail(R_leaf, T_leaf, soil_albedo, geoms, lai; + hotspot = hotspot_h, mode = :full) + +size(sail_hot.rsot) + +# Quick numeric readout at three reference angles: the hotspot +# (`signed_vza = -sza`), nadir, and the symmetric forward-scatter angle. +hotspot_idx = findfirst(==(-sza_deg), signed_vza) +nadir_idx = findfirst(==(0.0), signed_vza) +forward_idx = findfirst(==(sza_deg), signed_vza) + +for (iλ, λ) in pairs(lambda_nm) + println("λ = ", λ, " nm") + for (label, idx) in (("hotspot (-sza)", hotspot_idx), + ("nadir", nadir_idx), + ("forward (+sza)", forward_idx)) + idx === nothing && continue + println(" ", rpad(label, 16), + " rsot(h=0) = ", round(sail_nohot.rsot[iλ, idx], digits = 4), + " rsot(h=$(hotspot_h)) = ", round(sail_hot.rsot[iλ, idx], digits = 4), + " rsost = ", round(sail_hot.rsost[iλ, idx], digits = 4)) + end +end + +# ## Plot — full BRF and the hotspot increment +# +# Top row: BRF curves themselves (with and without hotspot, plus the +# direct/direct branch `rsost`). Bottom row: the hotspot increment +# `rsot(h) − rsot(0)` is concentrated near `vza = −sza`, exactly where +# the camera looks back along the sun direction. +fig = Figure(size = (1100, 700)) +colors = (sail_nohot = :dodgerblue3, + sail_hot = :seagreen4, + sail_ss = :darkorange3, + increment = :purple4) + +for (iλ, λ) in pairs(lambda_nm) + ax = Axis(fig[iλ, 1]; + xlabel = iλ == length(lambda_nm) ? + "signed VZA in solar principal plane (deg)" : "", + ylabel = "BRF", + title = "$(Int(round(λ))) nm — principal-plane BRF") + lines!(ax, signed_vza, sail_nohot.rsot[iλ, :]; + color = colors.sail_nohot, linewidth = 2.2, linestyle = :dash, + label = "rsot, h = 0") + lines!(ax, signed_vza, sail_hot.rsot[iλ, :]; + color = colors.sail_hot, linewidth = 2.5, + label = "rsot, h = $(hotspot_h)") + lines!(ax, signed_vza, sail_hot.rsost[iλ, :]; + color = colors.sail_ss, linewidth = 2.0, linestyle = :dot, + label = "rsost (direct/direct), h = $(hotspot_h)") + vlines!(ax, [-sza_deg]; color = :gray30, linestyle = :dashdot, linewidth = 1.2, + label = "hotspot direction") + vlines!(ax, [0.0]; color = :gray60, linestyle = :dot, linewidth = 1.0) + axislegend(ax; position = :lt, framevisible = false) + + axd = Axis(fig[iλ, 2]; + xlabel = iλ == length(lambda_nm) ? + "signed VZA in solar principal plane (deg)" : "", + ylabel = "Δ BRF", + title = "$(Int(round(λ))) nm — hotspot increment") + lines!(axd, signed_vza, sail_hot.rsot[iλ, :] .- sail_nohot.rsot[iλ, :]; + color = colors.increment, linewidth = 2.4, + label = "rsot(h) − rsot(0)") + hlines!(axd, [0.0]; color = :gray70, linewidth = 1.0) + vlines!(axd, [-sza_deg]; color = :gray30, linestyle = :dashdot, linewidth = 1.2) + axislegend(axd; position = :rt, framevisible = false) +end + +Label(fig[0, 1:2], + "4SAIL principal-plane BRF and hotspot increment " * + "(SZA = $(Int(sza_deg))°, LAI = $(lai), planophile LAD, soil albedo = $(soil_albedo))"; + fontsize = 16) +fig + +# ## Reading the figure +# +# * Without the hotspot (`h = 0`) the BRF is smooth across the principal +# plane. The slight back–forward asymmetry is just the LAD-weighted +# single-scattering shape. +# * Turning the hotspot on lifts the curve sharply near `vza = -sza`. +# The increment is essentially zero on the forward-scatter side, so +# the multiple-scattering remainder `rsodt` is hardly perturbed — the +# hotspot lives in the `rsost` branch. +# * `rsost` itself sits well below the full BRF in the NIR because most +# of the 800 nm signal arrives through diffuse multiple scattering, +# not single scattering. In the red band leaves absorb most of what +# they scatter once, so `rsost` and the full `rsot` are much closer.