Optical properties of leaves¶
J Gómez-Dans (NCEO & UCL)¶
In order to understand the reflectance signal acquired by a passive sensor that captures the sunlight scattered by a canopy, we need to understand how the different pigments that contribute to reflectance and transmittance of leaves, as well as understanding the geometrical disposition of leaves within the canopy. We will first look at an individual leaf, and try to come up with a simple model of how reflectance and transmittance can be calculated.
We will focus on the widely used PROSPECT model by Jacquemoud et al, 1990. Although other models do exist, PROSPECT is straightforward and works reasonably well.
Some first concepts¶
As a flux of photons hits an object, the photons can either be absorbed by the object ($A$), scattered away ($R$), or be transmitted ($T$) by the object. Conservation laws imply that the sum of absorptance, reflectance and transmittance is one:
$$
A + R + T = 1.
$$
Note that the value of these three magnitudes varies with wavelength ($\lambda$). Understanding the processes that give rise to these changes over the spectrum allows us to infer something about the leaf biochemical and structural composition.
We further consider that reflectance has a surface ($R_s$) component and a component due to the multiple scattering that takes place within the leaf, $R_d$. The first component is due to the air-leaf interface and the waxy leaf surface, and is fairly important. This surface reflectance is controlled thus by the surface properties and leaf surface microtopography. It has been observed that this component is far from Lambertian (isotropic), exhibiting a strong directional effect.
$R_d$ is mostly a consequence of air-cell wall interfaces. As a photon beam enters the leaf, it will traverse it suffering reflectance every time it encounters a discontinuity, following a unique pattern within the leaf, and either being reflected "upwards" through the leaf surface or "downwards" through the opposite surface of the leaf. These interactions are assumed to provide a reflected and transmitted fluxes that are broadly non-directional (Lambertian). The combination of $R_s$ and $R_d$ result in the observed directional and Lambertian properties of the leaf BRDF (bi-hemispherical reflectance distribution function)
The diffuse term thus is the result of photons traversing the leaf and interacting with different components of the mesophyll, and thus opens the possibiity to learn about the internal composition and structure of leaves by looking at how light is scattered by them.
We define the optical domain to lie between 400 and 2500$nm$. The following plot shows the spectral reflectance and transmittance from a clover leaf from the LOPEX'93 database:
In the visible region (400-700nm), there is a strong absorption of radiation by green leaves (note the peak in the green region, around 550nm). If chlorophyll concentration is low, absorption in this region is lower, and therefore, reflectance increases. In the NIR plateau (700-1100nm), absorption is controlled by cellulose and lignine, but also by the multiple scattering within the leaf, and thus, to the internal structure and interfaces of refraction within leaf layers. Towards the SWIR region, we see relatively strong absorption due to water, with some very important absorption peaks around 1450, 1950 and 2500 nm.
A simple plate model of leaf optical properties¶
A simple approximation to a monocot leaf can be a plate model, where we see the leaf as a layer with a particular refraction index $n_2$. The model can be depicted as
An incoming beam that enters the layer is partially reflected, partially transmitted and partially absorbed. The total reflectance can be calculated following the work of Airy, and after calculating an infinite summation yields
$$ R = r_{12}+\frac{t_{12}t_{21}r_{21}\tau^2}{1-r_{21}^2\tau^2}, $$where $r_{12}$ and $t_{12}$ are the average reflectivity and transmissivity for medium 1 into medium 2 (air to leaf) (similar arguments for $r_{21}$ and $t_{21}$. $\tau$ is the fraction of light transmitted through the medium. A similar expression can be derived from the leaf transmissivity, $T$ (not shown). $r_{12}$ is calculated from Fresnel's equations for an incidence angle $\theta$ and a refractive index $n$. We note that $r_{12} = 1- t_{12}$, $r_{21} = 1- t_{21}$ and that $t_{21}=t_{12}/n$. $\tau$ is related to the absorption coefficient of the plate, $k$, as is calculated through Beer's law. This means that in order to calculate the leaf transmittance and reflectance with this simple model, we only need the refraction index $n$ and the absorption coefficient $k$.
An improvement...¶
While the model presented above was successfully tested in a number of simple leaves and found to work reasonably well, it fails to account for more general leaves that exhibit a complex internal structure (such as senescent leaves, dicots...). To this end, the single plate model was extended to having a stack of $N$ plates and $N-1$ air layers between them, as shown below.
The solution to this new model was due to Stokes, and one finds that $T(N)$ and $R(N)$, respectively the transmittance and reflectance of a stack of $N$ layers is given by the following relationship
$$ \frac{R(N)}{b^{N}-b^{-N}}=\frac{T(N)}{a-a^{-1}}=\frac{1}{ab^{N}-a^{-1}b^{-N}}, $$where
$$ \begin{align} a &=\frac{1+R^{2}(1)-T^{2}(1) + \Delta}{2R(1)}\\ b &=\frac{1-R^{2}(1)+T^{2}(1) + \Delta}{2T(1)}\\ \Delta &= \sqrt { (T^{2}(1)-R^{2}(1)-1)-4R^{2}(1)}. \end{align} $$Here, $T(1)$ and $R(1)$ are the properties of a single layer, which can be calculated with a knowledge of the refraction index of the layer and the per layer absorption coefficient $k$. The latter is a summation of the weighted absorptions of the different leaf pigments.
The multilayer approach was further extended to consider the case where $N$ is a non-integer number, and resulted in the PROSPECT model. PROSPECT is able to reproduce measured leaf reflectance and transmittance spectra quite well, and has found a large user community. In order to parameterise the model, the initial version of PROSPECT used as input parameters the number of layers $N$ (remember this can be a non-integer number), the chlorophyll a+b concentration $[\mu gcm^{-2}]$ and the equivalent water thickness $[cm]$. In essence, one parameter controls the leaf internal structure, and the other two control the leaf optical properties in the visible (chlorophyll) and NIR and SWIR (water) regions.
Given that some compounds exhibit bonds that have absorption bands in the SWIR region (e.g. cellulose and lignin, proteins, ...), there was a concerted effort to add these contributions to the model, so in inversion studies, something about the carbon or nitrogen content of leaves could be retrieved. However, the lack of a very speicific signature for these components, and the fact that water absorption is an important part of the signal in the same region, results in very limited capabilities to infer these concentrations from either reflectance or transmittance. In practice, to account for these effects, a specific absorption spectrum of dry matter was developed, adding a new input to the prospect model, dry matter content $[gcm^{-2}]$. Finally, the chlorophyll contribution was split into chlorophyll and carotenoids, and an additional "brown pigment" was added (unitless). These constituents account for most of the variation in leaf optical properties, and thus make PROSPECT a good parsimonious model to model leaf transmittance and reflectance.
Some experiments with PROSPECT¶
In the next example, we'll try to see what the sensitivity of leaf reflectance and transmittance is to different PROSPECT parameters. A reasonably easy and intuitive way to achieve this is to do a "one at at time" study, where one of the inputs is sweeped and the output spectra are plotted. You can do this easily with a for loop. The PROSPECT model (version 5B) is part of the prosail package, which you should have installed. You can access the model with from prosail import prospect_5b. You can then check the (sparse!) documentation of the prospect_5b function by typing help( prospect_5b ). The parameter boundaries from Table 1 in Feret et al (2008) paper will also be useful:
| Table 1. | ||||
|---|---|---|---|---|
| Main characteristics of the datasets | ||||
| LOPEX | CALMIT | ANGERS | HAWAII | |
| Date | 1993 | 1996 | 2003 | 2007 |
| Number of samples | 64/580 | 106 | 276 | 41 |
| Number of species | 50 | 2 | 49 | 41 |
| Spectrophotometer/spectroradiometer | Perkin Elmer Lambda 19 | Shimadzu UV 2101 PC/LICOR LI-1800 | ASD FieldSpec | ASD FieldSpec |
| Spectral range | 400–2500 nm | 400–800 nm | 400–2450 nm | 400–2500 nm |
| Spectral sampling | 1–2 nm (400–1000 nm) | 1 nm | 1.4 nm (350–1050 nm) | 1.4 nm (350–1050 nm) |
| 4–5 nm (1000–2500 nm) | 2 nm (1000–2500 nm) | 2 nm (1000–2500 nm) | ||
| Solvent | Acetone 100% | Acetone 100% | Ethanol 95% | Acetone 100% |
| Method for pigments | Lichtenthaler (1987) | Lichtenthaler (1987) | Lichtenthaler (1987) | Lichtenthaler and Buschmann (2001) |
| Chlorophyll a (µg/cm2) | ||||
| Min | 0.3 | 1.0 | 0.4 | 13.3 |
| Max | 48.0 | 58.8 | 76.8 | 56.3 |
| Mean | 15.0 | 24.2 | 25.4 | 37.0 |
| Chlorophyll b (µg/cm2) | ||||
| Min | 0.2 | 0.2 | 0.3 | 5.4 |
| Max | 14.6 | 21.1 | 29.9 | 21.3 |
| Mean | 5 | 9 | 8 | 13 |
| Carotenoids (µg/cm2) | ||||
| Min | 0.6 | 2.0 | 0 | 5.1 |
| Max | 15.8 | 18.7 | 25.3 | 17.2 |
| Mean | 4.4 | 8.3 | 8.7 | 11.8 |
| Water (cm) | ||||
| Min | 0.0043 | n.a. | 0.0044 | 0.0127 |
| Max | 0.0439 | n.a. | 0.0340 | 0.0713 |
| Mean | 0.0113 | n.a. | 0.0116 | 0.0275 |
| Dry matter (g/cm2) | ||||
| Min | 0.0017 | n.a. | 0.0017 | 0.0064 |
| Max | 0.0152 | n.a. | 0.0331 | 0.0229 |
| Mean | 0.0053 | n.a. | 0.0052 | 0.0125 |
Questions¶
Use the prospect_sensivity_ssa function below to draw plots of the sensitivity of the different leaf constituents as a function of wavelength. The sensitivity, $S$, is calculated for each wavelength as
- Make sure the value of $\mathbf{x}_{0}$ you're using is sensible!
- Note where (spectrally) each constituent has a major impact.
- Are there any areas where a single constituent is the active?
- What are the consequences of areas where there is an overlap of effects from the different constituents?
- Make sure you try different locations in parameter space (e.g. $\mathbf{x}_0$). Are there any situations in which the plots change substantially from the original set?
from prosail_functions import *
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plot_config()
wv, sensitivity = prospect_sensitivity_ssa?
wv, sensitivity = prospect_sensitivity_ssa(x0=np.array([ 1.5 , 84, 5. , 0. , 0.011, 0.005]))
Exploring the parameter $N$¶
The number of leaf layers $N$ is an important parameter that stands in for the structural complexity of the internal leaf structure. The aim of the next experiment is to see what effect $N$ has in the ratio of reflecance to transmittance of the leaves.
- Comment on the form of the relationships found.
- Does the random variation of the other parameters seem to have much impact?
- What is the form of the relationship you see? Does the form vary with $N$? If so, how? What is the general trend with varying $N$?
wv = np.arange ( 400, 2501 )
xleafn, refl, trans = prospect_sensitivity_n( n_samples = 50)
fig, axs = plt.subplots ( nrows=5, ncols=2, figsize=(12,12), sharex=True )
axs=axs.flatten()
regr = []
for i,band in enumerate ( [ 560, 665, 705, 740, 783, 842, 865, 945, 1610, 2090 ] ):
axs[i].plot ( xleafn, np.mean(refl[:, wv==band]/trans[:, wv==band], axis=1),
'o',markerfacecolor="none")
y = np.clip ( np.mean(refl[:, wv==band]/trans[:, wv==band], axis=1), 0, 4 )
p = np.polyfit ( xleafn, y, 1 )
axs[i].plot ( xleafn, np.polyval (p, xleafn ), '-')
axs[i].set_title ("Band %d" % band )
axs[i].set_ylim (0, 4)
pretty_axes(axs[i])
Vegetation indices¶
Simple ratios of reflectance at different wavlenghts have been widely used to infer the properties of leaves and canopies. The goal of these indices is to maximise the sensitivity to a particular parameter of interest, while minimising the contribution of the other parameters. A typical form for these indices is a normalised vegetation index, e.g.
$$ VI = \frac{b_2 - b_1}{b_2 + b_1 }. $$The do_index function allows you to calculate an index like the one above over a sweep of a parameter of interest (e.g. $C_{ab}, C_{w}, \dots$), while randomly choosing the other parameters. The script allows you to select the centre wavelength for both bands 1 and 2, as well as the bandwidth (a band extends from $\lambda_c \pm \lambda_b$). It also allows you to calculate the index using either reflectance, transmittance or SSA
- Given your analysis of the sensitivity of PROSPECT, design VIs that maximise the impact of
- Chlorophyll
- Leaf water
- Dry matter
- Comment on the robustness of the relationship
- Consider reflectance, transmittance, single scattering albedo...
- What's the influence of spectral resolution? E.g. try to work out these indices for the typical spectral resolutions you find in e.g. Sentinel 2 or MODIS.
s,vi = do_index( sweep_param="cw", mode="ssa", band1=1610, band2=865, bwidth1=1, bwidth2=1, spaces=25, n_samples=50 )
Red edge model exploration¶
The leaf optical properties show a fast transition in the area known as the red edge, spectrally around 670-780 $nm$. In this transition region, chlorophyll goes from being strongly absorbing to the leaf internal structure being very reflective due to its internal structure. For high amounts of chlorophyll, the absorption feature in clorophyll around 680nm brings reflectance down, shifting the transition point between the chlorophyll absorption and the near infrared reflection to the right as a function of e.g. larger leaf chlorophyll concentrations. The shift of this transition zone has thus been related to "stress" on the leaves, but where is this "turning point"? One way to compress the data over this region is to fit an e.g. cubic model, and then finding where the turning point is. Remember that this turning point is given by the second derivative, $\partial f^{2}(\lambda)/\partial^2{\lambda}$. So if our model of reflectance is
$$ \rho = a + b\lambda + c\lambda^2 + d\lambda^3, $$then the second derivative is given by
$$ \frac{\partial^{2}\rho}{\partial^2\lambda}= 2c + 6d\lambda, $$and equating this to 0 yields $\lambda_{\textrm{red edge}} =-c/3d$. We have a function that fits this model to the data, and allows you to see how the location of $\lambda_{\textrm{red edge}}$ varies.
The red_edge function allows you to explore this concept. The code implements the rationale introduced above. You will see how the spectra, and the location of the inflexion point in the plots. You will also see the relationship of $\lambda_{\textrm{red edge}}$ with chlorophyll concentration.
- Based on the sensitivity analysis, try to explain what's happening in the red edge results that you are seeing.
- Is this method robust when we start considering slightly different spectral regions to fit the model to?
red_edge ( )
plt.hold(True)
red_edge( band1=620, band2=740)
pretty_axes()
Fitting some leaf transmittance and reflectance spectral using PROSPECT¶
While you might be familiar with statistical inference studies that aim to map magnitudes of interest from data, having a physical model that describes the different processes that give rise to your measurements is in many cases a far more robust way of infering the properties of (in this case) leaves than a purely statistical approach. We will consider a simple cost function minimisation for the time being, and use this explore and gain some intuition in the area of inverse problems.
The problem at hand is that we want to estimate, from a set of leaf reflectance and transmittance measurements over the solar reflective domain, the PROSPECT input parameters that result in the measurements. Since our model is able to predict the observations, the obvious thing to do is to minimise the difference between predictions and observations. So $M_{\rho}(\mathbf{x}, \lambda)$ and $M_{\tau}(\mathbf{x}, \lambda)$ represent the predicted reflectance and transmittance (respectively) calculated by PROSPECT with an input vector $\mathbf{x}$ (where $\mathbf{x}=\left[ N, Cab, Car, Csen, Cw,Cm\right]$) at a wavelength $\lambda$. The measurements will be denoted by $\rho(\lambda)$ and $\tau(\lambda)$, for refletance and transmittance (respectively). The least squares problem is then given by the minimisation of $J(\mathbf{x})$:
$$ J(\mathbf{x}) = \sum_{\lambda=400\,nm}^{2500\,nm} \left[ M_{rho}(\mathbf{x}, \lambda) - \rho(\lambda) \right]^{2} + \sum_{\lambda=400\,nm}^{2500\,nm} \left[ M_{tau}(\mathbf{x}, \lambda) - \tau(\lambda) \right]^{2}. $$We will be using some leaf spectra either from the field, or from the LOPEX'93 database. Below are two Python functions
read_lopex_samplereads a particular sample from the LOPEX database. Each sample has 5 replicates.optimise_random_startsminimises the misfit (e.g. $J(\mathbf{x})$ above) using a set of reflectance and transmittance measurements
You will need to provide your own implementation of the cost function (cost_function input parameters to optimise_random_starts).
refl, trans = read_lopex_sample ()
fwd_refl, fwd_trans, x, fun = optimise_random_starts (
the_cost_function, refl[2,:], trans[2,:] )
Some exercises...¶
- Try the above inversions for a couple of different samples, and try several of the replicates for each sample. Compare the results, what do you see?
- The five replicates are typically fairly noisy. We can take their difference as a statement of the measurement uncertainty. Can you suggest a way how you could modify the cost function in order to account for this?
- Uncertaint inputs certainly suggest uncertain outputs! Can you think of a way in which we could estimate the uncertainty in the outputs as a function of input uncertainty?
for samp in xrange(5):
fwd_refl, fwd_trans, x, fun = optimise_random_starts ( the_cost_function,
refl[samp,:], trans[samp,:], verbose=False, do_plots=False )
print samp+1, x, fun
from IPython.core.display import HTML
def css_styling():
styles = open("custom.css", "r").read()
return HTML(styles)
css_styling()