From 81fabc059f285c2abddeb4f2a2536110e5522901 Mon Sep 17 00:00:00 2001 From: David Orme Date: Mon, 14 Oct 2024 14:24:36 +0100 Subject: [PATCH] Revising calculations and docs for canopy --- docs/source/users/demography/canopy.md | 121 +++++++++++-------------- docs/source/users/demography/flora.md | 22 +++-- pyrealm/demography/canopy.py | 102 +++++++++++---------- 3 files changed, 125 insertions(+), 120 deletions(-) diff --git a/docs/source/users/demography/canopy.md b/docs/source/users/demography/canopy.md index 6894f11c..f875cfa7 100644 --- a/docs/source/users/demography/canopy.md +++ b/docs/source/users/demography/canopy.md @@ -179,11 +179,11 @@ ax1.set_xlabel("Profile radius (m)") ax1.set_ylabel("Vertical height (m)") # Plot the leaf area between heights for stems -ax2.plot(simple_canopy.layer_stem_leaf_area, hghts, color="red") +ax2.plot(simple_canopy.stem_leaf_area, hghts, color="red") ax2.set_xlabel("Leaf area (m2)") # Plot the fraction of light absorbed at different heights -ax3.plot(simple_canopy.layer_f_abs, hghts, color="red") +ax3.plot(simple_canopy.f_abs, hghts, color="red") ax3.set_xlabel("Light absorption fraction (-)") # Plot the light extinction profile through the canopy. @@ -254,8 +254,8 @@ community = Community( flora=flora, cell_area=32, cell_id=1, - cohort_dbh_values=np.array([0.05, 0.20, 0.5]), - cohort_n_individuals=np.array([7, 3, 1]), + cohort_dbh_values=np.array([0.1, 0.20, 0.5]), + cohort_n_individuals=np.array([7, 3, 2]), cohort_pft_names=np.array(["short", "short", "tall"]), ) @@ -284,6 +284,8 @@ print(canopy.extinction_profile[-1]) ``` ```{code-cell} ipython3 +:tags: [hide-input] + cols = ["r", "b", "g"] mpl.rcParams["axes.prop_cycle"] = mpl.cycler(color=cols) @@ -315,12 +317,12 @@ ax1.set_xlabel("Profile radius (m)") ax1.set_ylabel("Vertical height (m)") # Plot the leaf area between heights for stems -ax2.plot(canopy.layer_stem_leaf_area, hghts) +ax2.plot(canopy.stem_leaf_area, hghts) ax2.set_xlabel("Leaf area per stem (m2)") # Plot the fraction of light absorbed at different heights -ax3.plot(canopy.layer_f_abs, hghts, color="grey") -ax3.plot(1 - canopy.layer_cohort_f_trans, hghts) +ax3.plot(canopy.f_abs, hghts, color="grey") +ax3.plot(1 - canopy.cohort_f_trans, hghts) ax3.set_xlabel("Light absorption fraction (-)") # Plot the light extinction profile through the canopy. @@ -385,11 +387,11 @@ $$ canopy_ppa = Canopy(community=community, canopy_gap_fraction=0, fit_ppa=True) ``` -The `canopy_ppa.layer_heights` attribute now contains the heights at which the PPA +The `canopy_ppa.heights` attribute now contains the heights at which the PPA layers close: ```{code-cell} ipython3 -canopy_ppa.layer_heights +canopy_ppa.heights ``` And the final value in the canopy extinction profile still matches the expectation from @@ -426,27 +428,27 @@ to confirm that the calculated values coincide with the profile. Note here that total area at each closed layer height is omitting the community gap fraction. ```{code-cell} ipython3 +:tags: [hide-input] + fig, (ax1, ax2) = plt.subplots(ncols=2, sharey=True, figsize=(12, 6)) # Calculate the crown area at which each canopy layer closes. -closure_areas = ( - np.arange(1, len(canopy_ppa.layer_heights) + 1) * canopy.filled_community_area -) +closure_areas = np.arange(1, len(canopy_ppa.heights) + 1) * canopy.filled_community_area # LH plot - projected leaf area with height. # Add lines showing the canopy closure heights and closure areas. -for val in canopy_ppa.layer_heights: +for val in canopy_ppa.heights: ax1.axhline(val, color="red", linewidth=0.5, zorder=0) for val in closure_areas: ax1.axvline(val, color="red", linewidth=0.5, zorder=0) # Show the community projected crown area profile -ax1.plot(community_crown_area, canopy.layer_heights, zorder=1, label="Crown area") +ax1.plot(community_crown_area, canopy.heights, zorder=1, label="Crown area") ax1.plot( community_leaf_area, - canopy.layer_heights, + canopy.heights, zorder=1, linestyle="--", color="black", @@ -466,7 +468,7 @@ ax1.set_xlabel("Community-wide projected area (m2)") ax1.legend(frameon=False) # RH plot - light extinction -for val in canopy_ppa.layer_heights: +for val in canopy_ppa.heights: ax2.axhline(val, color="red", linewidth=0.5, zorder=0) for val in canopy_ppa.extinction_profile: @@ -489,111 +491,96 @@ ax2.set_xlabel("Light extinction (-)") ax2_rhs = ax2.twinx() ax2_rhs.set_ylim(ax2.get_ylim()) z_star_labels = [ - f"$z^*_{l + 1} = {val:.2f}$" - for l, val in enumerate(np.nditer(canopy_ppa.layer_heights)) + f"$z^*_{l + 1} = {val:.2f}$" for l, val in enumerate(np.nditer(canopy_ppa.heights)) ] -ax2_rhs.set_yticks(canopy_ppa.layer_heights.flatten()) +ax2_rhs.set_yticks(canopy_ppa.heights.flatten()) _ = ax2_rhs.set_yticklabels(z_star_labels) ``` ## Light allocation -In order to use the light extinction with the P Model, we need to calculate the -photosynthetic photon flux density (PPFD) captured within each layer and each cohort. + + +In order to use the light extinction with the P Model, we need to calculate the fraction +of absorbed photosynthetically active radiation $f_{APAR}$ within each layer for each +cohort. These values can be multiplied by the canopy-top photosynthetic photon flux +density (PPFD) to give the actual light absorbed for photosynthesis. + The steps below show this partitioning process for the PPA layers calculated above. 1. Calculate the fraction of light transmitted $f_{tr}$ through each layer for each cohort. The two arrays below show the extinction coefficients for the PFT of each cohort and then the cohort LAI ($L_H$, columns) components within each layer (rows). - The transmission through each component is then $f_{tr}=e^{-kL_H}$. + The transmission through each component is then $f_{tr}=e^{-kL_H}$ and + $f_{abs} = 1 - f_{tr}$ . ```{code-cell} ipython3 print("k = \n", community.stem_traits.par_ext, "\n") -print("L_H = \n", canopy_ppa.layer_cohort_lai) +print("L_H = \n", canopy_ppa.cohort_lai) ``` ```{code-cell} ipython3 -layer_cohort_f_tr = np.exp(-community.stem_traits.par_ext * canopy_ppa.layer_cohort_lai) +layer_cohort_f_tr = np.exp(-community.stem_traits.par_ext * canopy_ppa.cohort_lai) print(layer_cohort_f_tr) ``` -The fraction absorbed $f_{abs} = 1 - f_{tr}$. - ```{code-cell} ipython3 layer_cohort_f_abs = 1 - layer_cohort_f_tr print(layer_cohort_f_abs) ``` -These matrices show that there is complete transmission ($f_{abs} = 0, f_{tr} = 1$) -where a given stem has no leaf area within the layer but otherwise the leaves of each -stem absorb some light. - -If we want to calculate the total transmission within each layer, we need to sum the -individual cohort contributions within the exponent: - -```{code-cell} ipython3 -np.exp(np.sum(-community.stem_traits.par_ext * canopy_ppa.layer_cohort_lai, axis=1)) -``` + These matrices show that there is complete transmission ($f_{abs} = 0, f_{tr} = 1$) + where a given stem has no leaf area within the layer but otherwise the leaves of each + stem absorb some light. -Or alternatively, take the product within layers of the cohort components: +2. Calculate the total transmission across cohorts within each layer, as the product of + the individual cohort transmission within the layers, and then the absorption within + each layer ```{code-cell} ipython3 layer_f_tr = np.prod(layer_cohort_f_tr, axis=1) print(layer_f_tr) ``` -From that we can calculate $f_{abs} for each layer. - ```{code-cell} ipython3 layer_f_abs = 1 - layer_f_tr print(layer_f_abs) ``` -We can also calculate the cumulative fraction of light transmitted through the layers: +3. Calculate the transmission and extinction profiles through the layers as the + cumulative product of light transmitted. ```{code-cell} ipython3 transmission_profile = np.cumprod(layer_f_tr) print(transmission_profile) ``` -And then the extinction profile: - ```{code-cell} ipython3 extinction_profile = 1 - transmission_profile print(extinction_profile) ``` -The last thing to do is then calculate how much light is absorbed within each cohort. -The light can be partitioned into the light absorbed within each layer and reaching the -ground as follows: +4. Calculate the fraction of light transmitted through each each layer: ```{code-cell} ipython3 -ppfd = 1000 -ppfd_layer = -np.diff(ppfd * transmission_profile, prepend=ppfd, append=0) -print(ppfd_layer) +layer_fapar = -np.diff(transmission_profile, prepend=1) +print(layer_fapar) ``` -The cumulative sum across those quantities accounts for all of the incoming light and -matches the scaling of the extinction profile: +5. Calculate the relative absorbance across cohort within each layer and then use this + to partition the light absorbed in that layer across the cohorts: ```{code-cell} ipython3 -print(np.cumsum(ppfd_layer)) +cohort_fapar = ( + layer_cohort_f_abs / layer_cohort_f_abs.sum(axis=1)[:, None] +) * layer_fapar[:, None] +print(cohort_fapar) ``` -:::{admonition} TODO - -Now need to work out what the correct partition is of this within cohorts. It might -simply be weighted by relative $f_{abs}$ by cohort (i.e. `layer_cohort_f_abs`) but I'm -not 100% convinced! - -::: +6. Last, divide the cohort $f_{APAR}$ through by the number of individuals in each + cohort to given the $f_{APAR}$ for each stem at each height. -## Things to worry about later - -Herbivory - leaf fall (transmission increases, truncate at 0, replace from NSCs) vs leaf -turnover (transmission constant, increased GPP penalty) - -Leaf area dynamics in PlantFATE - acclimation to herbivory and transitory decreases in -transimission, need non-structural carbohydrates to recover from total defoliation. - -Leaf economics. +```{code-cell} ipython3 +stem_fapar = cohort_fapar / community.cohort_data["n_individuals"] +print(stem_fapar) +``` diff --git a/docs/source/users/demography/flora.md b/docs/source/users/demography/flora.md index 541ad5e2..8db0615d 100644 --- a/docs/source/users/demography/flora.md +++ b/docs/source/users/demography/flora.md @@ -10,6 +10,16 @@ kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 +language_info: + codemirror_mode: + name: ipython + version: 3 + file_extension: .py + mimetype: text/x-python + name: python + nbconvert_exporter: python + pygments_lexer: ipython3 + version: 3.11.9 --- # Plant Functional Types and Traits @@ -24,7 +34,7 @@ notes and initial demonstration code. This page introduces the main components of the {mod}`~pyrealm.demography` module that describe plant functional types (PFTs) and their traits. -```{code-cell} +```{code-cell} ipython3 from matplotlib import pyplot as plt import numpy as np import pandas as pd @@ -101,7 +111,7 @@ their maximum height. Note that the `q_m` and `z_max_prop` traits are calculated from the `m` and `n` traits and cannot be set directly. -```{code-cell} +```{code-cell} ipython3 short_pft = PlantFunctionalType(name="short", h_max=10) medium_pft = PlantFunctionalType(name="medium", h_max=20) tall_pft = PlantFunctionalType(name="tall", h_max=30) @@ -109,7 +119,7 @@ tall_pft = PlantFunctionalType(name="tall", h_max=30) The traits values set for a PFT instance can then be shown: -```{code-cell} +```{code-cell} ipython3 short_pft ``` @@ -130,13 +140,13 @@ that will be used in a demographic simulation. It can be created directly by pro the list of {class}`~pyrealm.demography.flora.PlantFunctionalType` instances. The only requirement is that each PFT instance uses a different name. -```{code-cell} +```{code-cell} ipython3 flora = Flora([short_pft, medium_pft, tall_pft]) flora ``` -```{code-cell} +```{code-cell} ipython3 pd.DataFrame({k: getattr(flora, k) for k in flora.trait_attrs}) ``` @@ -154,7 +164,7 @@ within {class}`~pyrealm.demography.community.Community` objects. A `StemTraits` instance can be created directly by providing arrays for each trait, but is more easily created from a `Flora` object by providing a list of PFT names: -```{code-cell} +```{code-cell} ipython3 # Get stem traits for a range of stems stem_pfts = ["short", "short", "short", "medium", "medium", "tall"] stem_traits = flora.get_stem_traits(pft_names=stem_pfts) diff --git a/pyrealm/demography/canopy.py b/pyrealm/demography/canopy.py index dd38854c..49f47176 100644 --- a/pyrealm/demography/canopy.py +++ b/pyrealm/demography/canopy.py @@ -203,36 +203,42 @@ def __init__( """Numerical tolerance for fitting the PPA model of canopy layer closure.""" # Define class attributes - self.total_community_crown_area: float - """Total crown area across all individuals in the community (m2).""" self.max_stem_height: float """Maximum height of any individual in the community (m).""" self.n_layers: int """Total number of canopy layers.""" self.n_cohorts: int """Total number of cohorts in the canopy.""" - self.layer_heights: NDArray[np.float32] - """Column vector of the heights of canopy layers.""" + self.heights: NDArray[np.float32] + """The vertical heights at which the canopy structure is calculated.""" + self.crown_profile: CrownProfile - """The crown profiles of stems at layer heights.""" - self.layer_stem_leaf_area: NDArray[np.float32] - """The leaf area of an individual stem per cohorts by layer.""" - self.layer_cohort_leaf_area: NDArray[np.float32] - """The total leaf areas within cohorts by layer.""" - self.layer_cohort_lai: NDArray[np.float32] + """The crown profiles of the community stems at the provided layer heights.""" + self.stem_leaf_area: NDArray[np.float32] + """The leaf area of the crown model for each cohort by layer.""" + self.cohort_lai: NDArray[np.float32] """The leaf area index for each cohort by layer.""" - self.layer_cohort_f_abs: NDArray[np.float32] + self.cohort_f_trans: NDArray[np.float32] + """The fraction of light transmitted by each cohort by layer.""" + self.cohort_f_abs: NDArray[np.float32] """The fraction of light absorbed by each cohort by layer.""" - self.layer_canopy_f_abs: NDArray[np.float32] - """The canopy-wide fraction of light absorbed by layer.""" - self.canopy_extinction_profile: NDArray[np.float32] - """The canopy-wide light extinction profile by layer.""" - self.fapar_profile: NDArray[np.float32] - """The canopy-wide fAPAR profile by layer.""" + self.f_trans: NDArray[np.float32] + """The fraction of light transmitted by the whole community by layer.""" + self.f_abs: NDArray[np.float32] + """The fraction of light absorbed by the whole community by layer.""" + self.transmission_profile: NDArray[np.float32] + """The light transmission profile for the whole community by layer.""" + self.extinction_profile: NDArray[np.float32] + """The light extinction profile for the whole community by layer.""" + self.fapar: NDArray[np.float32] + """The fraction of absorbed radiation for the whole community by layer.""" + self.cohort_fapar: NDArray[np.float32] + """The fraction of absorbed radiation for each cohort by layer.""" self.stem_fapar: NDArray[np.float32] - """The fAPAR for individual stems by layer.""" + """The fraction of absorbed radiation for each stem by layer.""" self.filled_community_area: float - """The area filled by crown.""" + """The area filled by crown after accounting for the crown gap fraction.""" + # Check operating mode if fit_ppa ^ (layer_heights is None): raise ValueError("Either set fit_ppa=True or provide layer heights.") @@ -246,9 +252,9 @@ def __init__( # Populate layer heights if layer_heights is not None: - self.layer_heights = layer_heights + self.heights = layer_heights else: - self.layer_heights = fit_perfect_plasticity_approximation( + self.heights = fit_perfect_plasticity_approximation( community=community, canopy_gap_fraction=canopy_gap_fraction, max_stem_height=self.max_stem_height, @@ -273,46 +279,48 @@ def _calculate_canopy(self, community: Community) -> None: self.crown_profile = CrownProfile( stem_traits=community.stem_traits, stem_allometry=community.stem_allometry, - z=self.layer_heights, + z=self.heights, ) # Partition the projected leaf area into the leaf area in each layer for each # stem and then scale up to the cohort leaf area in each layer. - self.layer_stem_leaf_area = np.diff( + self.stem_leaf_area = np.diff( self.crown_profile.projected_leaf_area, axis=0, prepend=0 ) # Calculate the leaf area index per layer per stem, using the stem # specific leaf area index values. LAI is a value per m2, so scale back down by # the available community area. - self.layer_cohort_lai = ( - self.layer_stem_leaf_area + self.cohort_lai = ( + self.stem_leaf_area * community.cohort_data["n_individuals"] * community.stem_traits.lai ) / community.cell_area # self.filled_community_area - # Calculate the Beer-Lambert light transmission component per layer and cohort - self.layer_cohort_f_trans = np.exp( - -community.stem_traits.par_ext * self.layer_cohort_lai - ) + # Calculate the Beer-Lambert light transmission and absorption components per + # layer and cohort + self.cohort_f_trans = np.exp(-community.stem_traits.par_ext * self.cohort_lai) + self.cohort_f_abs = 1 - self.cohort_f_trans + # Aggregate across cohorts into a layer wide transimissivity - self.layer_f_trans = self.layer_cohort_f_trans.prod(axis=1) + self.f_trans = self.cohort_f_trans.prod(axis=1) # Calculate the canopy wide light extinction per layer - self.layer_f_abs = 1 - self.layer_f_trans - - # Calculate cumulative light extinction across the canopy - self.extinction_profile = 1 - np.cumprod(self.layer_f_trans) - - # Calculate the fraction of radiation absorbed by each layer - # # TODO - include append=0 here to include ground or just backcalculate - self.fapar_profile = -np.diff( - np.cumprod(1 - self.layer_cohort_f_trans), - prepend=1, # append=0 - ) - - # # Partition the light back among the individual stems: simply weighting by per - # # cohort contribution to f_abs and divide through by the number of individuals - # self.stem_fapar = ( - # self.layer_cohort_f_trans * self.fapar_profile[:, None] - # ) / self.layer_cohort_f_trans.sum(axis=1)[:, None] + self.f_abs = 1 - self.f_trans + + # Calculate cumulative light transmission and extinction profiles + self.transmission_profile = np.cumprod(self.f_trans) + self.extinction_profile = 1 - self.transmission_profile + + # Calculate the fapar profile across cohorts and layers + # * The first part of the equation is calculating the relative absorption of + # each cohort within each layer + # * Each layer is then multiplied by fraction of the total light absorbed in the + # layer + # * The resulting matrix can be multiplied by a canopy top PPFD to generate the + # flux absorbed within each layer for each cohort. + self.fapar = -np.diff(self.transmission_profile, prepend=1) + self.cohort_fapar = ( + self.cohort_f_abs / self.cohort_f_abs.sum(axis=1)[:, None] + ) * self.fapar[:, None] + self.stem_fapar = self.cohort_fapar / community.cohort_data["n_individuals"]