From e3e5769308333a0bb90448b92e1adc760f812c40 Mon Sep 17 00:00:00 2001 From: David Orme Date: Wed, 2 Oct 2024 14:21:11 +0100 Subject: [PATCH 01/24] Initial functionality from light extinction doc draft --- pyrealm/demography/canopy.py | 84 ++++++++++++++++++++---------------- 1 file changed, 47 insertions(+), 37 deletions(-) diff --git a/pyrealm/demography/canopy.py b/pyrealm/demography/canopy.py index 05c35631..d27dd824 100644 --- a/pyrealm/demography/canopy.py +++ b/pyrealm/demography/canopy.py @@ -6,9 +6,7 @@ from pyrealm.demography.community import Community from pyrealm.demography.crown import ( - calculate_relative_crown_radius_at_z, - calculate_stem_projected_crown_area_at_z, - calculate_stem_projected_leaf_area_at_z, + CrownProfile, solve_community_projected_canopy_area, ) @@ -57,12 +55,8 @@ def __init__( """Total number of cohorts in the canopy.""" self.layer_heights: NDArray[np.float32] """Column vector of the heights of canopy layers.""" - self.stem_relative_radius: NDArray[np.float32] - """Relative radius values of stems at canopy layer heights.""" - self.stem_crown_area: NDArray[np.float32] - """Stem projected crown area at canopy layer heights.""" - self.stem_leaf_area: NDArray[np.float32] - """Stem projected leaf area at canopy layer heights.""" + self.crown_profile: CrownProfile + """The crown profiles of stems at layer heights.""" self._calculate_canopy(community=community) @@ -130,37 +124,53 @@ def _calculate_canopy(self, community: Community) -> None: self.layer_heights[layer] = starting_guess = solution.root - # Find relative canopy radius at the layer heights - # NOTE - here and in the calls below, validate=False is enforced because the - # Community class structures and code should guarantee valid inputs and so - # turning off the validation internally should simply speed up the code. - self.stem_relative_radius = calculate_relative_crown_radius_at_z( + # Calculate the crown profile at the layer heights + # TODO - reimpose validation + self.crown_profile = CrownProfile( + stem_traits=community.stem_traits, + stem_allometry=community.stem_allometry, z=self.layer_heights, - stem_height=community.stem_allometry.stem_height, - m=community.stem_traits.m, - n=community.stem_traits.n, - validate=False, ) - # Calculate projected crown area of a cohort stem at canopy closure heights. - self.stem_crown_area = calculate_stem_projected_crown_area_at_z( - z=self.layer_heights, - q_z=self.stem_relative_radius, - crown_area=community.stem_allometry.crown_area, - stem_height=community.stem_allometry.stem_height, - q_m=community.stem_traits.q_m, - z_max=community.stem_allometry.crown_z_max, - validate=False, + # self.canopy_projected_crown_area + # self.canopy_projected_leaf_area + + # 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.crown_profile.projected_leaf_area, axis=0, prepend=0 + ) + self.layer_cohort_leaf_area = ( + self.layer_stem_leaf_area * community.cohort_data["n_individuals"] ) - # Find the projected leaf area of a cohort stem at canopy closure heights. - self.stem_leaf_area = calculate_stem_projected_leaf_area_at_z( - z=self.layer_heights, - q_z=self.stem_relative_radius, - crown_area=community.stem_allometry.crown_area, - stem_height=community.stem_allometry.stem_height, - f_g=community.stem_traits.f_g, - q_m=community.stem_traits.q_m, - z_max=community.stem_allometry.crown_z_max, - validate=False, + # Calculate the leaf area index per layer per cohort, using the stem + # specific leaf area index values. LAI is a value per m2, so scale back down by + # the community area. + self.layer_cohort_lai = ( + self.layer_cohort_leaf_area * community.stem_traits.lai + ) / community.cell_area + + # Calculate the Beer-Lambert light extinction per layer and cohort + self.layer_cohort_f_abs = 1 - np.exp( + -community.stem_traits.par_ext * self.layer_cohort_lai ) + + # Calculate the canopy wide light extinction per layer + self.layer_canopy_f_abs = self.layer_cohort_f_abs.sum(axis=1) + + # Calculate cumulative light extinction across the canopy + self.canopy_extinction_profile = np.cumprod(self.layer_canopy_f_abs) + + # 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_canopy_f_abs), + 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_abs * self.fapar_profile[:, None] + ) / self.layer_cohort_f_abs.sum(axis=1)[:, None] From d260ebec6e69ba170645fbb83231fbaf60d2d669 Mon Sep 17 00:00:00 2001 From: David Orme Date: Wed, 2 Oct 2024 14:45:08 +0100 Subject: [PATCH 02/24] Moving canopy solver into canopy from crown, adding new canopy attributes and docstrings to __init__ --- pyrealm/demography/canopy.py | 91 +++++++++++++++++++++++++++++++++++- pyrealm/demography/crown.py | 69 --------------------------- 2 files changed, 89 insertions(+), 71 deletions(-) diff --git a/pyrealm/demography/canopy.py b/pyrealm/demography/canopy.py index d27dd824..7f484f73 100644 --- a/pyrealm/demography/canopy.py +++ b/pyrealm/demography/canopy.py @@ -7,10 +7,81 @@ from pyrealm.demography.community import Community from pyrealm.demography.crown import ( CrownProfile, - solve_community_projected_canopy_area, + _validate_z_qz_args, + calculate_relative_crown_radius_at_z, + calculate_stem_projected_crown_area_at_z, ) +def solve_canopy_area_filling_height( + z: float, + stem_height: NDArray[np.float32], + crown_area: NDArray[np.float32], + m: NDArray[np.float32], + n: NDArray[np.float32], + q_m: NDArray[np.float32], + z_max: NDArray[np.float32], + n_individuals: NDArray[np.float32], + target_area: float = 0, + validate: bool = True, +) -> NDArray[np.float32]: + """Solver function for finding the height where a canopy occupies a given area. + + This function takes the number of individuals in each cohort along with the stem + height and crown area and a given vertical height (:math:`z`). It then uses the + crown shape parameters associated with each cohort to calculate the community wide + projected crown area above that height (:math:`A_p(z)`). This is simply the sum of + the products of the individual stem crown projected area at :math:`z` and the number + of individuals in each cohort. + + The return value is the difference between the calculated :math:`A_p(z)` and a + user-specified target area, This allows the function to be used with a root solver + to find :math:`z` values that result in a given :math:`A_p(z)`. The default target + area is zero, so the default return value will be the actual total :math:`A_p(z)` + for the community. + + A typical use case for the target area would be to specify the area at which a given + canopy layer closes under the perfect plasticity approximation in order to find the + closure height. + + Args: + z: Vertical height on the z axis. + n_individuals: Number of individuals in each cohort + crown_area: Crown area of each cohort + stem_height: Stem height of each cohort + m: Crown shape parameter ``m``` for each cohort + n: Crown shape parameter ``n``` for each cohort + q_m: Crown shape parameter ``q_m``` for each cohort + z_max: Crown shape parameter ``z_m``` for each cohort + target_area: A target projected crown area. + validate: Boolean flag to suppress argument validation. + """ + # Convert z to array for validation and typing + z_arr = np.array(z) + + if validate: + _validate_z_qz_args( + z=z_arr, + stem_properties=[n_individuals, crown_area, stem_height, m, n, q_m, z_max], + ) + + q_z = calculate_relative_crown_radius_at_z( + z=z_arr, stem_height=stem_height, m=m, n=n, validate=False + ) + # Calculate A(p) for the stems in each cohort + A_p = calculate_stem_projected_crown_area_at_z( + z=z_arr, + q_z=q_z, + stem_height=stem_height, + crown_area=crown_area, + q_m=q_m, + z_max=z_max, + validate=False, + ) + + return (A_p * n_individuals).sum() - target_area + + class Canopy: """Model of the canopy for a plant community. @@ -57,6 +128,22 @@ def __init__( """Column vector of the heights of canopy layers.""" 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 leaf area index for each cohort by layer.""" + self.layer_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.stem_fapar: NDArray[np.float32] + """The fAPAR for individual stems by layer.""" self._calculate_canopy(community=community) @@ -101,7 +188,7 @@ def _calculate_canopy(self, community: Community) -> None: # TODO - the solution here is typically closer to the upper bracket, might # be a better algorithm to find the root (#293). solution = root_scalar( - solve_community_projected_canopy_area, + solve_canopy_area_filling_height, args=( community.stem_allometry.stem_height, community.stem_allometry.crown_area, diff --git a/pyrealm/demography/crown.py b/pyrealm/demography/crown.py index ed93fa32..bb3a01d8 100644 --- a/pyrealm/demography/crown.py +++ b/pyrealm/demography/crown.py @@ -217,75 +217,6 @@ def calculate_stem_projected_crown_area_at_z( return A_p -def solve_community_projected_canopy_area( - z: float, - stem_height: NDArray[np.float32], - crown_area: NDArray[np.float32], - m: NDArray[np.float32], - n: NDArray[np.float32], - q_m: NDArray[np.float32], - z_max: NDArray[np.float32], - n_individuals: NDArray[np.float32], - target_area: float = 0, - validate: bool = True, -) -> NDArray[np.float32]: - """Solver function for community wide projected canopy area. - - This function takes the number of individuals in each cohort along with the stem - height and crown area and a given vertical height (:math:`z`). It then uses the - crown shape parameters associated with each cohort to calculate the community wide - projected crown area above that height (:math:`A_p(z)`). This is simply the sum of - the products of the individual stem crown projected area at :math:`z` and the number - of individuals in each cohort. - - The return value is the difference between the calculated :math:`A_p(z)` and a - user-specified target area, This allows the function to be used with a root solver - to find :math:`z` values that result in a given :math:`A_p(z)`. The default target - area is zero, so the default return value will be the actual total :math:`A_p(z)` - for the community. - - A typical use case for the target area would be to specify the area at which a given - canopy layer closes under the perfect plasticity approximation in order to find the - closure height. - - Args: - z: Vertical height on the z axis. - n_individuals: Number of individuals in each cohort - crown_area: Crown area of each cohort - stem_height: Stem height of each cohort - m: Crown shape parameter ``m``` for each cohort - n: Crown shape parameter ``n``` for each cohort - q_m: Crown shape parameter ``q_m``` for each cohort - z_max: Crown shape parameter ``z_m``` for each cohort - target_area: A target projected crown area. - validate: Boolean flag to suppress argument validation. - """ - # Convert z to array for validation and typing - z_arr = np.array(z) - - if validate: - _validate_z_qz_args( - z=z_arr, - stem_properties=[n_individuals, crown_area, stem_height, m, n, q_m, z_max], - ) - - q_z = calculate_relative_crown_radius_at_z( - z=z_arr, stem_height=stem_height, m=m, n=n, validate=False - ) - # Calculate A(p) for the stems in each cohort - A_p = calculate_stem_projected_crown_area_at_z( - z=z_arr, - q_z=q_z, - stem_height=stem_height, - crown_area=crown_area, - q_m=q_m, - z_max=z_max, - validate=False, - ) - - return (A_p * n_individuals).sum() - target_area - - def calculate_stem_projected_leaf_area_at_z( z: NDArray[np.float32], q_z: NDArray[np.float32], From d0f2d72c49d4c483ea2653d9a53d6be649fd7de4 Mon Sep 17 00:00:00 2001 From: David Orme Date: Wed, 2 Oct 2024 15:35:28 +0100 Subject: [PATCH 03/24] Isolated the PPA model in its own function, made the canopy solver also accept an array of heights --- pyrealm/demography/canopy.py | 162 +++++++++++++++++++++-------------- 1 file changed, 98 insertions(+), 64 deletions(-) diff --git a/pyrealm/demography/canopy.py b/pyrealm/demography/canopy.py index 7f484f73..c4b5518a 100644 --- a/pyrealm/demography/canopy.py +++ b/pyrealm/demography/canopy.py @@ -82,6 +82,75 @@ def solve_canopy_area_filling_height( return (A_p * n_individuals).sum() - target_area +def fit_perfect_plasticity_approximation( + community: Community, + canopy_gap_fraction: float, + max_stem_height: float, + solver_tolerance: float, +) -> NDArray[np.float32]: + """Find canopy layer heights under the PPA model. + + Finds the closure heights of the canopy layers under the perfect plasticity + approximation by solving Ac(z) - L_n = 0 across the community where L is the + total cumulative crown area in layer n and above, discounted by the canopy gap + fraction. + + Args: + community: A community instance providing plant cohort data + canopy_gap_fraction: The canopy gap fraction + max_stem_height: The maximum stem height in the canopy, used as an upper bound + on finding the closure height of the topmost layer. + solver_tolerance: The absolute tolerance used with the root solver to find the + layer heights. + """ + + # Calculate the number of layers to contain the total community crown area + total_community_crown_area = ( + community.stem_allometry.crown_area * community.cohort_data["n_individuals"] + ).sum() + crown_area_per_layer = community.cell_area * (1 - canopy_gap_fraction) + n_layers = int(np.ceil(total_community_crown_area / crown_area_per_layer)) + + # Initialise the layer heights array and then loop over the layers indices, + # except for the final layer, which will be the partial remaining vegetation below + # the last closed layer. + layer_heights = np.zeros(n_layers, dtype=np.float32) + upper_bound = max_stem_height + + for layer in np.arange(n_layers - 1): + # Set the target area for this layer + target_area = (layer + 1) * crown_area_per_layer + + # TODO - the solution is typically closer to the upper bound of the bracket, + # there might be a better algorithm to find the root (#293). + solution = root_scalar( + solve_canopy_area_filling_height, + args=( + community.stem_allometry.stem_height, + community.stem_allometry.crown_area, + community.stem_traits.m, + community.stem_traits.n, + community.stem_traits.q_m, + community.stem_allometry.crown_z_max, + community.cohort_data["n_individuals"], + target_area, + False, # validate + ), + bracket=(0, upper_bound), + xtol=solver_tolerance, + ) + + if not solution.converged: + raise RuntimeError( + "Estimation of canopy layer closure heights failed to converge." + ) + + # Store the solution and update the upper bound for the next layer down. + layer_heights[layer] = upper_bound = solution.root + + return layer_heights + + class Canopy: """Model of the canopy for a plant community. @@ -106,20 +175,22 @@ class Canopy: def __init__( self, community: Community, + fit_ppa: bool = True, + layer_heights: NDArray[np.float32] | None = None, canopy_gap_fraction: float = 0.05, - layer_tolerance: float = 0.001, + solver_tolerance: float = 0.001, ) -> None: + # Store required init vars self.canopy_gap_fraction: float = canopy_gap_fraction """Canopy gap fraction.""" - self.layer_tolerance: float = layer_tolerance - """Numerical tolerance for solving canopy layer closure.""" + self.solver_tolerance: float = solver_tolerance + """Numerical tolerance for fitting the PPA model of canopy layer closure.""" + + # Define class attributes self.total_community_crown_area: float """Total crown area across individuals in the community (metres 2).""" self.max_stem_height: float """Maximum height of any individual in the community (metres).""" - self.crown_area_per_layer: float - """Total crown area permitted in a single canopy layer, given the available - cell area of the community and its canopy gap fraction.""" self.n_layers: int """Total number of canopy layers.""" self.n_cohorts: int @@ -145,72 +216,38 @@ def __init__( self.stem_fapar: NDArray[np.float32] """The fAPAR for individual stems by layer.""" + # Check operating mode + if fit_ppa ^ (layer_heights is None): + raise ValueError("Either set fit_ppa=True or provide layer heights.") + + # Set simple attributes + self.max_stem_height = community.stem_allometry.stem_height.max() + self.n_cohorts = community.number_of_cohorts + + # Populate layer heights + if layer_heights is not None: + self.layer_heights = layer_heights + else: + self.layer_heights = fit_perfect_plasticity_approximation( + community=community, + canopy_gap_fraction=canopy_gap_fraction, + max_stem_height=self.max_stem_height, + solver_tolerance=solver_tolerance, + ) + self._calculate_canopy(community=community) def _calculate_canopy(self, community: Community) -> None: """Calculate the canopy structure. This private method runs the calculations needed to populate the instance - attributes. + attributes, given the layer heights provided by the user or calculated using the + PPA model. Args: community: The Community object passed to the instance. """ - # Calculate community wide properties: total crown area, maximum height, crown - # area required to fill a layer and total number of canopy layers - self.total_community_crown_area = ( - community.stem_allometry.crown_area * community.cohort_data["n_individuals"] - ).sum() - - self.max_stem_height = community.stem_allometry.stem_height.max() - - self.crown_area_per_layer = community.cell_area * (1 - self.canopy_gap_fraction) - - self.n_layers = int( - np.ceil(self.total_community_crown_area / self.crown_area_per_layer) - ) - self.n_cohorts = community.number_of_cohorts - - # Find the closure heights of the canopy layers under the perfect plasticity - # approximation by solving Ac(z) - L_n = 0 across the community where L is the - # total cumulative crown area in layer n and above, discounted by the canopy gap - # fraction. - - self.layer_heights = np.zeros((self.n_layers, 1), dtype=np.float32) - - # Loop over the layers except for the final layer, which will be the partial - # remaining vegetation below the last closed layer. - starting_guess = self.max_stem_height - for layer in np.arange(self.n_layers - 1): - target_area = (layer + 1) * self.crown_area_per_layer - - # TODO - the solution here is typically closer to the upper bracket, might - # be a better algorithm to find the root (#293). - solution = root_scalar( - solve_canopy_area_filling_height, - args=( - community.stem_allometry.stem_height, - community.stem_allometry.crown_area, - community.stem_traits.m, - community.stem_traits.n, - community.stem_traits.q_m, - community.stem_allometry.crown_z_max, - community.cohort_data["n_individuals"], - target_area, - False, # validate - ), - bracket=(0, starting_guess), - xtol=self.layer_tolerance, - ) - - if not solution.converged: - raise RuntimeError( - "Estimation of canopy layer closure heights failed to converge." - ) - - self.layer_heights[layer] = starting_guess = solution.root - # Calculate the crown profile at the layer heights # TODO - reimpose validation self.crown_profile = CrownProfile( @@ -219,9 +256,6 @@ def _calculate_canopy(self, community: Community) -> None: z=self.layer_heights, ) - # self.canopy_projected_crown_area - # self.canopy_projected_leaf_area - # 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( From 53a6c76be60d18fd84a446c0a50b24912bf857eb Mon Sep 17 00:00:00 2001 From: David Orme Date: Thu, 3 Oct 2024 11:04:59 +0100 Subject: [PATCH 04/24] docstring updates --- pyrealm/demography/canopy.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/pyrealm/demography/canopy.py b/pyrealm/demography/canopy.py index c4b5518a..6bc9fe1b 100644 --- a/pyrealm/demography/canopy.py +++ b/pyrealm/demography/canopy.py @@ -148,17 +148,17 @@ def fit_perfect_plasticity_approximation( # Store the solution and update the upper bound for the next layer down. layer_heights[layer] = upper_bound = solution.root - return layer_heights + return layer_heights[:, None] class Canopy: - """Model of the canopy for a plant community. + """Calculate canopy characteristics for a plant community. This class generates a canopy structure for a community of trees using the - perfect-plasticity approximation model :cite:`purves:2008a`. In this approach, each - individual is assumed to arrange its canopy crown area plastically to take up space - in canopy layers and that new layers form below the canopy top as the available - space is occupied. + perfect-plasticity approximation (PPA) model :cite:`purves:2008a`. In this approach, + each individual is assumed to arrange its canopy crown area plastically to take up + space in canopy layers and that new layers form below the canopy top as the + available space is occupied. Real canopies contain canopy gaps, through process such as crown shyness. This is included in the model through the canopy gap fraction, which sets the proportion @@ -166,6 +166,10 @@ class Canopy: Args: community: A Community object that will be used to generate the canopy model. + layer_heights: A column array of vertical heights at which to calculate canopy + variables. + fit_ppa: Calculate layer heights as the canopy layer closure heights under the + PPA model. canopy_gap_fraction: The proportion of the available space unfilled by canopy (default: 0.05). layer_tolerance: The minimum precision used by the solver to find canopy layer @@ -175,8 +179,8 @@ class Canopy: def __init__( self, community: Community, - fit_ppa: bool = True, layer_heights: NDArray[np.float32] | None = None, + fit_ppa: bool = False, canopy_gap_fraction: float = 0.05, solver_tolerance: float = 0.001, ) -> None: From 9c1e31e30c38c101dbb73298d5ce9e5a78f2b70a Mon Sep 17 00:00:00 2001 From: David Orme Date: Thu, 3 Oct 2024 11:09:29 +0100 Subject: [PATCH 05/24] Crown doc changes --- docs/source/users/demography/crown.md | 437 ++++++++++++++++++++++++++ 1 file changed, 437 insertions(+) create mode 100644 docs/source/users/demography/crown.md diff --git a/docs/source/users/demography/crown.md b/docs/source/users/demography/crown.md new file mode 100644 index 00000000..eaa76c4c --- /dev/null +++ b/docs/source/users/demography/crown.md @@ -0,0 +1,437 @@ +--- +jupytext: + formats: md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# The tree crown model + +:::{admonition} Warning + +This area of `pyrealm` is in active development and this notebook currently contains +notes and initial demonstration code. + +::: + +```{code-cell} +from matplotlib import pyplot as plt +import numpy as np +import pandas as pd + +from pyrealm.demography.flora import ( + calculate_crown_q_m, + calculate_crown_z_max_proportion, + PlantFunctionalType, + Flora, +) + +from pyrealm.demography.t_model_functions import ( + calculate_dbh_from_height, + StemAllometry, +) + +from pyrealm.demography.crown import ( + CrownProfile, +) +``` + +The {mod}`pyrealm.demography` module uses three-dimensional model of crown shape to +define the vertical distribution of leaf area, following the implementation of crown +shape in the Plant-FATE model {cite}`joshi:2022a`. + +## Crown traits + +The crown model for a plant functional type (PFT) is driven by four traits within +the {class}`~pyrealm.demography.flora.PlantFunctionalType` class: + +* The `m` and `n` ($m, n$) traits set the vertical shape of the crown profile. +* The `ca_ratio` trait sets the area of the crown ($A_c$) relative to the stem size. +* The `f_g` ($f_g$) trait sets the crown gap fraction and sets the vertical + distribution of leaves within the crown. + +### Canopy shape + +For a stem of height $H$, the $m$ and $n$ traits are used to calculate the *relative* +crown radius $q(z)$ at a height $z$ of as: + +$$ +q(z)= m n \left(\dfrac{z}{H}\right) ^ {n -1} + \left( 1 - \left(\dfrac{z}{H}\right) ^ n \right)^{m-1} +$$ + +In order to align the arbitrary relative radius values with the predictions of the +T Model for the stem, we then need to find the height at which the relative radius +is at its maximum and a scaling factor that will convert the relative area at this +height to the expected crown area under the T Model. These can be calculated using +the following two additional traits, which are invariant for a PFT. + +* `z_max_prop` ($p_{zm}$) sets the proportion of the height of the stem at which + the maximum relative crown radius is found. +* `q_m` ($q_m$) is used to scale the size of the crown radius at $z_{max}$ to match + the expected crown area. + +$$ +\begin{align} +p_{zm} &= \left(\dfrac{n-1}{m n -1}\right)^ {\tfrac{1}{n}}\\[8pt] +q_m &= m n \left(\dfrac{n-1}{m n -1}\right)^ {1 - \tfrac{1}{n}} + \left(\dfrac{\left(m-1\right) n}{m n -1}\right)^ {m-1} +\end{align} +$$ + +For individual stems, with expected height $H$ and crown area $A_c$, we can then +calculate: + +* the height $z_m$ at which the maximum crown radius $r_m$ is found, +* a height-specific scaling factor $r_0$ such that $\pi q(z_m)^2 = A_c$, +* the actual crown radius $r(z)$ at a given height $z$. + +$$ +\begin{align} +z_m &= H p_{zm}\\[8pt] +r_0 &= \frac{1}{q_m}\sqrt{\frac{A_c}{\pi}}\\[8pt] +r(z) &= r_0 \; q(z) +\end{align} +$$ + +### Projected crown and leaf area + +From the crown radius profile, the model can then be used to calculate how crown area +and leaf area accumulates from the top of the crown towards the ground, giving +functions given height $z$ for: + +* the projected crown area $A_p(z)$, and +* the projected leaf area $\tilde{A}_{cp}(z)$. + +The projected crown area is calculated for a stem of known height and crown area is: + +$$ +A_p(z)= +\begin{cases} +A_c, & z \le z_m \\ +A_c \left(\dfrac{q(z)}{q_m}\right)^2, & H > z > z_m \\ +0, & z > H +\end{cases} +$$ + +That is, the projected crown area is zero above the top of the stem, increases to the +expected crown area at $z_max$ and is then constant to ground level. + +The projected leaf area $\tilde{A}_{cp}(z)$ models how the vertical distribution of +leaf area within the crown is modified by the crown gap fraction $f_g$. This trait +models how leaf gaps higher in the crown are filled by leaf area at lower heights: +it does not change the profile of the crown but allows leaf area in the crown to be +displaced downwards. When $f_g = 0$, the projected crown and leaf areas are identical, +but as $f_g \to 1$ the projected leaf area is pushed downwards. +The calculation of $\tilde{A}_{cp}(z)$ is defined as: + +$$ +\tilde{A}_{cp}(z)= +\begin{cases} +0, & z \gt H \\ +A_c \left(\dfrac{q(z)}{q_m}\right)^2 \left(1 - f_g\right), & H \gt z \gt z_m \\ +Ac - A_c \left(\dfrac{q(z)}{q_m}\right)^2 f_g, & zm \gt z +\end{cases} +$$ + +## Calculating crown model traits in `pyrealm` + +The {class}`~pyrealm.demography.flora.PlantFunctionalType` class is typically +used to set specific PFTs, but the functions to calculate $q_m$ and $p_{zm}$ +are used directly below to provides a demonstration of the impacts of each trait. + +```{code-cell} +# Set a range of values for m and n traits +m = n = np.arange(1.0, 5, 0.1) + +# Calculate the invariant q_m and z_max_prop traits +q_m = calculate_crown_q_m(m=m, n=n[:, None]) +z_max_prop = calculate_crown_z_max_proportion(m=m, n=n[:, None]) +``` + +```{code-cell} +fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10.9, 4)) + +# Plot q_m as a function of m and n +cntr_set1 = ax1.contourf(m, n, q_m, levels=10) +fig.colorbar(cntr_set1, ax=ax1, label="q_m") +ax1.set_xlabel("m") +ax1.set_ylabel("n") +ax1.set_aspect("equal") + +# Plot z_max_prop as a function of m and n +cntr_set2 = ax2.contourf(m, n, z_max_prop, levels=10) +fig.colorbar(cntr_set2, ax=ax2, label="z_max_prop") +ax2.set_xlabel("m") +ax2.set_ylabel("n") +ax2.set_aspect("equal") +``` + +## Crown calculations in `pyrealm` + +The examples below show the calculation of crown variables in `pyrealm`. + +### Calculting crown profiles + +The {class}`~pyrealm.demography.crown.CrownProfile` class is used to calculate crown +profiles for PFTs. It requires: + +* a {class}`~pyrealm.demography.flora.Flora` instance providing a set of PFTs to be + compared, +* a {class}`~pyrealm.demography.t_model_functions.StemAllometry` instance setting the + specific stem sizes for the profile, and +* a set of heights at which to estimate the profile variables + +The code below creates a set of PFTS with differing crown trait values and then creates +a `Flora` object using the PFTs. + +```{code-cell} +# A PFT with a small crown area and equal m and n values +narrow_pft = PlantFunctionalType(name="narrow", h_max=20, m=1.5, n=1.5, ca_ratio=20) +# A PFT with an intermediate crown area and m < n +medium_pft = PlantFunctionalType(name="medium", h_max=20, m=1.5, n=4, ca_ratio=500) +# A PFT with a wide crown area and m > n +wide_pft = PlantFunctionalType(name="wide", h_max=20, m=4, n=1.5, ca_ratio=2000) + +# Generate a Flora instance using those PFTs +flora = Flora([narrow_pft, medium_pft, wide_pft]) +flora +``` + +The Flora object can also be used to show a table of canopy variables: + +```{code-cell} +# TODO - add a Flora.to_pandas() method +flora_data = pd.DataFrame({k: getattr(flora, k) for k in flora.trait_attrs}) +flora_data[["name", "ca_ratio", "m", "n", "f_g", "q_m", "z_max_prop"]] +``` + +The next section of code generates the `StemAllometry` to use for the profiles. +The T Model requires DBH to define stem size - here the +{meth}`~pyrealm.demography.t_model_functions.calculate_dbh_from_height` function +is used to back-calculate the required DBH values to give three stems with similar +heights that are near the maximum height for each PFT. + +```{code-cell} +# Generate the expected stem allometries at similar heights for each PFT +stem_height = np.array([19, 17, 15]) +stem_dbh = calculate_dbh_from_height( + a_hd=flora.a_hd, h_max=flora.h_max, stem_height=stem_height +) +stem_dbh +``` + +```{code-cell} +# Calculate the stem allometries +allometry = StemAllometry(stem_traits=flora, at_dbh=stem_dbh) +``` + +We can again use {mod}`pandas` to get a table of those allometric predictions: + +```{code-cell} +pd.DataFrame({k: getattr(allometry, k) for k in allometry.allometry_attrs}) +``` + +Finally, we can define a set of vertical heights. In order to calculate the +variables for each PFT at each height, this needs to be provided as a column array, +that is with a shape `(N, 1)`. + +We can then calculate the crown profiles. + +```{code-cell} +# Create a set of vertical heights as a column array. +z = np.linspace(-1, 20.0, num=211)[:, None] + +# Calculate the crown profile across those heights for each PFT +crown_profiles = CrownProfile(stem_traits=flora, stem_allometry=allometry, z=z) +``` + +The `crown_profiles` object then provides the four crown profile attributes describe +above calculated at each height $z$: + +* The relative crown radius $q(z)$ +* The crown radius $r(z)$ +* The projected crown area +* The projected leaf area + +```{code-cell} +crown_profiles +``` + +### Visualising crown profiles + +The code below generates a plot of the vertical shape profiles of the crowns for each +stem. For each stem: + +* the dashed line shows how the relative crown radius $q(z)$ varies with height $z$, +* the solid line shows the actual crown radius $r(z)$ varies with height, and +* the dotted horizontal line shows the height at which the maximum crown radius is + found ($z_{max}$). + +Note that the equation for the relative radius $q(z)$ that defines canopy shape can +make predictions with actual values outside of the range of the actual stem and not +only where $0 \leq z \leq H$. + +```{code-cell} +fig, ax = plt.subplots(ncols=1) + +# Find the maximum of the actual and relative maximum crown widths +max_relative_crown_radius = np.nanmax(crown_profiles.relative_crown_radius, axis=0) +max_crown_radius = np.nanmax(crown_profiles.crown_radius, axis=0) +stem_max_width = np.maximum(max_crown_radius, max_relative_crown_radius) + + +for pft_idx, offset, colour in zip((0, 1, 2), (0, 5, 12), ("r", "g", "b")): + + # Plot relative radius either side of offset + stem_qz = crown_profiles.relative_crown_radius[:, pft_idx] + ax.plot(stem_qz + offset, z, color=colour, linestyle="--", linewidth=1) + ax.plot(-stem_qz + offset, z, color=colour, linestyle="--", linewidth=1) + + # Plot actual crown radius either side of offset + stem_rz = crown_profiles.crown_radius[:, pft_idx] + ax.plot(stem_rz + offset, z, color=colour) + ax.plot(-stem_rz + offset, z, color=colour) + + # Add the height of maximum crown radius + stem_rz_max = stem_max_width[pft_idx] + + ax.plot( + [offset - stem_rz_max, offset + stem_rz_max], + [allometry.crown_z_max[pft_idx]] * 2, + color=colour, + linewidth=1, + linestyle=":", + ) + + ax.set_xlabel("Crown profile") + ax.set_ylabel("Height above ground (m)") + +ax.set_aspect(aspect=1) +``` + +We can also use the `CanopyProfile` class with a single row of heights to calculate +the crown profile at the expected $z_{max}$ and show that this matches the expected +crown area from the T Model allometry. + +```{code-cell} +# Calculate the crown profile across those heights for each PFT +z_max = flora.z_max_prop * stem_height +profile_at_zmax = CrownProfile(stem_traits=flora, stem_allometry=allometry, z=z_max) + +print(profile_at_zmax.crown_radius**2 * np.pi) +print(allometry.crown_area) +``` + +### Visualising crown and leaf projected areas + +We can also use the crown profile to generate projected area plots. This is hard to see +using the PFTs defined above because they have very different crown areas, so the code +below generates new profiles for a new set of PFTs that have similar crown area ratios +but different shapes and gap fractions. + +```{code-cell} +no_gaps_pft = PlantFunctionalType( + name="no_gaps", h_max=20, m=1.5, n=1.5, f_g=0, ca_ratio=380 +) +few_gaps_pft = PlantFunctionalType( + name="few_gaps", h_max=20, m=1.5, n=4, f_g=0.05, ca_ratio=400 +) +many_gaps_pft = PlantFunctionalType( + name="many_gaps", h_max=20, m=4, n=1.5, f_g=0.2, ca_ratio=420 +) + +# Calculate allometries for each PFT at the same stem DBH +area_stem_dbh = np.array([0.4, 0.4, 0.4]) +area_flora = Flora([no_gaps_pft, few_gaps_pft, many_gaps_pft]) +area_allometry = StemAllometry(stem_traits=area_flora, at_dbh=area_stem_dbh) + +# Calculate the crown profiles across those heights for each PFT +area_z = np.linspace(0, area_allometry.stem_height.max(), 201)[:, None] +area_crown_profiles = CrownProfile( + stem_traits=area_flora, stem_allometry=area_allometry, z=area_z +) +``` + +The plot below then shows how projected crown area (solid lines) and leaf area (dashed +lines) change with height along the stem. + +* The projected crown area increases from zero at the top of the crown until it reaches + the maximum crown radius, at which point it remains constant down to ground level. The + total crown areas differs for each stem because of the slightly different values used + for the crown area ratio. + +* The projected leaf area is very similar but leaf area is displaced towards the ground + because of the crown gap fraction (`f_g`). Where `f_g = 0` (the red line), the two + lines are identical, but as `f_g` increases, more of the leaf area is displaced down + within the crown. + +```{code-cell} +fig, ax = plt.subplots(ncols=1) + +for pft_idx, offset, colour in zip((0, 1, 2), (0, 5, 10), ("r", "g", "b")): + + ax.plot(area_crown_profiles.projected_crown_area[:, pft_idx], area_z, color=colour) + ax.plot( + area_crown_profiles.projected_leaf_area[:, pft_idx], + area_z, + color=colour, + linestyle="--", + ) + ax.set_xlabel("Projected area (m2)") + ax.set_ylabel("Height above ground (m)") +``` + +We can also generate predictions for a single PFT with varying crown gap fraction. In +the plot below, note that all leaf area is above $z_{max}$ when $f_g=1$ and all leaf +area is *below* + +```{code-cell} +fig, ax = plt.subplots(ncols=1) + +# Loop over f_g values +for f_g in np.linspace(0, 1, num=11): + + label = None + colour = "gray" + + if f_g == 0: + label = "$f_g=0$" + colour = "red" + elif f_g == 1: + label = "$f_g=1$" + colour = "blue" + + # Create a flora with a single PFT with current f_g and then generate a + # stem allometry and crown profile + flora_f_g = Flora( + [PlantFunctionalType(name="example", h_max=20, m=2, n=2, f_g=f_g)] + ) + allometry_f_g = StemAllometry(stem_traits=flora_f_g, at_dbh=np.array([0.4])) + profile = CrownProfile( + stem_traits=flora_f_g, stem_allometry=allometry_f_g, z=area_z + ) + + # Plot the projected leaf area with height + ax.plot(profile.projected_leaf_area, area_z, color=colour, label=label, linewidth=1) + +# Add a horizontal line for z_max +ax.plot( + [-1, allometry_f_g.crown_area[0] + 1], + [allometry_f_g.crown_z_max, allometry_f_g.crown_z_max], + linestyle="--", + color="black", + label="$z_{max}$", + linewidth=1, +) + +ax.set_ylabel(r"Vertical height ($z$, m)") +ax.set_xlabel(r"Projected leaf area ($\tilde{A}_{cp}(z)$, m2)") +ax.legend(frameon=False) +``` From a71b38479c27b59d2924b9a4382b12b8bae812ad Mon Sep 17 00:00:00 2001 From: David Orme Date: Thu, 3 Oct 2024 16:13:18 +0100 Subject: [PATCH 06/24] Updating plot helper function get_crown_xy --- pyrealm/demography/crown.py | 130 ++++++++++++++++++++++++++++++++++-- 1 file changed, 125 insertions(+), 5 deletions(-) diff --git a/pyrealm/demography/crown.py b/pyrealm/demography/crown.py index f8833f57..9da04f51 100644 --- a/pyrealm/demography/crown.py +++ b/pyrealm/demography/crown.py @@ -3,6 +3,7 @@ """ # noqa: D205 from dataclasses import InitVar, dataclass, field +from typing import ClassVar import numpy as np from numpy.typing import NDArray @@ -306,11 +307,20 @@ class CrownProfile: z_max: A row array providing expected z_max height for each PFT. """ + var_attr_names: ClassVar[tuple[str, ...]] = ( + "relative_crown_radius", + "crown_radius", + "projected_crown_area", + "projected_leaf_area", + "projected_crown_radius", + "projected_leaf_radius", + ) + stem_traits: InitVar[StemTraits | Flora] """A Flora or StemTraits instance providing plant functional trait data.""" stem_allometry: InitVar[StemAllometry] """A StemAllometry instance setting the stem allometries for the crown profile.""" - z: InitVar[NDArray[np.float32]] + z: NDArray[np.float32] """An array of vertical height values at which to calculate crown profiles.""" relative_crown_radius: NDArray[np.float32] = field(init=False) @@ -332,12 +342,11 @@ def __post_init__( self, stem_traits: StemTraits | Flora, stem_allometry: StemAllometry, - z: NDArray[np.float32], ) -> None: """Populate crown profile attributes from the traits, allometry and height.""" # Calculate relative crown radius self.relative_crown_radius = calculate_relative_crown_radius_at_z( - z=z, + z=self.z, m=stem_traits.m, n=stem_traits.n, stem_height=stem_allometry.stem_height, @@ -350,7 +359,7 @@ def __post_init__( # Calculate projected crown area self.projected_crown_area = calculate_stem_projected_crown_area_at_z( - z=z, + z=self.z, q_z=self.relative_crown_radius, crown_area=stem_allometry.crown_area, q_m=stem_traits.q_m, @@ -360,7 +369,7 @@ def __post_init__( # Calculate projected leaf area self.projected_leaf_area = calculate_stem_projected_leaf_area_at_z( - z=z, + z=self.z, q_z=self.relative_crown_radius, f_g=stem_traits.f_g, q_m=stem_traits.q_m, @@ -383,3 +392,114 @@ def __repr__(self) -> str: f"CrownProfile: Prediction for {self._n_stems} stems " f"at {self._n_pred} observations." ) + + @property + def projected_crown_radius(self) -> NDArray[np.float32]: + """An array of the projected crown radius of stems at z heights.""" + return np.sqrt(self.projected_crown_area / np.pi) + + @property + def projected_leaf_radius(self) -> NDArray[np.float32]: + """An array of the projected leaf radius of stems at z heights.""" + return np.sqrt(self.projected_leaf_area / np.pi) + + +def get_crown_xy( + crown_profile: CrownProfile, + stem_allometry: StemAllometry, + attr: str, + stem_offsets: NDArray[np.float32] | None, + two_sided: bool = True, + as_xy: bool = False, +) -> list[tuple[NDArray, NDArray]] | list[NDArray]: + """Extract plotting data from crown profiles. + + A CrownProfile instance contains crown radius and projected area data for a set of + stems at given heights, but can contain predictions of these attributes above the + actual heights of some or all of the stems or indeed below ground. + + This function extracts plotting data for a given attribute for each crown that + includes only the predictions within the height range of the actual stem. It can + also mirror the values around the vertical midline to provide a two sided canopy + shape. + + The data are returned as a list with one entry per stem. The default value for each + entry a tuple of two arrays (height, attribute values) but the `as_xy=True` option + will return an `(N, 2)` dimensioned XY array suitable for use with + {class}`~matplotlib.patches.Polygon`. + + Args: + crown_profile: A crown profile instance + stem_allometry: The stem allometry instance used to create the crown profile + attr: The crown profile attribute to plot + stem_offsets: An optional array of offsets to add to the midline of stems. + two_sided: Should the plotting data show a two sided canopy. + as_xy: Should the plotting data be returned as a single XY array. + + """ + + # Input validation + if attr not in crown_profile.var_attr_names: + raise ValueError(f"Unknown crown profile attribute: {attr}") + + # TODO + # - more validation once the dimensioning has been thought through #317 + # - we're expecting a one d allometry and a 2D profile with multiple heights. + + # Get the attribute and flatten the heights from a column array to one dimensional + attr_values = getattr(crown_profile, attr) + z = crown_profile.z.flatten() + + # Orient the data so that lower heights always come first + if z[0] < z[-1]: + z = np.flip(z) + attr_values = np.flip(attr_values, axis=0) + + # Collect the per stem data + crown_plotting_data: list[tuple[NDArray, NDArray]] | list[NDArray] = [] + + for stem_index in np.arange(attr_values.shape[1]): + # Find the heights and values that fall within the individual stem + height_is_valid = np.logical_and( + z <= stem_allometry.stem_height[stem_index], z >= 0 + ) + valid_attr_values: NDArray = attr_values[height_is_valid, stem_index] + valid_heights: NDArray = z[height_is_valid] + + if two_sided: + # The values are extended to include the reverse profile as well as the zero + # value at the stem height + valid_heights = np.concatenate( + [ + np.flip(valid_heights), + stem_allometry.stem_height[[stem_index]], + valid_heights, + ] + ) + valid_attr_values = np.concatenate( + [-np.flip(valid_attr_values), [0], valid_attr_values] + ) + else: + # Only the zero value is added + valid_heights = np.concatenate( + [ + stem_allometry.stem_height[[stem_index]], + valid_heights, + ] + ) + valid_attr_values = np.concatenate([[0], valid_attr_values]) + + # Add offsets if provided + if stem_offsets is not None: + valid_attr_values += stem_offsets[stem_index] + + if as_xy: + # Combine the values into an XY array + crown_plotting_data.append( + np.hstack([valid_attr_values[:, None], valid_heights[:, None]]) + ) + else: + # Return the individual 1D arrays + crown_plotting_data.append((valid_heights, valid_attr_values)) + + return crown_plotting_data From 18a58da959e72af81063c5737f427c90bbbb33a8 Mon Sep 17 00:00:00 2001 From: David Orme Date: Thu, 3 Oct 2024 16:14:03 +0100 Subject: [PATCH 07/24] Updating crown.md --- docs/source/users/demography/crown.md | 50 +++++++++++++++++++++++++-- 1 file changed, 48 insertions(+), 2 deletions(-) diff --git a/docs/source/users/demography/crown.md b/docs/source/users/demography/crown.md index eaa76c4c..808b0b72 100644 --- a/docs/source/users/demography/crown.md +++ b/docs/source/users/demography/crown.md @@ -22,6 +22,7 @@ notes and initial demonstration code. ```{code-cell} from matplotlib import pyplot as plt +from matplotlib.patches import Polygon import numpy as np import pandas as pd @@ -39,6 +40,7 @@ from pyrealm.demography.t_model_functions import ( from pyrealm.demography.crown import ( CrownProfile, + get_crown_xy, ) ``` @@ -341,10 +343,10 @@ no_gaps_pft = PlantFunctionalType( name="no_gaps", h_max=20, m=1.5, n=1.5, f_g=0, ca_ratio=380 ) few_gaps_pft = PlantFunctionalType( - name="few_gaps", h_max=20, m=1.5, n=4, f_g=0.05, ca_ratio=400 + name="few_gaps", h_max=20, m=1.5, n=4, f_g=0.1, ca_ratio=400 ) many_gaps_pft = PlantFunctionalType( - name="many_gaps", h_max=20, m=4, n=1.5, f_g=0.2, ca_ratio=420 + name="many_gaps", h_max=20, m=4, n=1.5, f_g=0.3, ca_ratio=420 ) # Calculate allometries for each PFT at the same stem DBH @@ -435,3 +437,47 @@ ax.set_ylabel(r"Vertical height ($z$, m)") ax.set_xlabel(r"Projected leaf area ($\tilde{A}_{cp}(z)$, m2)") ax.legend(frameon=False) ``` + +```{code-cell} +# Set stem offsets for plotting +stem_offsets = np.array([0, 6, 12]) + +crown_radius_as_xy = get_crown_xy( + crown_profile=area_crown_profiles, + stem_allometry=area_allometry, + attr="crown_radius", + stem_offsets=stem_offsets, + as_xy=True, +) + +projected_crown_radius_xy = get_crown_xy( + crown_profile=area_crown_profiles, + stem_allometry=area_allometry, + attr="projected_crown_radius", + stem_offsets=stem_offsets, +) + +projected_leaf_radius_xy = get_crown_xy( + crown_profile=area_crown_profiles, + stem_allometry=area_allometry, + attr="projected_leaf_radius", + stem_offsets=stem_offsets, +) +``` + +```{code-cell} +fig, ax = plt.subplots() + +for cr_xy, (ch, cpr), (lh, lpr) in zip( + crown_radius_as_xy, projected_crown_radius_xy, projected_leaf_radius_xy +): + ax.add_patch(Polygon(cr_xy, color="lightgrey")) + ax.plot(cpr, ch, color="black", linewidth=1) + ax.plot(lpr, lh, color="red", linewidth=1) + +ax.set_aspect(0.5) +``` + +```{code-cell} + +``` From 4683fbfbb9e5d838099f4bfaa686fdd7e71aa9c1 Mon Sep 17 00:00:00 2001 From: David Orme Date: Thu, 3 Oct 2024 18:50:11 +0100 Subject: [PATCH 08/24] Add get_crown_xy example to crown doc --- docs/source/users/demography/crown.md | 60 ++++++++++++++++++++++----- pyrealm/demography/crown.py | 6 +++ 2 files changed, 56 insertions(+), 10 deletions(-) diff --git a/docs/source/users/demography/crown.md b/docs/source/users/demography/crown.md index 808b0b72..2653a5c0 100644 --- a/docs/source/users/demography/crown.md +++ b/docs/source/users/demography/crown.md @@ -22,7 +22,8 @@ notes and initial demonstration code. ```{code-cell} from matplotlib import pyplot as plt -from matplotlib.patches import Polygon +from matplotlib.lines import Line2D +from matplotlib.patches import Polygon, Patch import numpy as np import pandas as pd @@ -276,9 +277,20 @@ stem. For each stem: * the dotted horizontal line shows the height at which the maximum crown radius is found ($z_{max}$). -Note that the equation for the relative radius $q(z)$ that defines canopy shape can -make predictions with actual values outside of the range of the actual stem and not -only where $0 \leq z \leq H$. +:::{admonition} Note + +The predictions of the equation for the relative radius $q(z)$ are not limited to +height values within the range of the actual height of a given stem +($0 \leq z \leq H$). This is critical for calculating behaviour with height across +multiple stems when calculating canopy profiles for a community. The plot below +includes predictions of $q(z)$ below ground level and above stem height. + +The {meth}`~pyrealm.demography.crown.get_crown_xy` helper function can be used to +extract plotting structures for each stem within a `CrownProfile` that *are* +restricted to actual valid heights for that stem and is demonstrated in the +[code below](#plotting-tools-for-crown-shapes). + +::: ```{code-cell} fig, ax = plt.subplots(ncols=1) @@ -438,10 +450,30 @@ ax.set_xlabel(r"Projected leaf area ($\tilde{A}_{cp}(z)$, m2)") ax.legend(frameon=False) ``` +## Plotting tools for crown shapes + +The {meth}`~pyrealm.demography.crown.get_crown_xy` function makes it easier to extract +neat crown profiles from `CrownProfile` objects, for use in plotting crown data. The +function takes a paired `CrownProfile` and `StemAllometry` and extracts a particular +crown profile variable, and removes predictions for each stem that are outside of +the stem range for that stem. It converts the data for each stem into coordinates that +will plot as a complete two-sided crown outline. The returned value is a list with an +entry for each stem in one of two formats. + +* A pair of coordinate arrays: height and variable value. +* An single XY array with height and variable values in the columns, as used for + example in `matplotlib` Patch objects. + +The code below uses this function to generate plotting data for the crown radius, +projected crown radius and projected leaf radius. These last two variables do not +have direct computational use - the cumulative projected area is what matters - but +allow the projected variables to be visualised at the same scale as the crown radius. + ```{code-cell} -# Set stem offsets for plotting +# Set stem offsets for separating stems along the x axis stem_offsets = np.array([0, 6, 12]) +# Get the crown radius in XY format to plot as a polygon crown_radius_as_xy = get_crown_xy( crown_profile=area_crown_profiles, stem_allometry=area_allometry, @@ -450,6 +482,7 @@ crown_radius_as_xy = get_crown_xy( as_xy=True, ) +# Get the projected crown and leaf radii to plot as lines projected_crown_radius_xy = get_crown_xy( crown_profile=area_crown_profiles, stem_allometry=area_allometry, @@ -468,16 +501,23 @@ projected_leaf_radius_xy = get_crown_xy( ```{code-cell} fig, ax = plt.subplots() +# Bundle the three plotting structures and loop over the three stems. for cr_xy, (ch, cpr), (lh, lpr) in zip( crown_radius_as_xy, projected_crown_radius_xy, projected_leaf_radius_xy ): ax.add_patch(Polygon(cr_xy, color="lightgrey")) - ax.plot(cpr, ch, color="black", linewidth=1) + ax.plot(cpr, ch, color="0.4", linewidth=2) ax.plot(lpr, lh, color="red", linewidth=1) ax.set_aspect(0.5) -``` - -```{code-cell} - +plt.legend( + handles=[ + Patch(color="lightgrey", label="Crown profile"), + Line2D([0], [0], label="Projected crown", color="0.4", linewidth=2), + Line2D([0], [0], label="Projected leaf", color="red", linewidth=1), + ], + ncols=3, + loc="upper center", + bbox_to_anchor=(0.5, 1.15), +) ``` diff --git a/pyrealm/demography/crown.py b/pyrealm/demography/crown.py index 9da04f51..39311141 100644 --- a/pyrealm/demography/crown.py +++ b/pyrealm/demography/crown.py @@ -295,6 +295,12 @@ class CrownProfile: allometric predictions of stem height, crown area and z_max for an actual stem of a given size for each PFT. + In addition to the variables above, the class can also has properties the calculate + the projected crown radius and projected leaf radius. These are simply the radii + that would result in the two projected areas: the values are not directly meaningful + for calculating canopy models, but can be useful for exploring the behavour of + projected area on the same linear scale as the crown radius. + Args: stem_traits: A Flora or StemTraits instance providing plant functional trait data. From ffef774a9a5b8cdecccb2547bf4ae5c117d5eceb Mon Sep 17 00:00:00 2001 From: David Orme Date: Sat, 5 Oct 2024 11:37:00 +0100 Subject: [PATCH 09/24] Docstrings and comments --- pyrealm/demography/canopy.py | 20 ++++++++++++++++---- pyrealm/demography/crown.py | 2 +- 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/pyrealm/demography/canopy.py b/pyrealm/demography/canopy.py index 6bc9fe1b..1d57119d 100644 --- a/pyrealm/demography/canopy.py +++ b/pyrealm/demography/canopy.py @@ -88,12 +88,24 @@ def fit_perfect_plasticity_approximation( max_stem_height: float, solver_tolerance: float, ) -> NDArray[np.float32]: - """Find canopy layer heights under the PPA model. + r"""Find canopy layer heights under the PPA model. Finds the closure heights of the canopy layers under the perfect plasticity - approximation by solving Ac(z) - L_n = 0 across the community where L is the - total cumulative crown area in layer n and above, discounted by the canopy gap - fraction. + approximation by fidnding the set of heights that lead to complete closure of canopy + layers through the canopy. The function solves the following equation for integers + :math:`l \in (1,2,..., m)`: + + .. math:: + + \sum_{s=1}^{N_s}{ A_p(z^*_l)} = l A(1 - f_G) + + The right hand side sets out the total area needed to close a given layer :math:`l` + and all layers above it: :math:`l` times the total community area :math:`A` less + any canopy gap fraction (:math:`f_G`). The left hand side then calculates the + projected crown area for each stem :math:`s` :math:`A_p(z^*_l)_{[s]}` and sums those + areas across all stems in the community :math:`N_s`. The specific height + :math:`z^*_l` is then the height at which the two terms are equal and hence solves + the equation for layer :math:`l`. Args: community: A community instance providing plant cohort data diff --git a/pyrealm/demography/crown.py b/pyrealm/demography/crown.py index 39311141..9aaaee3f 100644 --- a/pyrealm/demography/crown.py +++ b/pyrealm/demography/crown.py @@ -500,7 +500,7 @@ def get_crown_xy( valid_attr_values += stem_offsets[stem_index] if as_xy: - # Combine the values into an XY array + # Combine the values into an (N,2) XY array crown_plotting_data.append( np.hstack([valid_attr_values[:, None], valid_heights[:, None]]) ) From 7602f6a12e8a2b6f96d486592e2509b3f36bf706 Mon Sep 17 00:00:00 2001 From: David Orme Date: Sat, 5 Oct 2024 11:46:18 +0100 Subject: [PATCH 10/24] Doc updates --- docs/source/users/demography/canopy.md | 72 +++++++++++++++++++++----- pyproject.toml | 16 +++--- 2 files changed, 68 insertions(+), 20 deletions(-) diff --git a/docs/source/users/demography/canopy.md b/docs/source/users/demography/canopy.md index f9c0688c..5fa428b1 100644 --- a/docs/source/users/demography/canopy.md +++ b/docs/source/users/demography/canopy.md @@ -1,17 +1,14 @@ --- jupytext: - formats: md:myst text_representation: - extension: .md - format_name: myst - format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- -# Canopy model +# The canopy model :::{admonition} Warning @@ -20,16 +17,39 @@ notes and initial demonstration code. ::: +The canopy is the representation of the combined crowns across all individual stems +within each cohort in a [community](./community.md). The main use of the canopy is +to calculate the light environment within the canopy, which allows the the productivity +of stems to be estimated across the vertical dimension. + +The key variables in calculating the canopy model are the crown projected area $A_p(z)$ +and leaf projected area $\tilde{A}_{cp}(z)$, which are calculated for a stem +of a given size at a vertical height $z$ using the [crown model](./crown.md). For a +given stem, the difference in projected leaf area values between two heights provides +the actual leaf area between two heights. This then gives an expression for the leaf +area index (LAI, $L$) between two heights $z_i$ and $z_j$: + +$$ + L_{i,j} = \frac{L\left(\tilde{A}_{cp}(z_j) - \tilde{A}_{cp}(z_i)\right)}{A}, +$$ + +where $z_i > z_j$ and $A$ is the area occupied by the community. The value of $L_{i,j}$ +is then the leaf area index per square metre for the stem between the two heights. The +light extinction between $z_i$ and $z_j$ then follows the Beer Lambert law, +with the fraction of light absorbed $f_{abs}$ calculated as: + +$$ +f_{abs[i,j]} = 1 - e^{-kL_{i,j}} +$$ + +given the extinction coefficient $k$ from the plant functional type of the stem. + The canopy model uses the perfect plasticity approximation (PPA) {cite}`purves:2008a`, which assumes that plants are always able to plastically arrange their crown within the broader canopy of the community to maximise their crown size and fill the available space $A$. When the area $A$ is filled, a new lower canopy layer is formed until all of the individual crown area has been distributed across within the canopy. -The key variables in calculating the canopy model are the crown projected area $A_p$ -and leaf projected area $\tilde{A}_{cp}(z)$, which are calculated for a stem -of a given size using the [crown model](./crown.md). - ```{code-cell} from matplotlib import pyplot as plt import numpy as np @@ -93,7 +113,7 @@ The code below creates a simple community: short_pft = PlantFunctionalType( name="short", h_max=15, m=1.5, n=1.5, f_g=0, ca_ratio=380 ) -tall_pft = PlantFunctionalType(name="tall", h_max=30, m=1.5, n=2, f_g=0.2, ca_ratio=500) +tall_pft = PlantFunctionalType(name="tall", h_max=30, m=3, n=1.5, f_g=0.2, ca_ratio=500) # Create the flora flora = Flora([short_pft, tall_pft]) @@ -107,7 +127,7 @@ community = Community( flora=flora, cell_area=32, cell_id=1, - cohort_dbh_values=np.array([0.02, 0.20, 0.5]), + cohort_dbh_values=np.array([0.05, 0.20, 0.5]), cohort_n_individuals=np.array([15, 5, 2]), cohort_pft_names=np.array(["short", "short", "tall"]), ) @@ -123,7 +143,35 @@ print("Ac = ", community.stem_allometry.crown_area) We can now calculate the canopy model for the community: ```{code-cell} -canopy = Canopy(community=community, canopy_gap_fraction=2 / 32) +:lines_to_next_cell: 2 + +hghts = np.linspace(26, 0, num=251)[:, None] +canopy = Canopy( + community=community, canopy_gap_fraction=2 / 32, layer_heights=hghts, fit_ppa=False +) +``` + +```{code-cell} +fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(ncols=5, sharey=True, figsize=(12, 6)) + +ax1.plot(canopy.fapar_profile, hghts) + +ax2.plot(np.cumsum(canopy.fapar_profile), hghts) +ax3.plot(canopy.layer_stem_leaf_area, hghts) +ax4.plot(canopy.crown_profile.projected_leaf_radius, hghts) + +ax5.set_prop_cycle("color", ["r", "b", "g"]) +ax5.plot(canopy.crown_profile.crown_radius, hghts) +ax5.plot(canopy.crown_profile.projected_leaf_radius, hghts, linestyle="--") + +# ax1.plot(canopy.crown_profile.projected_leaf_radius, hghts) +# ax1.plot(canopy.crown_profile.projected_crown_radius, hghts) +# ax1.set_aspect(1) +# plt.xscale('log') +``` + +```{code-cell} +canopy.stem_fapar.sum() ``` We can then look at three key properties of the canopy model: the layer closure diff --git a/pyproject.toml b/pyproject.toml index 2823f7cc..45aac833 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,15 +31,15 @@ repository = "https://github.com/ImperialCollegeLondon/pyrealm" version = "1.0.0" [tool.poetry.dependencies] -dacite = "^1.6.0" -numpy = "^2.0.0" +dacite = "^1.6.0" +numpy = "^2.0.0" python = ">=3.10" -scipy = "^1.7.3" -tabulate = "^0.8.10" +scipy = "^1.7.3" +tabulate = "^0.8.10" marshmallow = "^3.22.0" -pandas = "^2.2.2" marshmallow-dataclass = "^8.7.0" +pandas = "^2.2.2" pandas-stubs = "^2.2.2.240909" [tool.poetry.group.types.dependencies] pandas-stubs = "^2.2.0.240218" @@ -132,7 +132,7 @@ select = [ "I", # isort "UP", # pyupgrade "RUF", # RUF specific checks - "NPY201" + "NPY201", ] # On top of the Google convention, disable: @@ -147,6 +147,6 @@ convention = "google" [tool.jupytext] # Stop jupytext from removing mystnb and other settings in MyST Notebook YAML headers -notebook_metadata_filter = "-jupytext.text_representation.jupytext_version,settings,mystnb" +notebook_metadata_filter = "jupytext.text_representation.jupytext_version,settings,mystnb,language_info" # Also stop it from stripping cell metadata. -cell_metadata_filter = "all" \ No newline at end of file +cell_metadata_filter = "all" From 7f18eb8f6c283ff6b891ce962f0b4be2955308bd Mon Sep 17 00:00:00 2001 From: David Orme Date: Mon, 7 Oct 2024 09:54:09 +0100 Subject: [PATCH 11/24] Updating jupytext settings --- docs/source/users/demography/canopy.md | 7 +++++-- pyproject.toml | 4 ++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/docs/source/users/demography/canopy.md b/docs/source/users/demography/canopy.md index 5fa428b1..2dbc140a 100644 --- a/docs/source/users/demography/canopy.md +++ b/docs/source/users/demography/canopy.md @@ -1,6 +1,9 @@ --- jupytext: text_representation: + extension: .md + format_name: myst + format_version: 0.13 jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) @@ -103,6 +106,8 @@ The {class}`~pyrealm.demography.canopy.Canopy` class automatically finds the can closure heights, given a {class}`~pyrealm.demography.community.Community` instance and the required canopy gap fraction. +THIS IS A TEST CHANGE TO CHECK JUPYTEXT BEHAVIOUR. + The code below creates a simple community: ```{code-cell} @@ -143,8 +148,6 @@ print("Ac = ", community.stem_allometry.crown_area) We can now calculate the canopy model for the community: ```{code-cell} -:lines_to_next_cell: 2 - hghts = np.linspace(26, 0, num=251)[:, None] canopy = Canopy( community=community, canopy_gap_fraction=2 / 32, layer_heights=hghts, fit_ppa=False diff --git a/pyproject.toml b/pyproject.toml index 45aac833..b656af26 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -147,6 +147,6 @@ convention = "google" [tool.jupytext] # Stop jupytext from removing mystnb and other settings in MyST Notebook YAML headers -notebook_metadata_filter = "jupytext.text_representation.jupytext_version,settings,mystnb,language_info" +notebook_metadata_filter = "settings,mystnb,language_info" # Also stop it from stripping cell metadata. -cell_metadata_filter = "all" +cell_metadata_filter = "all,-trusted" From a1ad5f0c008964cc58856f9436d1182dfe6743dd Mon Sep 17 00:00:00 2001 From: David Orme Date: Mon, 7 Oct 2024 09:54:45 +0100 Subject: [PATCH 12/24] Testing jupytext settings --- docs/source/users/demography/canopy.md | 44 +++++++++++++++----------- 1 file changed, 26 insertions(+), 18 deletions(-) diff --git a/docs/source/users/demography/canopy.md b/docs/source/users/demography/canopy.md index 2dbc140a..9a9b9901 100644 --- a/docs/source/users/demography/canopy.md +++ b/docs/source/users/demography/canopy.md @@ -9,6 +9,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 --- # The canopy model @@ -53,7 +63,7 @@ broader canopy of the community to maximise their crown size and fill the availa space $A$. When the area $A$ is filled, a new lower canopy layer is formed until all of the individual crown area has been distributed across within the canopy. -```{code-cell} +```{code-cell} ipython3 from matplotlib import pyplot as plt import numpy as np import pandas as pd @@ -106,11 +116,9 @@ The {class}`~pyrealm.demography.canopy.Canopy` class automatically finds the can closure heights, given a {class}`~pyrealm.demography.community.Community` instance and the required canopy gap fraction. -THIS IS A TEST CHANGE TO CHECK JUPYTEXT BEHAVIOUR. - The code below creates a simple community: -```{code-cell} +```{code-cell} ipython3 # Two PFTs # - a shorter understory tree with a columnar canopy and no crown gaps # - a taller canopy tree with a top heavy canopy and more crown gaps @@ -140,21 +148,21 @@ community = Community( We can then look at the expected allometries for the stems in each cohort: -```{code-cell} +```{code-cell} ipython3 print("H = ", community.stem_allometry.stem_height) print("Ac = ", community.stem_allometry.crown_area) ``` We can now calculate the canopy model for the community: -```{code-cell} +```{code-cell} ipython3 hghts = np.linspace(26, 0, num=251)[:, None] canopy = Canopy( community=community, canopy_gap_fraction=2 / 32, layer_heights=hghts, fit_ppa=False ) ``` -```{code-cell} +```{code-cell} ipython3 fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(ncols=5, sharey=True, figsize=(12, 6)) ax1.plot(canopy.fapar_profile, hghts) @@ -173,7 +181,7 @@ ax5.plot(canopy.crown_profile.projected_leaf_radius, hghts, linestyle="--") # plt.xscale('log') ``` -```{code-cell} +```{code-cell} ipython3 canopy.stem_fapar.sum() ``` @@ -184,14 +192,14 @@ heights for each stem in the three cohorts. There are four canopy layers, with the top two very close together because of the large crown area in the two stems in the cohort of `tall` trees. -```{code-cell} +```{code-cell} ipython3 canopy.layer_heights ``` The `stem_crown_area` attribute then provides the crown area of each stem found in each layer. -```{code-cell} +```{code-cell} ipython3 canopy.stem_crown_area ``` @@ -200,7 +208,7 @@ the first two layers are taken up entirely by the two stems in the cohort of lar trees. We can confirm that the calculation is correct by calculating the total crown area across the cohorts at each height: -```{code-cell} +```{code-cell} ipython3 np.sum(canopy.stem_crown_area * community.cohort_data["n_individuals"], axis=1) ``` @@ -214,7 +222,7 @@ identical to the projected crown area for the first two cohorts because the crow fraction $f_g$ is zero for this PFT. The projected leaf area is however displaced towards the ground in the last cohort, because the `tall` PFT has a large gap fraction. -```{code-cell} +```{code-cell} ipython3 canopy.stem_leaf_area ``` @@ -225,7 +233,7 @@ community crown and leaf area profile across a range of height values. For each we calculate the sum of the product of stem projected area and the number of individuals in each cohort. -```{code-cell} +```{code-cell} ipython3 # Set of vertical height to calculate crown profiles at_z = np.linspace(0, 26, num=261)[:, None] @@ -249,7 +257,7 @@ superimpose the calculated $z^*_l$ values and the cumulative canopy area for eac to confirm that the calculated values coincide with the profile. Note here that the total area at each closed layer height is omitting the community gap fraction. -```{code-cell} +```{code-cell} ipython3 fig, ax = plt.subplots(ncols=1) # Calculate the crown area at which each canopy layer closes. @@ -307,7 +315,7 @@ $f_{abs} = 1 - e ^ {-kL}$, where $k$ is the light extinction coefficient ($k$) and $L$ is the leaf area index (LAI). The LAI can be calculated for each stem and layer: -```{code-cell} +```{code-cell} ipython3 # LAI = Acp_within_layer / canopy_area # print(LAI) ``` @@ -315,7 +323,7 @@ where $k$ is the light extinction coefficient ($k$) and $L$ is the leaf area ind This can be used to calculate the LAI of individual stems but also the LAI of each layer in the canopy: -```{code-cell} +```{code-cell} ipython3 # LAI_stem = LAI.sum(axis=0) # LAI_layer = LAI.sum(axis=1) @@ -326,7 +334,7 @@ in the canopy: The layer LAI values can now be used to calculate the light transmission of each layer and hence the cumulative light extinction profile through the canopy. -```{code-cell} +```{code-cell} ipython3 # f_abs = 1 - np.exp(-pft.traits.par_ext * LAI_layer) # ext = np.cumprod(f_abs) @@ -338,7 +346,7 @@ One issue that needs to be resolved is that the T Model implementation in `pyrea follows the original implementation of the T Model in having LAI as a fixed trait of a given plant functional type, so is constant for all stems of that PFT. -```{code-cell} +```{code-cell} ipython3 # print("f_abs = ", (1 - np.exp(-pft.traits.par_ext * pft.traits.lai))) ``` From bb8d3c9ca920f6503d32dfb3a772032ced8aeae0 Mon Sep 17 00:00:00 2001 From: David Orme Date: Wed, 9 Oct 2024 18:23:45 +0100 Subject: [PATCH 13/24] Working through implementation of canopy model with user docs --- docs/source/users/demography/canopy.md | 435 +++++++++++++++++-------- docs/source/users/demography/crown.md | 49 ++- pyrealm/demography/canopy.py | 31 +- pyrealm/demography/crown.py | 2 +- 4 files changed, 354 insertions(+), 163 deletions(-) diff --git a/docs/source/users/demography/canopy.md b/docs/source/users/demography/canopy.md index 9a9b9901..7068b33d 100644 --- a/docs/source/users/demography/canopy.md +++ b/docs/source/users/demography/canopy.md @@ -30,167 +30,374 @@ notes and initial demonstration code. ::: -The canopy is the representation of the combined crowns across all individual stems -within each cohort in a [community](./community.md). The main use of the canopy is -to calculate the light environment within the canopy, which allows the the productivity -of stems to be estimated across the vertical dimension. - -The key variables in calculating the canopy model are the crown projected area $A_p(z)$ -and leaf projected area $\tilde{A}_{cp}(z)$, which are calculated for a stem -of a given size at a vertical height $z$ using the [crown model](./crown.md). For a -given stem, the difference in projected leaf area values between two heights provides -the actual leaf area between two heights. This then gives an expression for the leaf -area index (LAI, $L$) between two heights $z_i$ and $z_j$: +```{code-cell} ipython3 +from matplotlib import pyplot as plt +import matplotlib as mpl +import numpy as np +import pandas as pd + +from pyrealm.demography.flora import PlantFunctionalType, Flora +from pyrealm.demography.community import Community +from pyrealm.demography.crown import CrownProfile, get_crown_xy +from pyrealm.demography.canopy import Canopy +from pyrealm.demography.t_model_functions import StemAllometry +``` + +The `canopy` module in `pyrealm` is used to calculate a model of the light environment +across all of cohorts within a plant [community](./community.md). Each cohort consists +of a number of identically-sized individuals $n_{ch}$ from a given [plant functional +type (PFT)](./flora.md). Because the individuals are identically sized, they all share +the same [stem allometry](./t_model.md) and the same [crown model](./crown.md). + +## Light extinction for a single stem + +The key variables in determining the light environment for a given +stem are as follows: + +* The projected crown area ${A}_{p}(z)$ sets how crown area accumulates, given + the crown shape, from the top of the stem to the ground. +* The projected leaf area $\tilde{A}_{cp}(z)$ modifies the crown area to allow + for the vertical displacement of crown area by the crown gap fraction. +* The leaf area index $L$ for the PFT is a simple factor that sets the leaf density + of the crown, allowing stems with identical crown area to vary in the density of + actual leaf surface for light capture. Values of $L$ are always expressed as the area + of leaf surface per square meter. +* The extinction coefficient $k$ for a PFT sets how much light is absorbed when passing + through the leaf surface of that PFT. + +For a single stem, the fraction of light absorbed through the entire crown is described +by the Beer-Lambert law: + +$$ +f_{abs} = 1 - e^{-kL} +$$ + +However, to calculate a vertical profile of light extinction through a crown with total +area $A_c$ and maximum stem height $H$, that equation needs to be expanded to calculate +the fraction of $L$ that falls between pairs of vertical heights $z_a > z_b$. The actual +area amount of leaf area $A_l$ for an individual stem falling between those two heights +is simply the diffence in projected leaf area between the two heights: $$ - L_{i,j} = \frac{L\left(\tilde{A}_{cp}(z_j) - \tilde{A}_{cp}(z_i)\right)}{A}, +A_{l[a,b]} = \tilde{A}_{cp}(z_a) - \tilde{A}_{cp}(z_b) $$ -where $z_i > z_j$ and $A$ is the area occupied by the community. The value of $L_{i,j}$ -is then the leaf area index per square metre for the stem between the two heights. The -light extinction between $z_i$ and $z_j$ then follows the Beer Lambert law, -with the fraction of light absorbed $f_{abs}$ calculated as: +Given that information, the calculation of $f_{abs}$ becomes: $$ -f_{abs[i,j]} = 1 - e^{-kL_{i,j}} +f_{abs[a,b]} = 1 - e^{\left(-k\dfrac{L A_{l[a,b]}}{A_c}\right)} $$ -given the extinction coefficient $k$ from the plant functional type of the stem. +When $z_a = H$ and $z_b=0$, then $A_{l[a,b]} = A_c$ and hence simplifies to +the original equation. -The canopy model uses the perfect plasticity approximation (PPA) {cite}`purves:2008a`, -which assumes that plants are always able to plastically arrange their crown within the -broader canopy of the community to maximise their crown size and fill the available -space $A$. When the area $A$ is filled, a new lower canopy layer is formed until all -of the individual crown area has been distributed across within the canopy. +The code below creates a simple example community containing a single cohort containing +a single stem and then calculates the light extinction profile down through the canopy. ```{code-cell} ipython3 -from matplotlib import pyplot as plt -import numpy as np -import pandas as pd +# Create a simple community with a single stem +simple_pft = PlantFunctionalType(name="defaults", m=2, n=2) +simple_flora = Flora([simple_pft]) +stem_dbh = np.array([0.5]) +simple_stem = StemAllometry(stem_traits=simple_flora, at_dbh=stem_dbh) + +# The total area is exactly the crown area +total_area = simple_stem.crown_area[0] + +# Define a simple community +simple_community = Community( + flora=simple_flora, + cell_area=total_area, + cell_id=1, + cohort_dbh_values=stem_dbh, + cohort_n_individuals=np.array([1]), + cohort_pft_names=np.array(["defaults"]), +) -from pyrealm.demography.flora import PlantFunctionalType, Flora -from pyrealm.demography.community import Community -from pyrealm.demography.crown import CrownProfile -from pyrealm.demography.canopy import Canopy +# Get the canopy model for the simple case from the canopy top +# to the ground +hghts = np.linspace(simple_stem.stem_height[0], 0, num=101)[:, None] +simple_canopy = Canopy(community=simple_community, layer_heights=hghts) ``` -## Canopy closure and canopy gap fraction +As a simple check that the calculation across height layers is correct, the canopy +instance returns a vertical light extinction profile. The last value in this profile +should equal the whole canopy $f_{abs}$ calculated using the simple Beer-Lambert +equation and the PFT trait values. -A simple method for finding the first canopy closure height is to find a height $z^*_1$ -at which the sum of crown projected area across all stems $N_s$ in a community equals $A$: +```{code-cell} ipython3 +:scrolled: true -$$ -\sum_1^{N_s}{ A_p(z^*_1)} = A -$$ +print(simple_canopy.extinction_profile[-1]) +``` + +```{code-cell} ipython3 +:scrolled: true + +print(1 - np.exp(-simple_pft.par_ext * simple_pft.lai)) +``` + +The plot below shows: + +1. The shape of crown radius for the stem (solid line) along with the projected leaf + radius (dashed line). The leaf radius does not show the actual expected boundary of + leaf area - which follows the crown - but is useful to visualise the displacement of + leaf area on the same scale as the crown radius. +2. The vertical profile of the actual leaf area $A_{l[a,b]}$ between two height. +3. The resulting light absorption at each height. +4. The light extinction profile through the canopy. + +```{code-cell} ipython3 +fig, (ax1, ax2, ax3, ax4) = plt.subplots(ncols=4, sharey=True, figsize=(12, 6)) + +# Generate plot structures for stem profiles +ch, crown_radius = get_crown_xy( + crown_profile=simple_canopy.crown_profile, + stem_allometry=simple_community.stem_allometry, + attr="crown_radius", + two_sided=False, +)[0] + +ch, projected_leaf_radius = get_crown_xy( + crown_profile=simple_canopy.crown_profile, + stem_allometry=simple_community.stem_allometry, + attr="projected_leaf_radius", + two_sided=False, +)[0] + +ax1.plot(crown_radius, ch, color="red") +ax1.plot(projected_leaf_radius, ch, linestyle="--", color="red") +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.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.set_xlabel("Light absorption fraction (-)") + +# Plot the light extinction profile through the canopy. +ax4.plot(simple_canopy.extinction_profile, hghts, color="red") +ax4.set_xlabel("Cumulative light\nabsorption fraction (-)") +``` + +## Light extinction within a community + +Within a community, the calculations above need to be modified to account for: + +* the number of cohorts $n$, +* the number of individuals $i_{h}$ within each cohort, +* differences in the LAI $L_{h}$ and light extinction coefficients $k_{h}$ between + cohorts, +* scaling LAI to the total area available to the community $A_T$ rather than the cohort + specific crown area $A_h$. -However, the canopy model is modified by a community-level -**canopy gap fraction** ($f_G$) that captures the overall proportion of the canopy area -that is left unfilled by canopy. This gap fraction, capturing processes such as crown -shyness, describes the proportion of open sky visible from the forest floor. This -gives the following definition of the height of canopy layer closure ($z^*_l$) for a -given canopy layer $l = 1, ..., l_m$: +Within the community, each cohort now requires a whole cohort LAI component $L_H$, +which consists of the total leaf area index across individuals divided by the total +community area to give an average leaf area index across the available space: $$ -\sum_1^{N_s}{ A_p(z^*_l)} = l A(1 - f_G) +L_H = \frac{i_h L_h A_h}{A_T} $$ -The set of heights $z^*$ can be found numerically by using a root solver to find -values of $z^*_l$ for $l = 1, ..., l_m$ that satisfy: +The Beer-Lambert equation across the cohorts is then: $$ -\sum_1^{N_s}{ A_p(z^*_l)} - l A(1 - f_G) = 0 +f_{abs} = 1 - e^{\left(\sum\limits_{m=1}^{n}-k_h{[m]} L_{H[m]}\right)} $$ -The total number of layers $l_m$ in a canopy, where the final layer may not be fully -closed, can be found given the total crown area across stems as: +or equivalently $$ -l_m = \left\lceil \frac{\sum_1^{N_s}{A_c}}{ A(1 - f_G)}\right\rceil +f_{abs} = 1 - \prod\limits_{m=1}^{n}e^{-k_{[m]} L_{H[m]}} $$ -+++ +This equation can be adjusted as before to partition light absorption within vertical +layers and the implementation is demonstrated below using a simple community containing +two plant functional types: -## Implementation in `pyrealm` +* a shorter understory tree with a columnar canopy and no crown gaps +* a taller canopy tree with a top heavy canopy and more crown gaps -The {class}`~pyrealm.demography.canopy.Canopy` class automatically finds the canopy -closure heights, given a {class}`~pyrealm.demography.community.Community` instance -and the required canopy gap fraction. +and then three cohorts: -The code below creates a simple community: +* 7 saplings of the short PFT +* 3 larger stems of the short PFT +* 1 large stems of tall PFT ```{code-cell} ipython3 -# Two PFTs -# - a shorter understory tree with a columnar canopy and no crown gaps -# - a taller canopy tree with a top heavy canopy and more crown gaps - +# Define PFTs short_pft = PlantFunctionalType( - name="short", h_max=15, m=1.5, n=1.5, f_g=0, ca_ratio=380 + name="short", + h_max=15, + m=1.5, + n=1.5, + f_g=0, + ca_ratio=380, + lai=4, ) tall_pft = PlantFunctionalType(name="tall", h_max=30, m=3, n=1.5, f_g=0.2, ca_ratio=500) # Create the flora flora = Flora([short_pft, tall_pft]) -# Create a simply community with three cohorts -# - 15 saplings of the short PFT -# - 5 larger stems of the short PFT -# - 2 large stems of tall PFT - +# Define a simply community with three cohorts 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([15, 5, 2]), + cohort_n_individuals=np.array([7, 3, 1]), cohort_pft_names=np.array(["short", "short", "tall"]), ) + +# Calculate the canopy profile across vertical heights +hghts = np.linspace(community.stem_allometry.stem_height.max(), 0, num=101)[:, None] +canopy = Canopy(community=community, canopy_gap_fraction=0, layer_heights=hghts) ``` -We can then look at the expected allometries for the stems in each cohort: +As before, we can verify that the cumulative light extinction at the bottom of the +vertical profile is equal to the expected value across the whole community. ```{code-cell} ipython3 -print("H = ", community.stem_allometry.stem_height) -print("Ac = ", community.stem_allometry.crown_area) +# Calculate L_h for each cohort +cohort_lai = ( + community.cohort_data["n_individuals"] + * community.stem_traits.lai + * community.stem_allometry.crown_area +) / community.cell_area + +# Calculate 1 - e ^ -k L +print(1 - np.exp(np.sum(-community.stem_traits.par_ext * cohort_lai))) ``` -We can now calculate the canopy model for the community: +```{code-cell} ipython3 +print(canopy.extinction_profile[-1]) +``` ```{code-cell} ipython3 -hghts = np.linspace(26, 0, num=251)[:, None] -canopy = Canopy( - community=community, canopy_gap_fraction=2 / 32, layer_heights=hghts, fit_ppa=False +cols = ["r", "b", "g"] + +mpl.rcParams["axes.prop_cycle"] = mpl.cycler(color=cols) + +fig, (ax1, ax2, ax3, ax4) = plt.subplots(ncols=4, sharey=True, figsize=(12, 6)) + +# Generate plot structures for stem profiles +crown_profile = get_crown_xy( + crown_profile=canopy.crown_profile, + stem_allometry=community.stem_allometry, + attr="crown_radius", + two_sided=False, +) + +leaf_profile = get_crown_xy( + crown_profile=canopy.crown_profile, + stem_allometry=community.stem_allometry, + attr="projected_leaf_radius", + two_sided=False, ) + +for (stem_rh, stem_cr), (stem_lh, stem_plr), col in zip( + crown_profile, leaf_profile, cols +): + ax1.plot(stem_cr, stem_rh, color=col) + ax1.plot(stem_plr, stem_lh, linestyle="--", color=col) + +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.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.set_xlabel("Light absorption fraction (-)") + +# Plot the light extinction profile through the canopy. +ax4.plot(canopy.extinction_profile, hghts, color="grey") +ax4.set_xlabel("Cumulative light\nabsorption fraction (-)") ``` +## Canopy closure and canopy gap fraction + +In addition to calculating profiles from a provided sequence of vertical heights, the +canopy model also implements the calculation of canopy layers, following the perfect +plasticity approximation (PPA) {cite}`purves:2008a`. This model divides the vertical +structure of the canopy into discrete closed layers. The model assumes that plants are +always able to plastically arrange their crown within the broader canopy of the +community to maximise their crown size and fill the available space $A$. When the area +$A$ is filled, a new lower canopy layer is formed until all of the individual crown area +has been distributed across within the canopy. + +A simple method for finding the first canopy closure height is to find a height $z^*_1$ +at which the sum of crown projected area across all stems $N_s$ in a community equals +$A$: + +$$ +\sum_1^{N_s}{ A_p(z^*_1)} = A +$$ + +However, the canopy model also allows for modification by a community-level **canopy gap +fraction** ($f_G$) that captures the overall proportion of the canopy area that is left +unfilled by canopy. This gap fraction, capturing processes such as crown shyness, +describes the proportion of open sky visible from the forest floor. This gives the +following definition of the height of canopy layer closure ($z^*_l$) for a given canopy +layer $l = 1, ..., l_m$: + +$$ +\sum_1^{N_s}{ A_p(z^*_l)} = l A(1 - f_G) +$$ + +The set of heights $z^*$ can be found numerically by using a root solver to find values +of $z^*_l$ for $l = 1, ..., l_m$ that satisfy: + +$$ +\sum_1^{N_s}{ A_p(z^*_l)} - l A(1 - f_G) = 0 +$$ + +The total number of layers $l_m$ in a canopy, where the final layer may not be fully +closed, can be found given the total crown area across stems as: + +$$ +l_m = \left\lceil \frac{\sum_1^{N_s}{A_c}}{ A(1 - f_G)}\right\rceil +$$ + ```{code-cell} ipython3 -fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(ncols=5, sharey=True, figsize=(12, 6)) +canopy = Canopy(community=community, canopy_gap_fraction=2 / 32, fit_ppa=True) +``` -ax1.plot(canopy.fapar_profile, hghts) +```{code-cell} ipython3 +canopy.layer_heights +``` -ax2.plot(np.cumsum(canopy.fapar_profile), hghts) -ax3.plot(canopy.layer_stem_leaf_area, hghts) -ax4.plot(canopy.crown_profile.projected_leaf_radius, hghts) +## Implementation in `pyrealm` -ax5.set_prop_cycle("color", ["r", "b", "g"]) -ax5.plot(canopy.crown_profile.crown_radius, hghts) -ax5.plot(canopy.crown_profile.projected_leaf_radius, hghts, linestyle="--") +The {class}`~pyrealm.demography.canopy.Canopy` class automatically finds the canopy +closure heights, given a {class}`~pyrealm.demography.community.Community` instance and +the required canopy gap fraction. -# ax1.plot(canopy.crown_profile.projected_leaf_radius, hghts) -# ax1.plot(canopy.crown_profile.projected_crown_radius, hghts) -# ax1.set_aspect(1) -# plt.xscale('log') +```{code-cell} ipython3 +canopy.crown_profile.projected_leaf_area[:, 1] ``` ```{code-cell} ipython3 -canopy.stem_fapar.sum() +canopy.layer_stem_leaf_area[:, 1] ``` -We can then look at three key properties of the canopy model: the layer closure -heights ($z^*_l$) and the projected crown areas and leaf areas at each of those -heights for each stem in the three cohorts. +```{code-cell} ipython3 +plt.plot(canopy.layer_stem_leaf_area[:, 1], "ro") +``` + +We can then look at three key properties of the canopy model: the layer closure heights +($z^*_l$) and the projected crown areas and leaf areas at each of those heights for each +stem in the three cohorts. -There are four canopy layers, with the top two very close together because of the -large crown area in the two stems in the cohort of `tall` trees. +There are four canopy layers, with the top two very close together because of the large +crown area in the two stems in the cohort of `tall` trees. ```{code-cell} ipython3 canopy.layer_heights @@ -261,7 +468,9 @@ total area at each closed layer height is omitting the community gap fraction. fig, ax = plt.subplots(ncols=1) # Calculate the crown area at which each canopy layer closes. -closure_areas = np.arange(1, canopy.n_layers + 1) * canopy.crown_area_per_layer +closure_areas = ( + np.arange(1, len(canopy.layer_heights) + 1) * community.cell_area * (30 / 32) +) # Add lines showing the canopy closure heights and closure areas. for val in canopy.layer_heights: @@ -305,49 +514,15 @@ ax.set_xlabel("Community-wide projected area (m2)") ax.legend(frameon=False) ``` -### Light transmission through the canopy - -Now we can use the leaf area by layer and the Beer-Lambert equation to calculate light -attenuation through the canopy layers. - -$f_{abs} = 1 - e ^ {-kL}$, - -where $k$ is the light extinction coefficient ($k$) and $L$ is the leaf area index -(LAI). The LAI can be calculated for each stem and layer: - -```{code-cell} ipython3 -# LAI = Acp_within_layer / canopy_area -# print(LAI) -``` - -This can be used to calculate the LAI of individual stems but also the LAI of each layer -in the canopy: - ```{code-cell} ipython3 -# LAI_stem = LAI.sum(axis=0) -# LAI_layer = LAI.sum(axis=1) +ppfd = 1000 -# print("LAI stem = ", LAI_stem) -# print("LAI layer = ", LAI_layer) +ppfd_layer = -np.diff(ppfd * np.cumprod(canopy.layer_f_trans), prepend=ppfd) +ppfd_layer ``` -The layer LAI values can now be used to calculate the light transmission of each layer and -hence the cumulative light extinction profile through the canopy. - -```{code-cell} ipython3 -# f_abs = 1 - np.exp(-pft.traits.par_ext * LAI_layer) -# ext = np.cumprod(f_abs) - -# print("f_abs = ", f_abs) -# print("extinction = ", ext) -``` - -One issue that needs to be resolved is that the T Model implementation in `pyrealm` -follows the original implementation of the T Model in having LAI as a fixed trait of -a given plant functional type, so is constant for all stems of that PFT. - ```{code-cell} ipython3 -# print("f_abs = ", (1 - np.exp(-pft.traits.par_ext * pft.traits.lai))) +(1 - canopy.layer_cohort_f_trans) * ppfd_layer[:, None] ``` ## Things to worry about later diff --git a/docs/source/users/demography/crown.md b/docs/source/users/demography/crown.md index 2653a5c0..69e64018 100644 --- a/docs/source/users/demography/crown.md +++ b/docs/source/users/demography/crown.md @@ -5,10 +5,21 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 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 --- # The tree crown model @@ -20,7 +31,7 @@ notes and initial demonstration code. ::: -```{code-cell} +```{code-cell} ipython3 from matplotlib import pyplot as plt from matplotlib.lines import Line2D from matplotlib.patches import Polygon, Patch @@ -149,7 +160,7 @@ The {class}`~pyrealm.demography.flora.PlantFunctionalType` class is typically used to set specific PFTs, but the functions to calculate $q_m$ and $p_{zm}$ are used directly below to provides a demonstration of the impacts of each trait. -```{code-cell} +```{code-cell} ipython3 # Set a range of values for m and n traits m = n = np.arange(1.0, 5, 0.1) @@ -158,7 +169,7 @@ q_m = calculate_crown_q_m(m=m, n=n[:, None]) z_max_prop = calculate_crown_z_max_proportion(m=m, n=n[:, None]) ``` -```{code-cell} +```{code-cell} ipython3 fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10.9, 4)) # Plot q_m as a function of m and n @@ -194,7 +205,7 @@ profiles for PFTs. It requires: The code below creates a set of PFTS with differing crown trait values and then creates a `Flora` object using the PFTs. -```{code-cell} +```{code-cell} ipython3 # A PFT with a small crown area and equal m and n values narrow_pft = PlantFunctionalType(name="narrow", h_max=20, m=1.5, n=1.5, ca_ratio=20) # A PFT with an intermediate crown area and m < n @@ -209,7 +220,7 @@ flora The Flora object can also be used to show a table of canopy variables: -```{code-cell} +```{code-cell} ipython3 # TODO - add a Flora.to_pandas() method flora_data = pd.DataFrame({k: getattr(flora, k) for k in flora.trait_attrs}) flora_data[["name", "ca_ratio", "m", "n", "f_g", "q_m", "z_max_prop"]] @@ -221,7 +232,7 @@ The T Model requires DBH to define stem size - here the is used to back-calculate the required DBH values to give three stems with similar heights that are near the maximum height for each PFT. -```{code-cell} +```{code-cell} ipython3 # Generate the expected stem allometries at similar heights for each PFT stem_height = np.array([19, 17, 15]) stem_dbh = calculate_dbh_from_height( @@ -230,14 +241,14 @@ stem_dbh = calculate_dbh_from_height( stem_dbh ``` -```{code-cell} +```{code-cell} ipython3 # Calculate the stem allometries allometry = StemAllometry(stem_traits=flora, at_dbh=stem_dbh) ``` We can again use {mod}`pandas` to get a table of those allometric predictions: -```{code-cell} +```{code-cell} ipython3 pd.DataFrame({k: getattr(allometry, k) for k in allometry.allometry_attrs}) ``` @@ -247,7 +258,7 @@ that is with a shape `(N, 1)`. We can then calculate the crown profiles. -```{code-cell} +```{code-cell} ipython3 # Create a set of vertical heights as a column array. z = np.linspace(-1, 20.0, num=211)[:, None] @@ -263,7 +274,7 @@ above calculated at each height $z$: * The projected crown area * The projected leaf area -```{code-cell} +```{code-cell} ipython3 crown_profiles ``` @@ -292,7 +303,7 @@ restricted to actual valid heights for that stem and is demonstrated in the ::: -```{code-cell} +```{code-cell} ipython3 fig, ax = plt.subplots(ncols=1) # Find the maximum of the actual and relative maximum crown widths @@ -334,7 +345,7 @@ We can also use the `CanopyProfile` class with a single row of heights to calcul the crown profile at the expected $z_{max}$ and show that this matches the expected crown area from the T Model allometry. -```{code-cell} +```{code-cell} ipython3 # Calculate the crown profile across those heights for each PFT z_max = flora.z_max_prop * stem_height profile_at_zmax = CrownProfile(stem_traits=flora, stem_allometry=allometry, z=z_max) @@ -350,7 +361,7 @@ using the PFTs defined above because they have very different crown areas, so th below generates new profiles for a new set of PFTs that have similar crown area ratios but different shapes and gap fractions. -```{code-cell} +```{code-cell} ipython3 no_gaps_pft = PlantFunctionalType( name="no_gaps", h_max=20, m=1.5, n=1.5, f_g=0, ca_ratio=380 ) @@ -386,7 +397,7 @@ lines) change with height along the stem. lines are identical, but as `f_g` increases, more of the leaf area is displaced down within the crown. -```{code-cell} +```{code-cell} ipython3 fig, ax = plt.subplots(ncols=1) for pft_idx, offset, colour in zip((0, 1, 2), (0, 5, 10), ("r", "g", "b")): @@ -406,7 +417,7 @@ We can also generate predictions for a single PFT with varying crown gap fractio the plot below, note that all leaf area is above $z_{max}$ when $f_g=1$ and all leaf area is *below* -```{code-cell} +```{code-cell} ipython3 fig, ax = plt.subplots(ncols=1) # Loop over f_g values @@ -469,7 +480,7 @@ projected crown radius and projected leaf radius. These last two variables do no have direct computational use - the cumulative projected area is what matters - but allow the projected variables to be visualised at the same scale as the crown radius. -```{code-cell} +```{code-cell} ipython3 # Set stem offsets for separating stems along the x axis stem_offsets = np.array([0, 6, 12]) @@ -498,7 +509,7 @@ projected_leaf_radius_xy = get_crown_xy( ) ``` -```{code-cell} +```{code-cell} ipython3 fig, ax = plt.subplots() # Bundle the three plotting structures and loop over the three stems. @@ -521,3 +532,7 @@ plt.legend( bbox_to_anchor=(0.5, 1.15), ) ``` + +```{code-cell} ipython3 + +``` diff --git a/pyrealm/demography/canopy.py b/pyrealm/demography/canopy.py index 1d57119d..7271fb93 100644 --- a/pyrealm/demography/canopy.py +++ b/pyrealm/demography/canopy.py @@ -277,37 +277,38 @@ def _calculate_canopy(self, community: Community) -> None: self.layer_stem_leaf_area = np.diff( self.crown_profile.projected_leaf_area, axis=0, prepend=0 ) - self.layer_cohort_leaf_area = ( - self.layer_stem_leaf_area * community.cohort_data["n_individuals"] - ) - # Calculate the leaf area index per layer per cohort, using the stem + # 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 community area. self.layer_cohort_lai = ( - self.layer_cohort_leaf_area * community.stem_traits.lai + self.layer_stem_leaf_area + * community.cohort_data["n_individuals"] + * community.stem_traits.lai ) / community.cell_area - # Calculate the Beer-Lambert light extinction per layer and cohort - self.layer_cohort_f_abs = 1 - np.exp( + # 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 ) + # Aggregate across cohorts into a layer wide transimissivity + self.layer_f_trans = self.layer_cohort_f_trans.prod(axis=1) # Calculate the canopy wide light extinction per layer - self.layer_canopy_f_abs = self.layer_cohort_f_abs.sum(axis=1) + self.layer_f_abs = 1 - self.layer_f_trans # Calculate cumulative light extinction across the canopy - self.canopy_extinction_profile = np.cumprod(self.layer_canopy_f_abs) + 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_canopy_f_abs), + 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_abs * self.fapar_profile[:, None] - ) / self.layer_cohort_f_abs.sum(axis=1)[:, None] + # # 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] diff --git a/pyrealm/demography/crown.py b/pyrealm/demography/crown.py index 9aaaee3f..7d5f7665 100644 --- a/pyrealm/demography/crown.py +++ b/pyrealm/demography/crown.py @@ -414,7 +414,7 @@ def get_crown_xy( crown_profile: CrownProfile, stem_allometry: StemAllometry, attr: str, - stem_offsets: NDArray[np.float32] | None, + stem_offsets: NDArray[np.float32] | None = None, two_sided: bool = True, as_xy: bool = False, ) -> list[tuple[NDArray, NDArray]] | list[NDArray]: From 5a56a2264ef7cfa7328da9bac23166175b0219a8 Mon Sep 17 00:00:00 2001 From: David Orme Date: Thu, 10 Oct 2024 11:08:29 +0100 Subject: [PATCH 14/24] Docs update and some bug catches --- docs/source/users/demography/canopy.md | 157 ++++++++----------------- pyrealm/demography/canopy.py | 16 ++- pyrealm/demography/crown.py | 6 +- 3 files changed, 61 insertions(+), 118 deletions(-) diff --git a/docs/source/users/demography/canopy.md b/docs/source/users/demography/canopy.md index 7068b33d..2a1fa595 100644 --- a/docs/source/users/demography/canopy.md +++ b/docs/source/users/demography/canopy.md @@ -45,23 +45,28 @@ from pyrealm.demography.t_model_functions import StemAllometry The `canopy` module in `pyrealm` is used to calculate a model of the light environment across all of cohorts within a plant [community](./community.md). Each cohort consists -of a number of identically-sized individuals $n_{ch}$ from a given [plant functional -type (PFT)](./flora.md). Because the individuals are identically sized, they all share -the same [stem allometry](./t_model.md) and the same [crown model](./crown.md). +of: + +* a number of identically-sized individuals, +* of the same [plant functional type (PFT)](./flora.md), +* that share the same [stem allometry](./t_model.md) and [crown model](./crown.md). + +The purpose of the `canopy` module is to estimate how light is absorbed down through the +canopy and allow the absorption of incoming light at different heights in the canopy to +be partitioned across stems within each cohort. ## Light extinction for a single stem -The key variables in determining the light environment for a given -stem are as follows: - -* The projected crown area ${A}_{p}(z)$ sets how crown area accumulates, given - the crown shape, from the top of the stem to the ground. -* The projected leaf area $\tilde{A}_{cp}(z)$ modifies the crown area to allow - for the vertical displacement of crown area by the crown gap fraction. -* The leaf area index $L$ for the PFT is a simple factor that sets the leaf density - of the crown, allowing stems with identical crown area to vary in the density of - actual leaf surface for light capture. Values of $L$ are always expressed as the area - of leaf surface per square meter. +The key variables in determining the light environment for a given stem are as follows: + +* The projected crown area ${A}_{p}(z)$ sets how crown area accumulates, given the crown + shape, from the top of the stem to the ground. +* The projected leaf area $\tilde{A}_{cp}(z)$ modifies the crown area to allow for the + vertical displacement of crown area by the crown gap fraction. +* The leaf area index $L$ for the PFT is a simple factor that sets the leaf density of + the crown, allowing stems with identical crown area to vary in the density of actual + leaf surface for light capture. Values of $L$ are always expressed as the area of leaf + surface per square meter. * The extinction coefficient $k$ for a PFT sets how much light is absorbed when passing through the leaf surface of that PFT. @@ -88,8 +93,8 @@ $$ f_{abs[a,b]} = 1 - e^{\left(-k\dfrac{L A_{l[a,b]}}{A_c}\right)} $$ -When $z_a = H$ and $z_b=0$, then $A_{l[a,b]} = A_c$ and hence simplifies to -the original equation. +When $z_a = H$ and $z_b=0$, then $A_{l[a,b]} = A_c$ and hence simplifies to the original +equation. The code below creates a simple example community containing a single cohort containing a single stem and then calculates the light extinction profile down through the canopy. @@ -117,7 +122,10 @@ simple_community = Community( # Get the canopy model for the simple case from the canopy top # to the ground hghts = np.linspace(simple_stem.stem_height[0], 0, num=101)[:, None] -simple_canopy = Canopy(community=simple_community, layer_heights=hghts) +simple_canopy = Canopy( + community=simple_community, + layer_heights=hghts, +) ``` As a simple check that the calculation across height layers is correct, the canopy @@ -126,14 +134,10 @@ should equal the whole canopy $f_{abs}$ calculated using the simple Beer-Lambert equation and the PFT trait values. ```{code-cell} ipython3 -:scrolled: true - print(simple_canopy.extinction_profile[-1]) ``` ```{code-cell} ipython3 -:scrolled: true - print(1 - np.exp(-simple_pft.par_ext * simple_pft.lai)) ``` @@ -205,13 +209,8 @@ $$ The Beer-Lambert equation across the cohorts is then: $$ -f_{abs} = 1 - e^{\left(\sum\limits_{m=1}^{n}-k_h{[m]} L_{H[m]}\right)} -$$ - -or equivalently - -$$ -f_{abs} = 1 - \prod\limits_{m=1}^{n}e^{-k_{[m]} L_{H[m]}} +f_{abs} = 1 - e^{\left(\sum\limits_{m=1}^{n}-k_{[m]} L_{H[m]}\right)} + = 1 - \prod\limits_{m=1}^{n}e^{-k_{[m]} L_{H[m]}} $$ This equation can be adjusted as before to partition light absorption within vertical @@ -255,7 +254,7 @@ community = Community( # Calculate the canopy profile across vertical heights hghts = np.linspace(community.stem_allometry.stem_height.max(), 0, num=101)[:, None] -canopy = Canopy(community=community, canopy_gap_fraction=0, layer_heights=hghts) +canopy = Canopy(community=community, layer_heights=hghts) ``` As before, we can verify that the cumulative light extinction at the bottom of the @@ -367,95 +366,31 @@ l_m = \left\lceil \frac{\sum_1^{N_s}{A_c}}{ A(1 - f_G)}\right\rceil $$ ```{code-cell} ipython3 -canopy = Canopy(community=community, canopy_gap_fraction=2 / 32, fit_ppa=True) -``` - -```{code-cell} ipython3 -canopy.layer_heights -``` - -## Implementation in `pyrealm` - -The {class}`~pyrealm.demography.canopy.Canopy` class automatically finds the canopy -closure heights, given a {class}`~pyrealm.demography.community.Community` instance and -the required canopy gap fraction. - -```{code-cell} ipython3 -canopy.crown_profile.projected_leaf_area[:, 1] -``` - -```{code-cell} ipython3 -canopy.layer_stem_leaf_area[:, 1] -``` - -```{code-cell} ipython3 -plt.plot(canopy.layer_stem_leaf_area[:, 1], "ro") -``` - -We can then look at three key properties of the canopy model: the layer closure heights -($z^*_l$) and the projected crown areas and leaf areas at each of those heights for each -stem in the three cohorts. - -There are four canopy layers, with the top two very close together because of the large -crown area in the two stems in the cohort of `tall` trees. - -```{code-cell} ipython3 -canopy.layer_heights -``` - -The `stem_crown_area` attribute then provides the crown area of each stem found in each -layer. - -```{code-cell} ipython3 -canopy.stem_crown_area +canopy_ppa = Canopy(community=community, canopy_gap_fraction=0, fit_ppa=True) ``` -Given the canopy gap fraction, the available crown area per layer is 30 m2, so -the first two layers are taken up entirely by the two stems in the cohort of large -trees. We can confirm that the calculation is correct by calculating the total crown area -across the cohorts at each height: - -```{code-cell} ipython3 -np.sum(canopy.stem_crown_area * community.cohort_data["n_individuals"], axis=1) -``` - -Those are equal to the layer closure areas of 30, 60 and 90 m2 and the last layer does -not quite close. The slight numerical differences result from the precision of the root -solver for finding $z^*_l$ and this can be adjusted by using the `layer_tolerance` -argument to the `Canopy` class - -The projected leaf area per stem is reported in the `stem_leaf_area` attribute. This is -identical to the projected crown area for the first two cohorts because the crown gap -fraction $f_g$ is zero for this PFT. The projected leaf area is however displaced -towards the ground in the last cohort, because the `tall` PFT has a large gap fraction. - ```{code-cell} ipython3 -canopy.stem_leaf_area +canopy_ppa.layer_heights ``` ### Visualizing layer closure heights and areas -We can use the {class}`~pyrealm.demography.crown.CrownProfile` class to calculate a -community crown and leaf area profile across a range of height values. For each height, +We can use the crown profile calculated for the previous canopy model to calculate a +whole community crown and leaf area profile for the community. For each height, we calculate the sum of the product of stem projected area and the number of individuals in each cohort. ```{code-cell} ipython3 -# Set of vertical height to calculate crown profiles -at_z = np.linspace(0, 26, num=261)[:, None] - -# Calculate the crown profile for the stem for each cohort -crown_profiles = CrownProfile( - stem_traits=community.stem_traits, stem_allometry=community.stem_allometry, z=at_z -) - # Calculate the total projected crown area across the community at each height community_crown_area = np.nansum( - crown_profiles.projected_crown_area * community.cohort_data["n_individuals"], axis=1 + canopy.crown_profile.projected_crown_area * community.cohort_data["n_individuals"], + axis=1, ) + # Do the same for the projected leaf area community_leaf_area = np.nansum( - crown_profiles.projected_leaf_area * community.cohort_data["n_individuals"], axis=1 + canopy.crown_profile.projected_leaf_area * community.cohort_data["n_individuals"], + axis=1, ) ``` @@ -469,21 +404,21 @@ fig, ax = plt.subplots(ncols=1) # Calculate the crown area at which each canopy layer closes. closure_areas = ( - np.arange(1, len(canopy.layer_heights) + 1) * community.cell_area * (30 / 32) + np.arange(1, len(canopy_ppa.layer_heights) + 1) * canopy.filled_community_area ) # Add lines showing the canopy closure heights and closure areas. -for val in canopy.layer_heights: +for val in canopy_ppa.layer_heights: ax.axhline(val, color="red", linewidth=0.5, zorder=0) for val in closure_areas: ax.axvline(val, color="red", linewidth=0.5, zorder=0) # Show the community projected crown area profile -ax.plot(community_crown_area, at_z, zorder=1, label="Crown area") +ax.plot(community_crown_area, canopy.layer_heights, zorder=1, label="Crown area") ax.plot( community_leaf_area, - at_z, + canopy.layer_heights, zorder=1, linestyle="--", color="black", @@ -497,9 +432,9 @@ ax_rhs = ax.twinx() ax_rhs.set_ylim(ax.get_ylim()) z_star_labels = [ f"$z^*_{l + 1} = {val:.2f}$" - for l, val in enumerate(np.nditer(canopy.layer_heights)) + for l, val in enumerate(np.nditer(canopy_ppa.layer_heights)) ] -ax_rhs.set_yticks(canopy.layer_heights.flatten()) +ax_rhs.set_yticks(canopy_ppa.layer_heights.flatten()) ax_rhs.set_yticklabels(z_star_labels) # Add cumulative canopy area at top @@ -514,15 +449,17 @@ ax.set_xlabel("Community-wide projected area (m2)") ax.legend(frameon=False) ``` +## Light allocation + ```{code-cell} ipython3 ppfd = 1000 -ppfd_layer = -np.diff(ppfd * np.cumprod(canopy.layer_f_trans), prepend=ppfd) +ppfd_layer = -np.diff(ppfd * np.cumprod(canopy_ppa.layer_f_trans), prepend=ppfd) ppfd_layer ``` ```{code-cell} ipython3 -(1 - canopy.layer_cohort_f_trans) * ppfd_layer[:, None] +(1 - canopy_ppa.layer_cohort_f_trans) * ppfd_layer[:, None] ``` ## Things to worry about later diff --git a/pyrealm/demography/canopy.py b/pyrealm/demography/canopy.py index 7271fb93..dd38854c 100644 --- a/pyrealm/demography/canopy.py +++ b/pyrealm/demography/canopy.py @@ -193,7 +193,7 @@ def __init__( community: Community, layer_heights: NDArray[np.float32] | None = None, fit_ppa: bool = False, - canopy_gap_fraction: float = 0.05, + canopy_gap_fraction: float = 0, solver_tolerance: float = 0.001, ) -> None: # Store required init vars @@ -204,9 +204,9 @@ def __init__( # Define class attributes self.total_community_crown_area: float - """Total crown area across individuals in the community (metres 2).""" + """Total crown area across all individuals in the community (m2).""" self.max_stem_height: float - """Maximum height of any individual in the community (metres).""" + """Maximum height of any individual in the community (m).""" self.n_layers: int """Total number of canopy layers.""" self.n_cohorts: int @@ -231,7 +231,8 @@ def __init__( """The canopy-wide fAPAR profile by layer.""" self.stem_fapar: NDArray[np.float32] """The fAPAR for individual stems by layer.""" - + self.filled_community_area: float + """The area filled by crown.""" # Check operating mode if fit_ppa ^ (layer_heights is None): raise ValueError("Either set fit_ppa=True or provide layer heights.") @@ -239,6 +240,9 @@ def __init__( # Set simple attributes self.max_stem_height = community.stem_allometry.stem_height.max() self.n_cohorts = community.number_of_cohorts + self.filled_community_area = community.cell_area * ( + 1 - self.canopy_gap_fraction + ) # Populate layer heights if layer_heights is not None: @@ -280,12 +284,12 @@ def _calculate_canopy(self, community: Community) -> None: # 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 community area. + # the available community area. self.layer_cohort_lai = ( self.layer_stem_leaf_area * community.cohort_data["n_individuals"] * community.stem_traits.lai - ) / community.cell_area + ) / 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( diff --git a/pyrealm/demography/crown.py b/pyrealm/demography/crown.py index 7d5f7665..ac9dac50 100644 --- a/pyrealm/demography/crown.py +++ b/pyrealm/demography/crown.py @@ -77,7 +77,9 @@ def _validate_z_qz_args( # Now test q_z congruence with z if provided if q_z is not None: - if ((z.size == 1) or (z.ndim == 1)) and (q_z.shape != stem_shape): + if q_z.shape == z.shape: + pass + elif ((z.size == 1) or (z.ndim == 1)) and (q_z.shape != stem_shape): raise ValueError( f"The q_z argument (shape: {q_z.shape}) is not a row array " f"matching stem properties (shape: {stem_shape})" @@ -198,7 +200,7 @@ def calculate_stem_projected_crown_area_at_z( crown_area: Crown area of each stem stem_height: Stem height of each stem q_m: Canopy shape parameter ``q_m``` for each stem - z_max: Height of maximum crown radous for each stem + z_max: Height of maximum crown radi `us for each stem validate: Boolean flag to suppress argument validation. """ From 42d8106dfc66ca14abbe1f6cad6251854abccfc4 Mon Sep 17 00:00:00 2001 From: David Orme Date: Thu, 10 Oct 2024 15:58:57 +0100 Subject: [PATCH 15/24] Doc build fixes --- pyrealm/demography/crown.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyrealm/demography/crown.py b/pyrealm/demography/crown.py index ac9dac50..00f0019b 100644 --- a/pyrealm/demography/crown.py +++ b/pyrealm/demography/crown.py @@ -200,7 +200,7 @@ def calculate_stem_projected_crown_area_at_z( crown_area: Crown area of each stem stem_height: Stem height of each stem q_m: Canopy shape parameter ``q_m``` for each stem - z_max: Height of maximum crown radi `us for each stem + z_max: Height of maximum crown radius for each stem validate: Boolean flag to suppress argument validation. """ From a57287ce41c8938079a82e9ed917d7c4df49f612 Mon Sep 17 00:00:00 2001 From: David Orme Date: Fri, 11 Oct 2024 13:40:03 +0100 Subject: [PATCH 16/24] Updating canopy docs --- docs/source/users/demography/canopy.md | 176 +++++++++++++++++++++---- 1 file changed, 151 insertions(+), 25 deletions(-) diff --git a/docs/source/users/demography/canopy.md b/docs/source/users/demography/canopy.md index 2a1fa595..6894f11c 100644 --- a/docs/source/users/demography/canopy.md +++ b/docs/source/users/demography/canopy.md @@ -19,6 +19,8 @@ language_info: nbconvert_exporter: python pygments_lexer: ipython3 version: 3.11.9 +settings: + output_matplotlib_strings: remove --- # The canopy model @@ -152,6 +154,8 @@ The plot below shows: 4. The light extinction profile through the canopy. ```{code-cell} ipython3 +:tags: [hide-input] + fig, (ax1, ax2, ax3, ax4) = plt.subplots(ncols=4, sharey=True, figsize=(12, 6)) # Generate plot structures for stem profiles @@ -233,11 +237,14 @@ short_pft = PlantFunctionalType( h_max=15, m=1.5, n=1.5, + par_ext=0.4, f_g=0, ca_ratio=380, lai=4, ) -tall_pft = PlantFunctionalType(name="tall", h_max=30, m=3, n=1.5, f_g=0.2, ca_ratio=500) +tall_pft = PlantFunctionalType( + name="tall", h_max=30, m=3, n=1.5, par_ext=0.6, f_g=0.2, ca_ratio=500 +) # Create the flora flora = Flora([short_pft, tall_pft]) @@ -318,11 +325,20 @@ ax3.set_xlabel("Light absorption fraction (-)") # Plot the light extinction profile through the canopy. ax4.plot(canopy.extinction_profile, hghts, color="grey") -ax4.set_xlabel("Cumulative light\nabsorption fraction (-)") +_ = ax4.set_xlabel("Cumulative light\nabsorption fraction (-)") ``` ## Canopy closure and canopy gap fraction +:::{admonition} TODO + +Need to work out how to include the gap fraction in the calculation of light extinction +because at the moment, the gap fraction in the PPA calculates the layer closure heights +accounting for that, but the LAI is not accounting for it so there is no shift in the +light extinction profile. + +::: + In addition to calculating profiles from a provided sequence of vertical heights, the canopy model also implements the calculation of canopy layers, following the perfect plasticity approximation (PPA) {cite}`purves:2008a`. This model divides the vertical @@ -369,10 +385,20 @@ $$ 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 +layers close: + ```{code-cell} ipython3 canopy_ppa.layer_heights ``` +And the final value in the canopy extinction profile still matches the expectation from +above: + +```{code-cell} ipython3 +print(canopy_ppa.extinction_profile[-1]) +``` + ### Visualizing layer closure heights and areas We can use the crown profile calculated for the previous canopy model to calculate a @@ -400,23 +426,25 @@ 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 -fig, ax = plt.subplots(ncols=1) +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 ) +# LH plot - projected leaf area with height. + # Add lines showing the canopy closure heights and closure areas. for val in canopy_ppa.layer_heights: - ax.axhline(val, color="red", linewidth=0.5, zorder=0) + ax1.axhline(val, color="red", linewidth=0.5, zorder=0) for val in closure_areas: - ax.axvline(val, color="red", linewidth=0.5, zorder=0) + ax1.axvline(val, color="red", linewidth=0.5, zorder=0) # Show the community projected crown area profile -ax.plot(community_crown_area, canopy.layer_heights, zorder=1, label="Crown area") -ax.plot( +ax1.plot(community_crown_area, canopy.layer_heights, zorder=1, label="Crown area") +ax1.plot( community_leaf_area, canopy.layer_heights, zorder=1, @@ -426,42 +454,140 @@ ax.plot( label="Leaf area", ) +# Add cumulative canopy area at top +ax1_top = ax1.twiny() +ax1_top.set_xlim(ax1.get_xlim()) +area_labels = [f"$A_{l + 1}$ = {z:.1f}" for l, z in enumerate(np.nditer(closure_areas))] +ax1_top.set_xticks(closure_areas) +ax1_top.set_xticklabels(area_labels, rotation=90) + +ax1.set_ylabel("Vertical height ($z$, m)") +ax1.set_xlabel("Community-wide projected area (m2)") +ax1.legend(frameon=False) + +# RH plot - light extinction +for val in canopy_ppa.layer_heights: + ax2.axhline(val, color="red", linewidth=0.5, zorder=0) + +for val in canopy_ppa.extinction_profile: + ax2.axvline(val, color="red", linewidth=0.5, zorder=0) + +ax2.plot(canopy.extinction_profile, hghts) + +ax2_top = ax2.twiny() +ax2_top.set_xlim(ax2.get_xlim()) +extinction_labels = [ + f"$f_{{abs{l + 1}}}$ = {z:.3f}" + for l, z in enumerate(np.nditer(canopy_ppa.extinction_profile)) +] +ax2_top.set_xticks(canopy_ppa.extinction_profile) +ax2_top.set_xticklabels(extinction_labels, rotation=90) + +ax2.set_xlabel("Light extinction (-)") # Add z* values on the righthand axis -ax_rhs = ax.twinx() -ax_rhs.set_ylim(ax.get_ylim()) +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)) ] -ax_rhs.set_yticks(canopy_ppa.layer_heights.flatten()) -ax_rhs.set_yticklabels(z_star_labels) +ax2_rhs.set_yticks(canopy_ppa.layer_heights.flatten()) +_ = ax2_rhs.set_yticklabels(z_star_labels) +``` -# Add cumulative canopy area at top -ax_top = ax.twiny() -ax_top.set_xlim(ax.get_xlim()) -area_labels = [f"$A_{l + 1}$ = {z:.1f}" for l, z in enumerate(np.nditer(closure_areas))] -ax_top.set_xticks(closure_areas) -ax_top.set_xticklabels(area_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. +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}$. -ax.set_ylabel("Vertical height ($z$, m)") -ax.set_xlabel("Community-wide projected area (m2)") -ax.legend(frameon=False) +```{code-cell} ipython3 +print("k = \n", community.stem_traits.par_ext, "\n") +print("L_H = \n", canopy_ppa.layer_cohort_lai) ``` -## Light allocation +```{code-cell} ipython3 +layer_cohort_f_tr = np.exp(-community.stem_traits.par_ext * canopy_ppa.layer_cohort_lai) +print(layer_cohort_f_tr) +``` + +The fraction absorbed $f_{abs} = 1 - f_{tr}$. ```{code-cell} ipython3 -ppfd = 1000 +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)) +``` + +Or alternatively, take the product within layers of the cohort components: + +```{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) +``` -ppfd_layer = -np.diff(ppfd * np.cumprod(canopy_ppa.layer_f_trans), prepend=ppfd) -ppfd_layer +We can also calculate the cumulative fraction of light transmitted through the layers: + +```{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: + +```{code-cell} ipython3 +ppfd = 1000 +ppfd_layer = -np.diff(ppfd * transmission_profile, prepend=ppfd, append=0) +print(ppfd_layer) ``` +The cumulative sum across those quantities accounts for all of the incoming light and +matches the scaling of the extinction profile: + ```{code-cell} ipython3 -(1 - canopy_ppa.layer_cohort_f_trans) * ppfd_layer[:, None] +print(np.cumsum(ppfd_layer)) ``` +:::{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! + +::: + ## Things to worry about later Herbivory - leaf fall (transmission increases, truncate at 0, replace from NSCs) vs leaf From 9d920b1e9002177fe21f7999270484d1a28da90d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 11 Oct 2024 13:59:32 +0000 Subject: [PATCH 17/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- docs/source/api/constants_api.md | 1 + docs/source/api/core_api.md | 1 + docs/source/api/demography_api.md | 1 + docs/source/api/pmodel_api.md | 1 + docs/source/api/splash_api.md | 1 + docs/source/api/tmodel_api.md | 1 + docs/source/development/code_qa_and_typing.md | 1 + docs/source/development/code_testing.md | 1 + docs/source/development/documentation.md | 1 + docs/source/development/github_actions.md | 1 + docs/source/development/overview.md | 1 + .../development/profiling_and_benchmarking.md | 1 + docs/source/development/pyrealm_build_data.md | 1 + docs/source/development/release_process.md | 1 + docs/source/index.md | 1 + docs/source/users/constants.md | 1 + docs/source/users/demography/community.md | 1 + docs/source/users/demography/flora.md | 1 + .../users/demography/module_overview.md | 1 + docs/source/users/demography/t_model.md | 1 + docs/source/users/hygro.md | 1 + docs/source/users/pmodel/c3c4model.md | 1 + .../users/pmodel/isotopic_discrimination.md | 1 + docs/source/users/pmodel/module_overview.md | 1 + .../pmodel_details/envt_variation_outputs.md | 1 + .../pmodel/pmodel_details/extreme_values.md | 6 +---- .../pmodel/pmodel_details/jmax_limitation.md | 1 + .../pmodel/pmodel_details/optimal_chi.md | 1 + .../photosynthetic_environment.md | 1 + .../pmodel/pmodel_details/pmodel_overview.md | 1 + .../pmodel/pmodel_details/quantum_yield.md | 1 + .../users/pmodel/pmodel_details/rpmodel.md | 1 + .../pmodel/pmodel_details/soil_moisture.md | 1 + .../pmodel/pmodel_details/worked_examples.md | 25 +------------------ .../pmodel/subdaily_details/acclimation.md | 1 + .../subdaily_details/subdaily_calculations.md | 24 +----------------- .../subdaily_model_and_missing_data.md | 1 + .../subdaily_details/subdaily_overview.md | 1 + .../pmodel/subdaily_details/worked_example.md | 21 +--------------- docs/source/users/splash.md | 1 + docs/source/users/tmodel/canopy.md | 1 + docs/source/users/tmodel/tmodel.md | 1 + 42 files changed, 42 insertions(+), 72 deletions(-) diff --git a/docs/source/api/constants_api.md b/docs/source/api/constants_api.md index 1b992418..80bd4ba2 100644 --- a/docs/source/api/constants_api.md +++ b/docs/source/api/constants_api.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/api/core_api.md b/docs/source/api/core_api.md index 1324720e..89843589 100644 --- a/docs/source/api/core_api.md +++ b/docs/source/api/core_api.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/api/demography_api.md b/docs/source/api/demography_api.md index 9290a0b0..6700bad9 100644 --- a/docs/source/api/demography_api.md +++ b/docs/source/api/demography_api.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/api/pmodel_api.md b/docs/source/api/pmodel_api.md index 9bb4f669..0d385f1f 100644 --- a/docs/source/api/pmodel_api.md +++ b/docs/source/api/pmodel_api.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/api/splash_api.md b/docs/source/api/splash_api.md index fe55c923..e8177ab3 100644 --- a/docs/source/api/splash_api.md +++ b/docs/source/api/splash_api.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/api/tmodel_api.md b/docs/source/api/tmodel_api.md index b7fa9312..d74b7989 100644 --- a/docs/source/api/tmodel_api.md +++ b/docs/source/api/tmodel_api.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/development/code_qa_and_typing.md b/docs/source/development/code_qa_and_typing.md index 848f7a5c..b951cd8a 100644 --- a/docs/source/development/code_qa_and_typing.md +++ b/docs/source/development/code_qa_and_typing.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/development/code_testing.md b/docs/source/development/code_testing.md index 43829394..9d126f5a 100644 --- a/docs/source/development/code_testing.md +++ b/docs/source/development/code_testing.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/development/documentation.md b/docs/source/development/documentation.md index 0ce6b2f8..cbe8cc3f 100644 --- a/docs/source/development/documentation.md +++ b/docs/source/development/documentation.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/development/github_actions.md b/docs/source/development/github_actions.md index d26a7452..fda8ad3d 100644 --- a/docs/source/development/github_actions.md +++ b/docs/source/development/github_actions.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/development/overview.md b/docs/source/development/overview.md index 837b1efe..72207d53 100644 --- a/docs/source/development/overview.md +++ b/docs/source/development/overview.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/development/profiling_and_benchmarking.md b/docs/source/development/profiling_and_benchmarking.md index d705d43f..ccafd4f3 100644 --- a/docs/source/development/profiling_and_benchmarking.md +++ b/docs/source/development/profiling_and_benchmarking.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/development/pyrealm_build_data.md b/docs/source/development/pyrealm_build_data.md index b59ec2b2..0f9dd7c1 100644 --- a/docs/source/development/pyrealm_build_data.md +++ b/docs/source/development/pyrealm_build_data.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/development/release_process.md b/docs/source/development/release_process.md index e9ff6b00..9396484c 100644 --- a/docs/source/development/release_process.md +++ b/docs/source/development/release_process.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/index.md b/docs/source/index.md index 43aef59b..10234635 100644 --- a/docs/source/index.md +++ b/docs/source/index.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/constants.md b/docs/source/users/constants.md index 899076b3..7997e3a9 100644 --- a/docs/source/users/constants.md +++ b/docs/source/users/constants.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/demography/community.md b/docs/source/users/demography/community.md index 0961b8e2..15dd5954 100644 --- a/docs/source/users/demography/community.md +++ b/docs/source/users/demography/community.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python diff --git a/docs/source/users/demography/flora.md b/docs/source/users/demography/flora.md index 4b9bf3f7..541ad5e2 100644 --- a/docs/source/users/demography/flora.md +++ b/docs/source/users/demography/flora.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python diff --git a/docs/source/users/demography/module_overview.md b/docs/source/users/demography/module_overview.md index e2aaceab..f4d2e47c 100644 --- a/docs/source/users/demography/module_overview.md +++ b/docs/source/users/demography/module_overview.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python diff --git a/docs/source/users/demography/t_model.md b/docs/source/users/demography/t_model.md index 23510daa..f23fcbc7 100644 --- a/docs/source/users/demography/t_model.md +++ b/docs/source/users/demography/t_model.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python diff --git a/docs/source/users/hygro.md b/docs/source/users/hygro.md index a235875e..39322beb 100644 --- a/docs/source/users/hygro.md +++ b/docs/source/users/hygro.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/c3c4model.md b/docs/source/users/pmodel/c3c4model.md index ac4ba204..4a400e50 100644 --- a/docs/source/users/pmodel/c3c4model.md +++ b/docs/source/users/pmodel/c3c4model.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/isotopic_discrimination.md b/docs/source/users/pmodel/isotopic_discrimination.md index 574b2f67..5020de17 100644 --- a/docs/source/users/pmodel/isotopic_discrimination.md +++ b/docs/source/users/pmodel/isotopic_discrimination.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/module_overview.md b/docs/source/users/pmodel/module_overview.md index 04c29687..26d69734 100644 --- a/docs/source/users/pmodel/module_overview.md +++ b/docs/source/users/pmodel/module_overview.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/pmodel_details/envt_variation_outputs.md b/docs/source/users/pmodel/pmodel_details/envt_variation_outputs.md index e7364837..eae62e95 100644 --- a/docs/source/users/pmodel/pmodel_details/envt_variation_outputs.md +++ b/docs/source/users/pmodel/pmodel_details/envt_variation_outputs.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/pmodel_details/extreme_values.md b/docs/source/users/pmodel/pmodel_details/extreme_values.md index 7119ac7f..b9d28fcc 100644 --- a/docs/source/users/pmodel/pmodel_details/extreme_values.md +++ b/docs/source/users/pmodel/pmodel_details/extreme_values.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python @@ -53,7 +54,6 @@ Note that the default values for C3 photosynthesis give **non-zero values below ```{code-cell} :tags: [hide-input] -:trusted: true from matplotlib import pyplot import numpy as np @@ -93,7 +93,6 @@ that again, $\Gamma^_$ has non-zero values for sub-zero temperatures. ```{code-cell} :tags: [hide-input] -:trusted: true # Calculate gammastar at different pressures tc_1d = np.linspace(-80, 100, n_pts) @@ -119,7 +118,6 @@ temperature and atmospheric pressure and again behaves smoothly with extreme val ```{code-cell} :tags: [hide-input] -:trusted: true fig, ax = pyplot.subplots(1, 1) @@ -144,7 +142,6 @@ issue with low temperatures arising from the equations for the density of water. ```{code-cell} :tags: [hide-input] -:trusted: true fig, ax = pyplot.subplots(1, 1) @@ -169,7 +166,6 @@ predictions below about -30 °C. ```{code-cell} :tags: [hide-input] -:trusted: true fig, ax = pyplot.subplots(1, 1) diff --git a/docs/source/users/pmodel/pmodel_details/jmax_limitation.md b/docs/source/users/pmodel/pmodel_details/jmax_limitation.md index 5e07802f..84b22e97 100644 --- a/docs/source/users/pmodel/pmodel_details/jmax_limitation.md +++ b/docs/source/users/pmodel/pmodel_details/jmax_limitation.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python diff --git a/docs/source/users/pmodel/pmodel_details/optimal_chi.md b/docs/source/users/pmodel/pmodel_details/optimal_chi.md index fb22b49e..4c794fe0 100644 --- a/docs/source/users/pmodel/pmodel_details/optimal_chi.md +++ b/docs/source/users/pmodel/pmodel_details/optimal_chi.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/pmodel_details/photosynthetic_environment.md b/docs/source/users/pmodel/pmodel_details/photosynthetic_environment.md index cedd27fe..adbbca1a 100644 --- a/docs/source/users/pmodel/pmodel_details/photosynthetic_environment.md +++ b/docs/source/users/pmodel/pmodel_details/photosynthetic_environment.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/pmodel_details/pmodel_overview.md b/docs/source/users/pmodel/pmodel_details/pmodel_overview.md index 9626ee98..1343a4b5 100644 --- a/docs/source/users/pmodel/pmodel_details/pmodel_overview.md +++ b/docs/source/users/pmodel/pmodel_details/pmodel_overview.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python diff --git a/docs/source/users/pmodel/pmodel_details/quantum_yield.md b/docs/source/users/pmodel/pmodel_details/quantum_yield.md index 92a3ae8c..dd6e3600 100644 --- a/docs/source/users/pmodel/pmodel_details/quantum_yield.md +++ b/docs/source/users/pmodel/pmodel_details/quantum_yield.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python diff --git a/docs/source/users/pmodel/pmodel_details/rpmodel.md b/docs/source/users/pmodel/pmodel_details/rpmodel.md index 7e691cc2..f2b4eb93 100644 --- a/docs/source/users/pmodel/pmodel_details/rpmodel.md +++ b/docs/source/users/pmodel/pmodel_details/rpmodel.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/pmodel_details/soil_moisture.md b/docs/source/users/pmodel/pmodel_details/soil_moisture.md index 2b4b8a53..afe55613 100644 --- a/docs/source/users/pmodel/pmodel_details/soil_moisture.md +++ b/docs/source/users/pmodel/pmodel_details/soil_moisture.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/pmodel_details/worked_examples.md b/docs/source/users/pmodel/pmodel_details/worked_examples.md index f3db6947..ccd1d0c8 100644 --- a/docs/source/users/pmodel/pmodel_details/worked_examples.md +++ b/docs/source/users/pmodel/pmodel_details/worked_examples.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python @@ -39,8 +40,6 @@ The example shows the steps required using a single site with: ### Estimate photosynthetic environment ```{code-cell} -:trusted: true - from importlib import resources from matplotlib import pyplot as plt @@ -61,14 +60,10 @@ terse - just the shape of the data - but the more detailed summary of the attributes. ```{code-cell} -:trusted: true - env ``` ```{code-cell} -:trusted: true - env.summarize() ``` @@ -78,8 +73,6 @@ Next, the P Model can be fitted to the photosynthetic environment using the ({class}`~pyrealm.pmodel.pmodel.PModel`) class: ```{code-cell} -:trusted: true - model = PModel(env) ``` @@ -87,8 +80,6 @@ The returned model object holds a lot of information. The representation of the model object shows a terse display of the settings used to run the model: ```{code-cell} -:trusted: true - model ``` @@ -99,8 +90,6 @@ photosynthetic efficiency: the intrinsic water use efficiency (``iwue``) and the use efficiency (``lue``). ```{code-cell} -:trusted: true - model.summarize() ``` @@ -113,8 +102,6 @@ This object also has a {meth}`~pyrealm.pmodel.optimal_chi.OptimalChiABC.summariz method: ```{code-cell} -:trusted: true - model.optchi.summarize() ``` @@ -131,8 +118,6 @@ Here we are using: * a PPFD of 834 µmol m-2 s-1. ```{code-cell} -:trusted: true - model.estimate_productivity(fapar=0.91, ppfd=834) model.summarize() ``` @@ -165,8 +150,6 @@ to be the same size so some of the variables have repeated data across dimension * Elevation is constant across months, so the data for each month is repeated. ```{code-cell} -:trusted: true - # Load an example dataset containing the forcing variables. data_path = resources.files("pyrealm_build_data.rpmodel") / "pmodel_global.nc" ds = xarray.load_dataset(data_path) @@ -186,8 +169,6 @@ data to atmospheric pressure, and then this is used to set the photosynthetic environment for the model: ```{code-cell} -:trusted: true - # Convert elevation to atmospheric pressure patm = calc_patm(elev) @@ -206,8 +187,6 @@ That environment can then be run to calculate the P model predictions for light efficiency: ```{code-cell} -:trusted: true - # Run the P model model = PModel(env) @@ -221,8 +200,6 @@ Finally, the light use efficiency can be used to calculate GPP given the photosynthetic photon flux density and fAPAR. ```{code-cell} -:trusted: true - # Scale the outputs from values per unit iabs to realised values model.estimate_productivity(fapar, ppfd) diff --git a/docs/source/users/pmodel/subdaily_details/acclimation.md b/docs/source/users/pmodel/subdaily_details/acclimation.md index 193b38f9..09512e49 100644 --- a/docs/source/users/pmodel/subdaily_details/acclimation.md +++ b/docs/source/users/pmodel/subdaily_details/acclimation.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/subdaily_details/subdaily_calculations.md b/docs/source/users/pmodel/subdaily_details/subdaily_calculations.md index 1a134c84..f16c6831 100644 --- a/docs/source/users/pmodel/subdaily_details/subdaily_calculations.md +++ b/docs/source/users/pmodel/subdaily_details/subdaily_calculations.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python @@ -21,7 +22,6 @@ are handled internally by the model fitting in `pyrealm`. ```{code-cell} :tags: [hide-input] -:trusted: true from importlib import resources @@ -48,8 +48,6 @@ site](https://fluxnet.org/doi/FLUXNET2015/BE-Vie), which was also used as a demonstration in {cite:t}`mengoli:2022a`. ```{code-cell} -:trusted: true - data_path = resources.files("pyrealm_build_data.subdaily") / "subdaily_BE_Vie_2014.csv" data = pandas.read_csv(str(data_path)) @@ -71,8 +69,6 @@ subdaily timescale. The code below also estimates GPP under the standard P Model slow responses for comparison. ```{code-cell} -:trusted: true - # Calculate the photosynthetic environment subdaily_env = PModelEnvironment( tc=temp_subdaily, @@ -99,8 +95,6 @@ conditions at the observation closest to noon, or the mean environmental conditi window around noon. ```{code-cell} -:trusted: true - # Create the fast slow scaler fsscaler = SubdailyScaler(datetime_subdaily) @@ -122,7 +116,6 @@ pmodel_subdaily = SubdailyPModel( ```{code-cell} :tags: [hide-input] -:trusted: true idx = np.arange(48 * 120, 48 * 130) plt.figure(figsize=(10, 4)) @@ -146,8 +139,6 @@ inputs to the standard P Model to calculate the optimal behaviour of plants unde conditions. ```{code-cell} -:trusted: true - # Get the daily acclimation conditions for the forcing variables temp_acclim = fsscaler.get_daily_means(temp_subdaily) co2_acclim = fsscaler.get_daily_means(co2_subdaily) @@ -179,8 +170,6 @@ at 25°C. This is acheived by multiplying by the reciprocal of the exponential p the Arrhenius equation ($h^{-1}$ in {cite}`mengoli:2022a`). ```{code-cell} -:trusted: true - # Are these any of the existing values in the constants? ha_vcmax25 = 65330 ha_jmax25 = 43900 @@ -196,8 +185,6 @@ The memory effect can now be applied to the three parameters with slow responses to calculate realised values, here using the default 15 day window. ```{code-cell} -:trusted: true - # Calculation of memory effect in xi, vcmax25 and jmax25 xi_real = memory_effect(pmodel_acclim.optchi.xi, alpha=1 / 15) vcmax25_real = memory_effect(vcmax25_acclim, alpha=1 / 15, allow_holdover=True) @@ -210,7 +197,6 @@ application of the memory effect. ```{code-cell} :tags: [hide-input] -:trusted: true fig, axes = plt.subplots(1, 3, figsize=(16, 5)) @@ -244,8 +230,6 @@ temperature at fast scales: responses of $J_{max}$ and $V_{cmax}$. ```{code-cell} -:trusted: true - tk_subdaily = subdaily_env.tc + pmodel_subdaily.env.core_const.k_CtoK # Fill the realised jmax and vcmax from subdaily to daily @@ -266,8 +250,6 @@ optimal $\chi$, rather than calculating the instantaneously optimal values of $\ as is the case in the standard P Model. ```{code-cell} -:trusted: true - # Interpolate xi to subdaily scale xi_subdaily = fsscaler.fill_daily_to_subdaily(xi_real) @@ -287,8 +269,6 @@ include the slow responses of $V_{cmax25}$ and $J_{max25}$ and fast responses to temperature. ```{code-cell} -:trusted: true - # Calculate Ac Ac_subdaily = ( vcmax_subdaily @@ -319,7 +299,5 @@ print(np.nanmin(diff), np.nanmax(diff)) ``` ```{code-cell} -:trusted: true - ``` diff --git a/docs/source/users/pmodel/subdaily_details/subdaily_model_and_missing_data.md b/docs/source/users/pmodel/subdaily_details/subdaily_model_and_missing_data.md index 881f71b8..5f1e995b 100644 --- a/docs/source/users/pmodel/subdaily_details/subdaily_model_and_missing_data.md +++ b/docs/source/users/pmodel/subdaily_details/subdaily_model_and_missing_data.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/subdaily_details/subdaily_overview.md b/docs/source/users/pmodel/subdaily_details/subdaily_overview.md index 3278b056..31234000 100644 --- a/docs/source/users/pmodel/subdaily_details/subdaily_overview.md +++ b/docs/source/users/pmodel/subdaily_details/subdaily_overview.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/subdaily_details/worked_example.md b/docs/source/users/pmodel/subdaily_details/worked_example.md index fac12b6a..4d60302d 100644 --- a/docs/source/users/pmodel/subdaily_details/worked_example.md +++ b/docs/source/users/pmodel/subdaily_details/worked_example.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python @@ -14,8 +15,6 @@ kernelspec: # Worked example of the Subdaily P Model ```{code-cell} -:trusted: true - from importlib import resources import xarray @@ -61,8 +60,6 @@ The test data use some UK WFDE data for three sites in order to compare predicti over a time series. ```{code-cell} -:trusted: true - # Loading the example dataset: dpath = ( resources.files("pyrealm_build_data.uk_data") / "UK_WFDE5_FAPAR_2018_JuneJuly.nc" @@ -84,8 +81,6 @@ The WFDE data need some conversion for use in the PModel, along with the definit the atmospheric CO2 concentration. ```{code-cell} -:trusted: true - # Variable set up # Air temperature in °C from Tair in Kelvin tc = (ds["Tair"] - 273.15).to_numpy() @@ -105,8 +100,6 @@ co2 = np.ones_like(tc) * 400 The code below then calculates the photosynthetic environment. ```{code-cell} -:trusted: true - # Generate and check the PModelEnvironment pm_env = PModelEnvironment(tc=tc, patm=patm, vpd=vpd, co2=co2) pm_env.summarize() @@ -118,8 +111,6 @@ The standard implementation of the P Model used below assumes that plants can instantaneously adopt optimal behaviour. ```{code-cell} -:trusted: true - # Standard PModels pmodC3 = PModel( env=pm_env, method_kphio="fixed", reference_kphio=1 / 8, method_optchi="prentice14" @@ -129,8 +120,6 @@ pmodC3.summarize() ``` ```{code-cell} -:trusted: true - pmodC4 = PModel( env=pm_env, method_kphio="fixed", reference_kphio=1 / 8, method_optchi="c4_no_gamma" ) @@ -147,8 +136,6 @@ calculations: essentially the plant does not acclimate until the optimal values calculated again to update those realised estimates. ```{code-cell} -:trusted: true - # Set the acclimation window to an hour either side of noon fsscaler = SubdailyScaler(datetimes) fsscaler.set_window( @@ -188,8 +175,6 @@ shown above and plots the instantaneous predictions against predictions includin photosynthetic responses. ```{code-cell} -:trusted: true - # Store the predictions in the xarray Dataset to use indexing ds["GPP_pmodC3"] = (ds["Tair"].dims, pmodC3.gpp) ds["GPP_subdailyC3"] = (ds["Tair"].dims, subdailyC3.gpp) @@ -247,8 +232,6 @@ The subdaily models can also be obtained directly from the standard models, usin `convert_pmodel_to_subdaily` method: ```{code-cell} -:trusted: true - # Convert standard C3 model converted_C3 = convert_pmodel_to_subdaily( pmodel=pmodC3, @@ -270,8 +253,6 @@ This produces the same outputs as the `SubdailyPModel` class, but is convenient compact when the two models are going to be compared. ```{code-cell} -:trusted: true - # Models have identical GPP - maximum absolute difference is zero. print(np.nanmax(abs(subdailyC3.gpp.flatten() - converted_C3.gpp.flatten()))) print(np.nanmax(abs(subdailyC4.gpp.flatten() - converted_C4.gpp.flatten()))) diff --git a/docs/source/users/splash.md b/docs/source/users/splash.md index e5c4d283..e852f49c 100644 --- a/docs/source/users/splash.md +++ b/docs/source/users/splash.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/tmodel/canopy.md b/docs/source/users/tmodel/canopy.md index 3fab2df8..b5e49112 100644 --- a/docs/source/users/tmodel/canopy.md +++ b/docs/source/users/tmodel/canopy.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python diff --git a/docs/source/users/tmodel/tmodel.md b/docs/source/users/tmodel/tmodel.md index fda39de0..ca6c5eec 100644 --- a/docs/source/users/tmodel/tmodel.md +++ b/docs/source/users/tmodel/tmodel.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python From 0339465c5976943a4bb7cbfee784c44043c53a22 Mon Sep 17 00:00:00 2001 From: David Orme Date: Fri, 11 Oct 2024 15:16:13 +0100 Subject: [PATCH 18/24] Fix mypy issue --- docs/source/api/constants_api.md | 1 + docs/source/api/core_api.md | 1 + docs/source/api/demography_api.md | 1 + docs/source/api/pmodel_api.md | 1 + docs/source/api/splash_api.md | 1 + docs/source/api/tmodel_api.md | 1 + docs/source/development/code_qa_and_typing.md | 1 + docs/source/development/code_testing.md | 1 + docs/source/development/documentation.md | 1 + docs/source/development/github_actions.md | 1 + docs/source/development/overview.md | 1 + .../development/profiling_and_benchmarking.md | 1 + docs/source/development/pyrealm_build_data.md | 1 + docs/source/development/release_process.md | 1 + docs/source/index.md | 1 + docs/source/users/constants.md | 1 + docs/source/users/demography/community.md | 1 + docs/source/users/demography/flora.md | 1 + .../users/demography/module_overview.md | 1 + docs/source/users/demography/t_model.md | 1 + docs/source/users/hygro.md | 1 + docs/source/users/pmodel/c3c4model.md | 1 + .../users/pmodel/isotopic_discrimination.md | 1 + docs/source/users/pmodel/module_overview.md | 1 + .../pmodel_details/envt_variation_outputs.md | 1 + .../pmodel/pmodel_details/extreme_values.md | 6 +---- .../pmodel/pmodel_details/jmax_limitation.md | 1 + .../pmodel/pmodel_details/optimal_chi.md | 1 + .../photosynthetic_environment.md | 1 + .../pmodel/pmodel_details/pmodel_overview.md | 1 + .../pmodel/pmodel_details/quantum_yield.md | 1 + .../users/pmodel/pmodel_details/rpmodel.md | 1 + .../pmodel/pmodel_details/soil_moisture.md | 1 + .../pmodel/pmodel_details/worked_examples.md | 25 +------------------ .../pmodel/subdaily_details/acclimation.md | 1 + .../subdaily_details/subdaily_calculations.md | 24 +----------------- .../subdaily_model_and_missing_data.md | 1 + .../subdaily_details/subdaily_overview.md | 1 + .../pmodel/subdaily_details/worked_example.md | 21 +--------------- docs/source/users/splash.md | 1 + docs/source/users/tmodel/canopy.md | 1 + docs/source/users/tmodel/tmodel.md | 1 + pyrealm/demography/crown.py | 2 +- 43 files changed, 43 insertions(+), 73 deletions(-) diff --git a/docs/source/api/constants_api.md b/docs/source/api/constants_api.md index 1b992418..80bd4ba2 100644 --- a/docs/source/api/constants_api.md +++ b/docs/source/api/constants_api.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/api/core_api.md b/docs/source/api/core_api.md index 1324720e..89843589 100644 --- a/docs/source/api/core_api.md +++ b/docs/source/api/core_api.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/api/demography_api.md b/docs/source/api/demography_api.md index 9290a0b0..6700bad9 100644 --- a/docs/source/api/demography_api.md +++ b/docs/source/api/demography_api.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/api/pmodel_api.md b/docs/source/api/pmodel_api.md index 9bb4f669..0d385f1f 100644 --- a/docs/source/api/pmodel_api.md +++ b/docs/source/api/pmodel_api.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/api/splash_api.md b/docs/source/api/splash_api.md index fe55c923..e8177ab3 100644 --- a/docs/source/api/splash_api.md +++ b/docs/source/api/splash_api.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/api/tmodel_api.md b/docs/source/api/tmodel_api.md index b7fa9312..d74b7989 100644 --- a/docs/source/api/tmodel_api.md +++ b/docs/source/api/tmodel_api.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/development/code_qa_and_typing.md b/docs/source/development/code_qa_and_typing.md index 848f7a5c..b951cd8a 100644 --- a/docs/source/development/code_qa_and_typing.md +++ b/docs/source/development/code_qa_and_typing.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/development/code_testing.md b/docs/source/development/code_testing.md index 43829394..9d126f5a 100644 --- a/docs/source/development/code_testing.md +++ b/docs/source/development/code_testing.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/development/documentation.md b/docs/source/development/documentation.md index 0ce6b2f8..cbe8cc3f 100644 --- a/docs/source/development/documentation.md +++ b/docs/source/development/documentation.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/development/github_actions.md b/docs/source/development/github_actions.md index d26a7452..fda8ad3d 100644 --- a/docs/source/development/github_actions.md +++ b/docs/source/development/github_actions.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/development/overview.md b/docs/source/development/overview.md index 837b1efe..72207d53 100644 --- a/docs/source/development/overview.md +++ b/docs/source/development/overview.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/development/profiling_and_benchmarking.md b/docs/source/development/profiling_and_benchmarking.md index d705d43f..ccafd4f3 100644 --- a/docs/source/development/profiling_and_benchmarking.md +++ b/docs/source/development/profiling_and_benchmarking.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/development/pyrealm_build_data.md b/docs/source/development/pyrealm_build_data.md index b59ec2b2..0f9dd7c1 100644 --- a/docs/source/development/pyrealm_build_data.md +++ b/docs/source/development/pyrealm_build_data.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/development/release_process.md b/docs/source/development/release_process.md index e9ff6b00..9396484c 100644 --- a/docs/source/development/release_process.md +++ b/docs/source/development/release_process.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/index.md b/docs/source/index.md index 43aef59b..10234635 100644 --- a/docs/source/index.md +++ b/docs/source/index.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/constants.md b/docs/source/users/constants.md index 899076b3..7997e3a9 100644 --- a/docs/source/users/constants.md +++ b/docs/source/users/constants.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/demography/community.md b/docs/source/users/demography/community.md index 0961b8e2..15dd5954 100644 --- a/docs/source/users/demography/community.md +++ b/docs/source/users/demography/community.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python diff --git a/docs/source/users/demography/flora.md b/docs/source/users/demography/flora.md index 4b9bf3f7..541ad5e2 100644 --- a/docs/source/users/demography/flora.md +++ b/docs/source/users/demography/flora.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python diff --git a/docs/source/users/demography/module_overview.md b/docs/source/users/demography/module_overview.md index e2aaceab..f4d2e47c 100644 --- a/docs/source/users/demography/module_overview.md +++ b/docs/source/users/demography/module_overview.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python diff --git a/docs/source/users/demography/t_model.md b/docs/source/users/demography/t_model.md index 23510daa..f23fcbc7 100644 --- a/docs/source/users/demography/t_model.md +++ b/docs/source/users/demography/t_model.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python diff --git a/docs/source/users/hygro.md b/docs/source/users/hygro.md index a235875e..39322beb 100644 --- a/docs/source/users/hygro.md +++ b/docs/source/users/hygro.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/c3c4model.md b/docs/source/users/pmodel/c3c4model.md index ac4ba204..4a400e50 100644 --- a/docs/source/users/pmodel/c3c4model.md +++ b/docs/source/users/pmodel/c3c4model.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/isotopic_discrimination.md b/docs/source/users/pmodel/isotopic_discrimination.md index 574b2f67..5020de17 100644 --- a/docs/source/users/pmodel/isotopic_discrimination.md +++ b/docs/source/users/pmodel/isotopic_discrimination.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/module_overview.md b/docs/source/users/pmodel/module_overview.md index 04c29687..26d69734 100644 --- a/docs/source/users/pmodel/module_overview.md +++ b/docs/source/users/pmodel/module_overview.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/pmodel_details/envt_variation_outputs.md b/docs/source/users/pmodel/pmodel_details/envt_variation_outputs.md index e7364837..eae62e95 100644 --- a/docs/source/users/pmodel/pmodel_details/envt_variation_outputs.md +++ b/docs/source/users/pmodel/pmodel_details/envt_variation_outputs.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/pmodel_details/extreme_values.md b/docs/source/users/pmodel/pmodel_details/extreme_values.md index 7119ac7f..b9d28fcc 100644 --- a/docs/source/users/pmodel/pmodel_details/extreme_values.md +++ b/docs/source/users/pmodel/pmodel_details/extreme_values.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python @@ -53,7 +54,6 @@ Note that the default values for C3 photosynthesis give **non-zero values below ```{code-cell} :tags: [hide-input] -:trusted: true from matplotlib import pyplot import numpy as np @@ -93,7 +93,6 @@ that again, $\Gamma^_$ has non-zero values for sub-zero temperatures. ```{code-cell} :tags: [hide-input] -:trusted: true # Calculate gammastar at different pressures tc_1d = np.linspace(-80, 100, n_pts) @@ -119,7 +118,6 @@ temperature and atmospheric pressure and again behaves smoothly with extreme val ```{code-cell} :tags: [hide-input] -:trusted: true fig, ax = pyplot.subplots(1, 1) @@ -144,7 +142,6 @@ issue with low temperatures arising from the equations for the density of water. ```{code-cell} :tags: [hide-input] -:trusted: true fig, ax = pyplot.subplots(1, 1) @@ -169,7 +166,6 @@ predictions below about -30 °C. ```{code-cell} :tags: [hide-input] -:trusted: true fig, ax = pyplot.subplots(1, 1) diff --git a/docs/source/users/pmodel/pmodel_details/jmax_limitation.md b/docs/source/users/pmodel/pmodel_details/jmax_limitation.md index 5e07802f..84b22e97 100644 --- a/docs/source/users/pmodel/pmodel_details/jmax_limitation.md +++ b/docs/source/users/pmodel/pmodel_details/jmax_limitation.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python diff --git a/docs/source/users/pmodel/pmodel_details/optimal_chi.md b/docs/source/users/pmodel/pmodel_details/optimal_chi.md index fb22b49e..4c794fe0 100644 --- a/docs/source/users/pmodel/pmodel_details/optimal_chi.md +++ b/docs/source/users/pmodel/pmodel_details/optimal_chi.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/pmodel_details/photosynthetic_environment.md b/docs/source/users/pmodel/pmodel_details/photosynthetic_environment.md index cedd27fe..adbbca1a 100644 --- a/docs/source/users/pmodel/pmodel_details/photosynthetic_environment.md +++ b/docs/source/users/pmodel/pmodel_details/photosynthetic_environment.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/pmodel_details/pmodel_overview.md b/docs/source/users/pmodel/pmodel_details/pmodel_overview.md index 9626ee98..1343a4b5 100644 --- a/docs/source/users/pmodel/pmodel_details/pmodel_overview.md +++ b/docs/source/users/pmodel/pmodel_details/pmodel_overview.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python diff --git a/docs/source/users/pmodel/pmodel_details/quantum_yield.md b/docs/source/users/pmodel/pmodel_details/quantum_yield.md index 92a3ae8c..dd6e3600 100644 --- a/docs/source/users/pmodel/pmodel_details/quantum_yield.md +++ b/docs/source/users/pmodel/pmodel_details/quantum_yield.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python diff --git a/docs/source/users/pmodel/pmodel_details/rpmodel.md b/docs/source/users/pmodel/pmodel_details/rpmodel.md index 7e691cc2..f2b4eb93 100644 --- a/docs/source/users/pmodel/pmodel_details/rpmodel.md +++ b/docs/source/users/pmodel/pmodel_details/rpmodel.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/pmodel_details/soil_moisture.md b/docs/source/users/pmodel/pmodel_details/soil_moisture.md index 2b4b8a53..afe55613 100644 --- a/docs/source/users/pmodel/pmodel_details/soil_moisture.md +++ b/docs/source/users/pmodel/pmodel_details/soil_moisture.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/pmodel_details/worked_examples.md b/docs/source/users/pmodel/pmodel_details/worked_examples.md index f3db6947..ccd1d0c8 100644 --- a/docs/source/users/pmodel/pmodel_details/worked_examples.md +++ b/docs/source/users/pmodel/pmodel_details/worked_examples.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python @@ -39,8 +40,6 @@ The example shows the steps required using a single site with: ### Estimate photosynthetic environment ```{code-cell} -:trusted: true - from importlib import resources from matplotlib import pyplot as plt @@ -61,14 +60,10 @@ terse - just the shape of the data - but the more detailed summary of the attributes. ```{code-cell} -:trusted: true - env ``` ```{code-cell} -:trusted: true - env.summarize() ``` @@ -78,8 +73,6 @@ Next, the P Model can be fitted to the photosynthetic environment using the ({class}`~pyrealm.pmodel.pmodel.PModel`) class: ```{code-cell} -:trusted: true - model = PModel(env) ``` @@ -87,8 +80,6 @@ The returned model object holds a lot of information. The representation of the model object shows a terse display of the settings used to run the model: ```{code-cell} -:trusted: true - model ``` @@ -99,8 +90,6 @@ photosynthetic efficiency: the intrinsic water use efficiency (``iwue``) and the use efficiency (``lue``). ```{code-cell} -:trusted: true - model.summarize() ``` @@ -113,8 +102,6 @@ This object also has a {meth}`~pyrealm.pmodel.optimal_chi.OptimalChiABC.summariz method: ```{code-cell} -:trusted: true - model.optchi.summarize() ``` @@ -131,8 +118,6 @@ Here we are using: * a PPFD of 834 µmol m-2 s-1. ```{code-cell} -:trusted: true - model.estimate_productivity(fapar=0.91, ppfd=834) model.summarize() ``` @@ -165,8 +150,6 @@ to be the same size so some of the variables have repeated data across dimension * Elevation is constant across months, so the data for each month is repeated. ```{code-cell} -:trusted: true - # Load an example dataset containing the forcing variables. data_path = resources.files("pyrealm_build_data.rpmodel") / "pmodel_global.nc" ds = xarray.load_dataset(data_path) @@ -186,8 +169,6 @@ data to atmospheric pressure, and then this is used to set the photosynthetic environment for the model: ```{code-cell} -:trusted: true - # Convert elevation to atmospheric pressure patm = calc_patm(elev) @@ -206,8 +187,6 @@ That environment can then be run to calculate the P model predictions for light efficiency: ```{code-cell} -:trusted: true - # Run the P model model = PModel(env) @@ -221,8 +200,6 @@ Finally, the light use efficiency can be used to calculate GPP given the photosynthetic photon flux density and fAPAR. ```{code-cell} -:trusted: true - # Scale the outputs from values per unit iabs to realised values model.estimate_productivity(fapar, ppfd) diff --git a/docs/source/users/pmodel/subdaily_details/acclimation.md b/docs/source/users/pmodel/subdaily_details/acclimation.md index 193b38f9..09512e49 100644 --- a/docs/source/users/pmodel/subdaily_details/acclimation.md +++ b/docs/source/users/pmodel/subdaily_details/acclimation.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/subdaily_details/subdaily_calculations.md b/docs/source/users/pmodel/subdaily_details/subdaily_calculations.md index 1a134c84..f16c6831 100644 --- a/docs/source/users/pmodel/subdaily_details/subdaily_calculations.md +++ b/docs/source/users/pmodel/subdaily_details/subdaily_calculations.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python @@ -21,7 +22,6 @@ are handled internally by the model fitting in `pyrealm`. ```{code-cell} :tags: [hide-input] -:trusted: true from importlib import resources @@ -48,8 +48,6 @@ site](https://fluxnet.org/doi/FLUXNET2015/BE-Vie), which was also used as a demonstration in {cite:t}`mengoli:2022a`. ```{code-cell} -:trusted: true - data_path = resources.files("pyrealm_build_data.subdaily") / "subdaily_BE_Vie_2014.csv" data = pandas.read_csv(str(data_path)) @@ -71,8 +69,6 @@ subdaily timescale. The code below also estimates GPP under the standard P Model slow responses for comparison. ```{code-cell} -:trusted: true - # Calculate the photosynthetic environment subdaily_env = PModelEnvironment( tc=temp_subdaily, @@ -99,8 +95,6 @@ conditions at the observation closest to noon, or the mean environmental conditi window around noon. ```{code-cell} -:trusted: true - # Create the fast slow scaler fsscaler = SubdailyScaler(datetime_subdaily) @@ -122,7 +116,6 @@ pmodel_subdaily = SubdailyPModel( ```{code-cell} :tags: [hide-input] -:trusted: true idx = np.arange(48 * 120, 48 * 130) plt.figure(figsize=(10, 4)) @@ -146,8 +139,6 @@ inputs to the standard P Model to calculate the optimal behaviour of plants unde conditions. ```{code-cell} -:trusted: true - # Get the daily acclimation conditions for the forcing variables temp_acclim = fsscaler.get_daily_means(temp_subdaily) co2_acclim = fsscaler.get_daily_means(co2_subdaily) @@ -179,8 +170,6 @@ at 25°C. This is acheived by multiplying by the reciprocal of the exponential p the Arrhenius equation ($h^{-1}$ in {cite}`mengoli:2022a`). ```{code-cell} -:trusted: true - # Are these any of the existing values in the constants? ha_vcmax25 = 65330 ha_jmax25 = 43900 @@ -196,8 +185,6 @@ The memory effect can now be applied to the three parameters with slow responses to calculate realised values, here using the default 15 day window. ```{code-cell} -:trusted: true - # Calculation of memory effect in xi, vcmax25 and jmax25 xi_real = memory_effect(pmodel_acclim.optchi.xi, alpha=1 / 15) vcmax25_real = memory_effect(vcmax25_acclim, alpha=1 / 15, allow_holdover=True) @@ -210,7 +197,6 @@ application of the memory effect. ```{code-cell} :tags: [hide-input] -:trusted: true fig, axes = plt.subplots(1, 3, figsize=(16, 5)) @@ -244,8 +230,6 @@ temperature at fast scales: responses of $J_{max}$ and $V_{cmax}$. ```{code-cell} -:trusted: true - tk_subdaily = subdaily_env.tc + pmodel_subdaily.env.core_const.k_CtoK # Fill the realised jmax and vcmax from subdaily to daily @@ -266,8 +250,6 @@ optimal $\chi$, rather than calculating the instantaneously optimal values of $\ as is the case in the standard P Model. ```{code-cell} -:trusted: true - # Interpolate xi to subdaily scale xi_subdaily = fsscaler.fill_daily_to_subdaily(xi_real) @@ -287,8 +269,6 @@ include the slow responses of $V_{cmax25}$ and $J_{max25}$ and fast responses to temperature. ```{code-cell} -:trusted: true - # Calculate Ac Ac_subdaily = ( vcmax_subdaily @@ -319,7 +299,5 @@ print(np.nanmin(diff), np.nanmax(diff)) ``` ```{code-cell} -:trusted: true - ``` diff --git a/docs/source/users/pmodel/subdaily_details/subdaily_model_and_missing_data.md b/docs/source/users/pmodel/subdaily_details/subdaily_model_and_missing_data.md index 881f71b8..5f1e995b 100644 --- a/docs/source/users/pmodel/subdaily_details/subdaily_model_and_missing_data.md +++ b/docs/source/users/pmodel/subdaily_details/subdaily_model_and_missing_data.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/subdaily_details/subdaily_overview.md b/docs/source/users/pmodel/subdaily_details/subdaily_overview.md index 3278b056..31234000 100644 --- a/docs/source/users/pmodel/subdaily_details/subdaily_overview.md +++ b/docs/source/users/pmodel/subdaily_details/subdaily_overview.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/pmodel/subdaily_details/worked_example.md b/docs/source/users/pmodel/subdaily_details/worked_example.md index fac12b6a..4d60302d 100644 --- a/docs/source/users/pmodel/subdaily_details/worked_example.md +++ b/docs/source/users/pmodel/subdaily_details/worked_example.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python @@ -14,8 +15,6 @@ kernelspec: # Worked example of the Subdaily P Model ```{code-cell} -:trusted: true - from importlib import resources import xarray @@ -61,8 +60,6 @@ The test data use some UK WFDE data for three sites in order to compare predicti over a time series. ```{code-cell} -:trusted: true - # Loading the example dataset: dpath = ( resources.files("pyrealm_build_data.uk_data") / "UK_WFDE5_FAPAR_2018_JuneJuly.nc" @@ -84,8 +81,6 @@ The WFDE data need some conversion for use in the PModel, along with the definit the atmospheric CO2 concentration. ```{code-cell} -:trusted: true - # Variable set up # Air temperature in °C from Tair in Kelvin tc = (ds["Tair"] - 273.15).to_numpy() @@ -105,8 +100,6 @@ co2 = np.ones_like(tc) * 400 The code below then calculates the photosynthetic environment. ```{code-cell} -:trusted: true - # Generate and check the PModelEnvironment pm_env = PModelEnvironment(tc=tc, patm=patm, vpd=vpd, co2=co2) pm_env.summarize() @@ -118,8 +111,6 @@ The standard implementation of the P Model used below assumes that plants can instantaneously adopt optimal behaviour. ```{code-cell} -:trusted: true - # Standard PModels pmodC3 = PModel( env=pm_env, method_kphio="fixed", reference_kphio=1 / 8, method_optchi="prentice14" @@ -129,8 +120,6 @@ pmodC3.summarize() ``` ```{code-cell} -:trusted: true - pmodC4 = PModel( env=pm_env, method_kphio="fixed", reference_kphio=1 / 8, method_optchi="c4_no_gamma" ) @@ -147,8 +136,6 @@ calculations: essentially the plant does not acclimate until the optimal values calculated again to update those realised estimates. ```{code-cell} -:trusted: true - # Set the acclimation window to an hour either side of noon fsscaler = SubdailyScaler(datetimes) fsscaler.set_window( @@ -188,8 +175,6 @@ shown above and plots the instantaneous predictions against predictions includin photosynthetic responses. ```{code-cell} -:trusted: true - # Store the predictions in the xarray Dataset to use indexing ds["GPP_pmodC3"] = (ds["Tair"].dims, pmodC3.gpp) ds["GPP_subdailyC3"] = (ds["Tair"].dims, subdailyC3.gpp) @@ -247,8 +232,6 @@ The subdaily models can also be obtained directly from the standard models, usin `convert_pmodel_to_subdaily` method: ```{code-cell} -:trusted: true - # Convert standard C3 model converted_C3 = convert_pmodel_to_subdaily( pmodel=pmodC3, @@ -270,8 +253,6 @@ This produces the same outputs as the `SubdailyPModel` class, but is convenient compact when the two models are going to be compared. ```{code-cell} -:trusted: true - # Models have identical GPP - maximum absolute difference is zero. print(np.nanmax(abs(subdailyC3.gpp.flatten() - converted_C3.gpp.flatten()))) print(np.nanmax(abs(subdailyC4.gpp.flatten() - converted_C4.gpp.flatten()))) diff --git a/docs/source/users/splash.md b/docs/source/users/splash.md index e5c4d283..e852f49c 100644 --- a/docs/source/users/splash.md +++ b/docs/source/users/splash.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/docs/source/users/tmodel/canopy.md b/docs/source/users/tmodel/canopy.md index 3fab2df8..b5e49112 100644 --- a/docs/source/users/tmodel/canopy.md +++ b/docs/source/users/tmodel/canopy.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python diff --git a/docs/source/users/tmodel/tmodel.md b/docs/source/users/tmodel/tmodel.md index fda39de0..ca6c5eec 100644 --- a/docs/source/users/tmodel/tmodel.md +++ b/docs/source/users/tmodel/tmodel.md @@ -5,6 +5,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 + jupytext_version: 1.16.4 kernelspec: display_name: Python 3 language: python diff --git a/pyrealm/demography/crown.py b/pyrealm/demography/crown.py index 00f0019b..0f1d43ed 100644 --- a/pyrealm/demography/crown.py +++ b/pyrealm/demography/crown.py @@ -464,7 +464,7 @@ def get_crown_xy( attr_values = np.flip(attr_values, axis=0) # Collect the per stem data - crown_plotting_data: list[tuple[NDArray, NDArray]] | list[NDArray] = [] + crown_plotting_data: list = [] for stem_index in np.arange(attr_values.shape[1]): # Find the heights and values that fall within the individual stem From 848a96ef5ce97d90e2dcc7ed7c79fee167280022 Mon Sep 17 00:00:00 2001 From: David Orme Date: Fri, 11 Oct 2024 15:37:28 +0100 Subject: [PATCH 19/24] Fixing tests --- tests/unit/demography/conftest.py | 19 +++++++++ tests/unit/demography/test_canopy.py | 44 ++++++++++++++++++++- tests/unit/demography/test_crown.py | 58 +--------------------------- 3 files changed, 62 insertions(+), 59 deletions(-) diff --git a/tests/unit/demography/conftest.py b/tests/unit/demography/conftest.py index 5605db74..735cb4c1 100644 --- a/tests/unit/demography/conftest.py +++ b/tests/unit/demography/conftest.py @@ -63,6 +63,25 @@ def rtmodel_flora(): ) +@pytest.fixture +def fixture_community(): + """A fixture providing a simple community.""" + from pyrealm.demography.community import Community + from pyrealm.demography.flora import Flora, PlantFunctionalType + + # A simple community containing one sample stem, with an initial crown gap fraction + # of zero. + flora = Flora([PlantFunctionalType(name="test", f_g=0.0)]) + return Community( + cell_id=1, + cell_area=100, + flora=flora, + cohort_n_individuals=np.repeat([1], 4), + cohort_pft_names=np.repeat(["test"], 4), + cohort_dbh_values=np.array([0.2, 0.4, 0.6, 0.8]), + ) + + @pytest.fixture def rtmodel_data(): """Loads some simple predictions from the R implementation for testing.""" diff --git a/tests/unit/demography/test_canopy.py b/tests/unit/demography/test_canopy.py index 615f963a..33d457a5 100644 --- a/tests/unit/demography/test_canopy.py +++ b/tests/unit/demography/test_canopy.py @@ -1,6 +1,7 @@ """Testing the Canopy object.""" import numpy as np +import pytest def test_Canopy__init__(): @@ -31,7 +32,7 @@ def test_Canopy__init__(): ) canopy_gap_fraction = 0.05 - canopy = Canopy(community, canopy_gap_fraction=canopy_gap_fraction) + canopy = Canopy(community, canopy_gap_fraction=canopy_gap_fraction, fit_ppa=True) # Simply check that the shape of the stem leaf area matrix is the right shape n_layers_from_crown_area = int( @@ -46,4 +47,43 @@ def test_Canopy__init__(): / community.cell_area ) ) - assert canopy.stem_leaf_area.shape == (n_layers_from_crown_area, canopy.n_cohorts) + assert canopy.layer_stem_leaf_area.shape == ( + n_layers_from_crown_area, + canopy.n_cohorts, + ) + + +def test_solve_canopy_area_filling_height(fixture_community): + """Test solve_community_projected_canopy_area. + + The logic of this test is that given the cumulative sum of the crown areas in the + fixture from tallest to shortest as the target, providing the z_max of each stem as + the height _should_ always return zero, as this is exactly the height at which that + cumulative area would close: crown 1 closes at z_max 1, crown 1 + 2 closes at z_max + 2 and so on. + """ + + from pyrealm.demography.canopy import ( + solve_canopy_area_filling_height, + ) + + for ( + this_height, + this_target, + ) in zip( + np.flip(fixture_community.stem_allometry.crown_z_max), + np.cumsum(np.flip(fixture_community.stem_allometry.crown_area)), + ): + solved = solve_canopy_area_filling_height( + z=this_height, + stem_height=fixture_community.stem_allometry.stem_height, + crown_area=fixture_community.stem_allometry.crown_area, + n_individuals=fixture_community.cohort_data["n_individuals"], + m=fixture_community.stem_traits.m, + n=fixture_community.stem_traits.n, + q_m=fixture_community.stem_traits.q_m, + z_max=fixture_community.stem_allometry.crown_z_max, + target_area=this_target, + ) + + assert solved == pytest.approx(0) diff --git a/tests/unit/demography/test_crown.py b/tests/unit/demography/test_crown.py index 221eac2c..a7cd10ec 100644 --- a/tests/unit/demography/test_crown.py +++ b/tests/unit/demography/test_crown.py @@ -1,4 +1,4 @@ -"""test the functions in canopy_functions.py.""" +"""Test the functions in crown.py.""" from collections import namedtuple from contextlib import nullcontext as does_not_raise @@ -6,26 +6,6 @@ import numpy as np import pytest - -@pytest.fixture -def fixture_community(): - """A fixture providing a simple community.""" - from pyrealm.demography.community import Community - from pyrealm.demography.flora import Flora, PlantFunctionalType - - # A simple community containing one sample stem, with an initial crown gap fraction - # of zero. - flora = Flora([PlantFunctionalType(name="test", f_g=0.0)]) - return Community( - cell_id=1, - cell_area=100, - flora=flora, - cohort_n_individuals=np.repeat([1], 4), - cohort_pft_names=np.repeat(["test"], 4), - cohort_dbh_values=np.array([0.2, 0.4, 0.6, 0.8]), - ) - - ZQZInput = namedtuple( "ZQZInput", ["z", "stem", "more_stem", "q_z", "outcome", "excep_msg", "output_shape"], @@ -466,42 +446,6 @@ def test_calculate_stem_projected_crown_area_at_z_values( ) -def test_solve_community_projected_canopy_area(fixture_community): - """Test solve_community_projected_canopy_area. - - The logic of this test is that given the cumulative sum of the crown areas in the - fixture from tallest to shortest as the target, providing the z_max of each stem as - the height _should_ always return zero, as this is exactly the height at which that - cumulative area would close: crown 1 closes at z_max 1, crown 1 + 2 closes at z_max - 2 and so on. - """ - - from pyrealm.demography.crown import ( - solve_community_projected_canopy_area, - ) - - for ( - this_height, - this_target, - ) in zip( - np.flip(fixture_community.stem_allometry.crown_z_max), - np.cumsum(np.flip(fixture_community.stem_allometry.crown_area)), - ): - solved = solve_community_projected_canopy_area( - z=this_height, - stem_height=fixture_community.stem_allometry.stem_height, - crown_area=fixture_community.stem_allometry.crown_area, - n_individuals=fixture_community.cohort_data["n_individuals"], - m=fixture_community.stem_traits.m, - n=fixture_community.stem_traits.n, - q_m=fixture_community.stem_traits.q_m, - z_max=fixture_community.stem_allometry.crown_z_max, - target_area=this_target, - ) - - assert solved == pytest.approx(0) - - @pytest.mark.parametrize( argnames="fixture_z_qz_stem_properties", argvalues=[ From 69819e637e2a1133b260ac80d8b328078e826ae3 Mon Sep 17 00:00:00 2001 From: David Orme Date: Fri, 11 Oct 2024 16:01:43 +0100 Subject: [PATCH 20/24] Force matplotlib to build font cache early in conf.py --- docs/source/conf.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/docs/source/conf.py b/docs/source/conf.py index 9c3a2970..2de60f86 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -12,6 +12,9 @@ from dataclasses import dataclass, field from datetime import datetime +# Import Matplotlib to avoid this message in notebooks: +# "Matplotlib is building the font cache; this may take a moment." +import matplotlib.pyplot # noqa: F401 import sphinxcontrib.bibtex.plugin from sphinxcontrib.bibtex.style.referencing import BracketStyle from sphinxcontrib.bibtex.style.referencing.author_year import AuthorYearReferenceStyle From 81fabc059f285c2abddeb4f2a2536110e5522901 Mon Sep 17 00:00:00 2001 From: David Orme Date: Mon, 14 Oct 2024 14:24:36 +0100 Subject: [PATCH 21/24] 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"] From 802dc05fb775925e23c1f5ed03537a342a429c42 Mon Sep 17 00:00:00 2001 From: David Orme Date: Mon, 14 Oct 2024 14:37:14 +0100 Subject: [PATCH 22/24] Fix test --- tests/unit/demography/test_canopy.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit/demography/test_canopy.py b/tests/unit/demography/test_canopy.py index 33d457a5..a078fbfe 100644 --- a/tests/unit/demography/test_canopy.py +++ b/tests/unit/demography/test_canopy.py @@ -47,7 +47,7 @@ def test_Canopy__init__(): / community.cell_area ) ) - assert canopy.layer_stem_leaf_area.shape == ( + assert canopy.stem_leaf_area.shape == ( n_layers_from_crown_area, canopy.n_cohorts, ) From 4504b3bd8ff61f985337c7d9c8699966d57ee2ed Mon Sep 17 00:00:00 2001 From: David Orme Date: Wed, 16 Oct 2024 11:27:15 +0100 Subject: [PATCH 23/24] Simple community drawing --- docs/source/users/demography/canopy.md | 37 ++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/docs/source/users/demography/canopy.md b/docs/source/users/demography/canopy.md index f875cfa7..234abb56 100644 --- a/docs/source/users/demography/canopy.md +++ b/docs/source/users/demography/canopy.md @@ -35,6 +35,7 @@ notes and initial demonstration code. ```{code-cell} ipython3 from matplotlib import pyplot as plt import matplotlib as mpl +from matplotlib.patches import Polygon import numpy as np import pandas as pd @@ -228,7 +229,7 @@ and then three cohorts: * 7 saplings of the short PFT * 3 larger stems of the short PFT -* 1 large stems of tall PFT +* 2 large stems of tall PFT ```{code-cell} ipython3 # Define PFTs @@ -264,6 +265,33 @@ hghts = np.linspace(community.stem_allometry.stem_height.max(), 0, num=101)[:, N canopy = Canopy(community=community, layer_heights=hghts) ``` +The plot below then shows a simplistic 2D representation of the community. + +```{code-cell} ipython3 +fig, ax = plt.subplots(ncols=1) + +# Extract the crown profiles as XY arrays for plotting +profiles = get_crown_xy( + crown_profile=canopy.crown_profile, + stem_allometry=community.stem_allometry, + attr="crown_radius", + as_xy=True, +) + +for idx, crown in enumerate(profiles): + + # Get spaced but slightly randomized stem locations + n_stems = community.cohort_data["n_individuals"][idx] + stem_locations = np.linspace(0, 10, num=n_stems) + np.random.normal(size=n_stems) + + # Plot the crown model for each stem + for stem_loc in stem_locations: + ax.add_patch(Polygon(crown + np.array([stem_loc, 0]), color="#00550055")) + +ax.autoscale_view() +ax.set_aspect(1) +``` + As before, we can verify that the cumulative light extinction at the bottom of the vertical profile is equal to the expected value across the whole community. @@ -384,7 +412,7 @@ l_m = \left\lceil \frac{\sum_1^{N_s}{A_c}}{ A(1 - f_G)}\right\rceil $$ ```{code-cell} ipython3 -canopy_ppa = Canopy(community=community, canopy_gap_fraction=0, fit_ppa=True) +canopy_ppa = Canopy(community=community, canopy_gap_fraction=2 / 32, fit_ppa=True) ``` The `canopy_ppa.heights` attribute now contains the heights at which the PPA @@ -497,6 +525,11 @@ ax2_rhs.set_yticks(canopy_ppa.heights.flatten()) _ = ax2_rhs.set_yticklabels(z_star_labels) ``` +```{code-cell} ipython3 +t = 0.6 +(1 - 0.2) - t +``` + ## Light allocation From 6270aee574ed14e1f6dba0f5351b1168bc38372f Mon Sep 17 00:00:00 2001 From: David Orme Date: Mon, 21 Oct 2024 14:40:39 +0100 Subject: [PATCH 24/24] Link to crown profile attributes --- pyrealm/demography/crown.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyrealm/demography/crown.py b/pyrealm/demography/crown.py index 0f1d43ed..90481507 100644 --- a/pyrealm/demography/crown.py +++ b/pyrealm/demography/crown.py @@ -439,7 +439,8 @@ def get_crown_xy( Args: crown_profile: A crown profile instance stem_allometry: The stem allometry instance used to create the crown profile - attr: The crown profile attribute to plot + attr: The crown profile attribute to plot (see + :class:`~pyrealm.demography.crown.CrownProfile`) stem_offsets: An optional array of offsets to add to the midline of stems. two_sided: Should the plotting data show a two sided canopy. as_xy: Should the plotting data be returned as a single XY array.