diff --git a/docs/source/_toc.yml b/docs/source/_toc.yml index b91919bc..423064a0 100644 --- a/docs/source/_toc.yml +++ b/docs/source/_toc.yml @@ -36,8 +36,6 @@ subtrees: - file: users/demography/crown.md - file: users/demography/community.md - file: users/demography/canopy.md - - file: users/tmodel/tmodel.md - - file: users/tmodel/canopy.md - file: users/splash.md - file: users/constants.md - file: users/hygro.md diff --git a/docs/source/api/demography_api.md b/docs/source/api/demography_api.md index 0630e5a9..a5b08c28 100644 --- a/docs/source/api/demography_api.md +++ b/docs/source/api/demography_api.md @@ -46,10 +46,10 @@ language_info: :members: ``` -## The {mod}`~pyrealm.demography.t_model_functions` module +## The {mod}`~pyrealm.demography.tmodel` module ```{eval-rst} -.. automodule:: pyrealm.demography.t_model_functions +.. automodule:: pyrealm.demography.tmodel :autosummary: :members: ``` diff --git a/docs/source/api/tmodel_api.md b/docs/source/api/tmodel_api.md deleted file mode 100644 index de3c0158..00000000 --- a/docs/source/api/tmodel_api.md +++ /dev/null @@ -1,33 +0,0 @@ ---- -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 - 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 {mod}`~pyrealm.tmodel` module - -```{eval-rst} -.. automodule:: pyrealm.tmodel - :autosummary: - :members: - :special-members: __init__ - -``` diff --git a/docs/source/index.md b/docs/source/index.md index 32ec1e69..f8138181 100644 --- a/docs/source/index.md +++ b/docs/source/index.md @@ -30,18 +30,27 @@ following modules: ## Module overview +The `core` module +: Contains fundamental utilities and physics functionality shared across the + package, including the [hygro](users/hygro) and the utilities submodules. + The `pmodel` module : Fitting the [P Model](users/pmodel/module_overview), which is an ecophysiological model of optimal carbon dioxide uptake by plants {cite:p}`Prentice:2014bc, Wang:2017go,Stocker:2020dh`, along with various extensions. -The `tmodel` module -: Estimating plant allocation of gross primary productivity to growth and respiration, - using the [T Model](users/tmodel/tmodel) {cite:p}`Li:2014bc`. - -The `core` module -: Contains fundamental utilities and physics functionality shared across the - package, including the [hygro](users/hygro) and the utilities submodules. +The `splash` module +: Fits the [SPLASH v1 model](users/splash.md), which can be used to + estimate soil moisture, actual evapotranspiration and soil runoff from daily + temperature, precipitation and sunshine data {cite:p}`davis:2017a`. + +The `demography` module +: Provides functionality for [modelling plant allocation and growth and + demography](users/demography/module_overview.md), including classes to represent plant + functional types, cohorts and communities. This module includes an implementation of + the T Model for estimating plant allocation of gross primary productivity to growth + and respiration {cite:p}`Li:2014bc`. This module is still in active development but a + lot of initial functionality is present. ## Indices and tables diff --git a/docs/source/users/demography/canopy.md b/docs/source/users/demography/canopy.md index c606f0f5..ba4a31dd 100644 --- a/docs/source/users/demography/canopy.md +++ b/docs/source/users/demography/canopy.md @@ -43,7 +43,7 @@ from pyrealm.demography.flora import PlantFunctionalType, Flora from pyrealm.demography.community import Cohorts, Community from pyrealm.demography.crown import CrownProfile, get_crown_xy from pyrealm.demography.canopy import Canopy -from pyrealm.demography.t_model_functions import StemAllometry +from pyrealm.demography.tmodel import StemAllometry ``` The `canopy` module in `pyrealm` is used to calculate a model of the light environment diff --git a/docs/source/users/demography/crown.md b/docs/source/users/demography/crown.md index a2de5922..3e6a81b8 100644 --- a/docs/source/users/demography/crown.md +++ b/docs/source/users/demography/crown.md @@ -45,7 +45,7 @@ from pyrealm.demography.flora import ( Flora, ) -from pyrealm.demography.t_model_functions import ( +from pyrealm.demography.tmodel import ( calculate_dbh_from_height, StemAllometry, ) @@ -198,7 +198,7 @@ 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 +* a {class}`~pyrealm.demography.tmodel.StemAllometry` instance setting the specific stem sizes for the profile, and * a set of heights at which to estimate the profile variables @@ -226,7 +226,7 @@ flora.to_pandas()[["name", "h_max", "ca_ratio", "m", "n", "f_g", "q_m", "z_max_p 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 +{meth}`~pyrealm.demography.tmodel.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. diff --git a/docs/source/users/demography/t_model.md b/docs/source/users/demography/t_model.md index 44c0fd73..cb5708a8 100644 --- a/docs/source/users/demography/t_model.md +++ b/docs/source/users/demography/t_model.md @@ -36,7 +36,7 @@ import numpy as np import pandas as pd from pyrealm.demography.flora import PlantFunctionalType, Flora -from pyrealm.demography.t_model_functions import StemAllocation, StemAllometry +from pyrealm.demography.tmodel import StemAllocation, StemAllometry ``` To generate predictions under the T Model, we need a Flora object providing the @@ -55,7 +55,7 @@ flora = Flora([short_pft, medium_pft, tall_pft]) ## Stem allometry We can visualise how the stem size, canopy size and various masses of PFTs change with -stem diameter by using the {class}`~pyrealm.demography.t_model_functions.StemAllometry` +stem diameter by using the {class}`~pyrealm.demography.tmodel.StemAllometry` class. Creating a `StemAllometry` instance needs an existing `Flora` instance and an array of values for diameter at breast height (DBH, metres). The returned class contains the predictions of the T Model for: @@ -77,7 +77,7 @@ for each PFT. This then calculates a single estimate at the given size for each single_allometry = StemAllometry(stem_traits=flora, at_dbh=np.array([0.1, 0.1, 0.1])) ``` -The {meth}`~pyrealm.demography.t_model_functions.StemAllometry` class provides the +The {meth}`~pyrealm.demography.tmodel.StemAllometry` class provides the {meth}`~pyrealm.demography.core.PandasExporter.to_pandas()` method to export the stem data for data exploration. @@ -124,7 +124,7 @@ for ax, (var, ylab) in zip(axes.flatten(), plot_details): ``` The {meth}`~pyrealm.demography.core.PandasExporter.to_pandas()` method of the -{meth}`~pyrealm.demography.t_model_functions.StemAllometry` class can still be used, but +{meth}`~pyrealm.demography.tmodel.StemAllometry` class can still be used, but the values are stacked into columns along with a index showing the different cohorts. ```{code-cell} ipython3 @@ -135,7 +135,7 @@ allometries.to_pandas() The T Model also predicts how potential GPP will be allocated to respiration, turnover and growth for stems with a given PFT and allometry using the -{meth}`~pyrealm.demography.t_model_functions.StemAllometry` class. Again, a single +{meth}`~pyrealm.demography.tmodel.StemAllometry` class. Again, a single value can be provided to get a single estimate of the allocation model for each stem: ```{code-cell} ipython3 @@ -146,7 +146,7 @@ single_allocation ``` The {meth}`~pyrealm.demography.core.PandasExporter.to_pandas()` method of the -{meth}`~pyrealm.demography.t_model_functions.StemAllocation` class can be used to +{meth}`~pyrealm.demography.tmodel.StemAllocation` class can be used to export data for exploration. ```{code-cell} ipython3 @@ -229,7 +229,7 @@ fig.delaxes(axes[-1]) ``` As before, the {meth}`~pyrealm.demography.core.PandasExporter.to_pandas()` method of the -{meth}`~pyrealm.demography.t_model_functions.StemAllometry` classs can be used to export +{meth}`~pyrealm.demography.tmodel.StemAllometry` classs can be used to export the data for each stem: ```{code-cell} ipython3 diff --git a/docs/source/users/tmodel/canopy.md b/docs/source/users/tmodel/canopy.md deleted file mode 100644 index d760ca77..00000000 --- a/docs/source/users/tmodel/canopy.md +++ /dev/null @@ -1,688 +0,0 @@ ---- -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 -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 ---- - -# Canopy model - -This notebook walks through the steps in generating the canopy model as (hopefully) used -in Plant-FATE. - -## The T Model - -The T Model provides a numerical description of how tree geometry scales with stem -diameter, and an allocation model of how GPP predicts changes in stem diameter. - -The implementation in `pyrealm` provides a class representing a particular plant -functional type, using a set of traits. The code below creates a PFT with the default -set of trait values. - -```{warning} -This sketch: - -* Assumes a single individual of each stem diameter, but in practice - we are going to want to include a number of individuals to capture cohorts. -* Assumes a single PFT, where we will need to provide a mixed community. -* Consequently handles forest inventory properties in a muddled way: we will - likely package all of the stem data into a single class, probably a community - object. -``` - -```{code-cell} ipython3 -import numpy as np -import matplotlib.pyplot as plt -from scipy.optimize import root_scalar - -from pyrealm.tmodel import TTree - -np.set_printoptions(precision=3) - -# Plant functional type with default parameterization -pft = TTree(diameters=np.array([0.1, 0.15, 0.2, 0.25, 0.38, 0.4, 1.0])) -pft.traits -``` - -The scaling of a set of trees is automatically calculated using the initial diameters to -the `TTree` instance. This automatically calculates the other dimensions, such as -height, using the underlying scaling equations of the T Model. - -```{code-cell} ipython3 -pft.height -``` - -```{code-cell} ipython3 -:lines_to_next_cell: 2 - -pft.crown_area -``` - -+++ {"lines_to_next_cell": 2} - -### Crown shape - -Jaideep's extension of the T Model adds a crown shape model, driven by two parameters -($m$ and $n$) that provide a very flexible description of the vertical crown profile. -These are used to derive two further invariant parameters: $q_m$ scales the canopy -radius to match the predictions of crown area under the T model and $p_{zm}$ is the -proportion of the total stem height at which the maximum crown radius occurs. Both these -values are constant for a plant functional type. - -$$ -\begin{align} -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} \\[8pt] -p_{zm} &= \left(\dfrac{n-1}{m n -1}\right)^ {\tfrac{1}{n}} -\end{align} -$$ - -For individual stems, with expected height $H$ and crown area $A_c$, we can then -estimate: - -* the height $z_m$ at which the maximum crown radius $r_m$ is found, and -* a slope $r_0$ that scales the relative canopy radius so that the $r_m$ matches - the allometric prediction of $A_c$ from the T Model. - -$$ -\begin{align} -z_m &= H p_{zm}\\[8pt] -r_0 &= \frac{1}{q_m}\sqrt{\frac{A_c}{\pi}} -\end{align} -$$ - -```{code-cell} ipython3 -:lines_to_next_cell: 2 - -def calculate_qm(m, n): - - # Constant q_m - return ( - m - * n - * ((n - 1) / (m * n - 1)) ** (1 - 1 / n) - * (((m - 1) * n) / (m * n - 1)) ** (m - 1) - ) - - -def calculate_stem_canopy_factors(pft, m, n): - - # Height of maximum crown radius - zm = pft.height * ((n - 1) / (m * n - 1)) ** (1 / n) - - # Slope to give Ac at zm - r0 = 1 / qm * np.sqrt(pft.crown_area / np.pi) - - return zm, r0 - - -# Shape parameters for a fairly top heavy crown profile -m = 2 -n = 5 -qm = calculate_qm(m=m, n=n) -zm, r0 = calculate_stem_canopy_factors(pft=pft, m=m, n=n) - -print("qm = ", np.round(qm, 4)) -print("zm = ", zm) -print("r0 = ", r0) -``` - -+++ {"lines_to_next_cell": 2} - -The following functions then provide the value at height $z$ of relative $q(z)$ and -actual $r(z)$ canopy radius: - -$$ -\begin{align} -q(z) &= m n \left(\dfrac{z}{H}\right) ^ {n -1} - \left( 1 - \left(\dfrac{z}{H}\right) ^ n \right)^{m-1}\\[8pt] -r(z) &= r_0 \; q(z) -\end{align} -$$ - -```{code-cell} ipython3 -def calculate_relative_canopy_radius_at_z(z, H, m, n): - """Calculate q(z)""" - - z_over_H = z / H - - return m * n * z_over_H ** (n - 1) * (1 - z_over_H**n) ** (m - 1) -``` - -```{code-cell} ipython3 -# Validate that zm and r0 generate the predicted maximum crown area -q_zm = calculate_relative_canopy_radius_at_z(zm, pft.height, m, n) -rm = r0 * q_zm -print("rm = ", rm) -``` - -```{code-cell} ipython3 -np.allclose(rm**2 * np.pi, pft.crown_area) -``` - -Vertical crown radius profiles can now be calculated for each stem: - -```{code-cell} ipython3 -# Create an interpolation from ground to maximum stem height, with 5 cm resolution. -# Also append a set of values _fractionally_ less than the exact height of stems -# so that the height at the top of each stem is included but to avoid floating -# point issues with exact heights. - -zres = 0.05 -z = np.arange(0, pft.height.max() + 1, zres) -z = np.sort(np.concatenate([z, pft.height - 0.00001])) - -# Convert the heights into a column matrix to broadcast against the stems -# and then calculate r(z) = r0 * q(z) -rz = r0 * calculate_relative_canopy_radius_at_z(z[:, None], pft.height, m, n) - -# When z > H, rz < 0, so set radius to 0 where rz < 0 -rz[np.where(rz < 0)] = 0 - -np.cumsum(np.convolve(rm, np.ones(2), "valid") + 0.1) -``` - -Those can be plotted out to show the vertical crown radius profiles - -```{code-cell} ipython3 -# Separate the stems along the x axis for plotting -stem_x = np.concatenate( - [np.array([0]), np.cumsum(np.convolve(rm, np.ones(2), "valid") + 0.4)] -) - -# Get the canopy sections above and below zm -rz_above_zm = np.where(np.logical_and(rz > 0, np.greater.outer(z, zm)), rz, np.nan) -rz_below_zm = np.where(np.logical_and(rz > 0, np.less_equal.outer(z, zm)), rz, np.nan) - -# Plot the canopy parts -plt.plot(stem_x + rz_below_zm, z, color="khaki") -plt.plot(stem_x - rz_below_zm, z, color="khaki") -plt.plot(stem_x + rz_above_zm, z, color="forestgreen") -plt.plot(stem_x - rz_above_zm, z, color="forestgreen") - -# Add the maximum radius -plt.plot(np.vstack((stem_x - rm, stem_x + rm)), np.vstack((zm, zm)), color="firebrick") - -# Plot the stem centre lines -plt.vlines(stem_x, 0, pft.height, linestyles="-", color="grey") - -plt.gca().set_aspect("equal") -``` - -## Canopy structure - -The canopy structure model uses the perfect plasticity approximation (PPA), which -assumes that plants can arrange their canopies to fill the available space $A$. -It takes the **projected area of stems** $Ap(z)$ within the canopy and finds the heights -at which each canopy layer closes ($z^*_l$ for $l = 1, 2, 3 ...$) where the total projected -area of the canopy equals $lA$. - -### Canopy projected area - -The projected area $A_p$ for a stem with height $H$, a maximum crown area $A_c$ at a -height $z_m$ and $m$, $n$ and $q_m$ for the associated plant functional type 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} -$$ - -```{code-cell} ipython3 -Stems = float | np.ndarray - - -def calculate_projected_area( - z: float, - pft, - m: Stems, - n: Stems, - qm: Stems, - zm: Stems, -) -> np.ndarray: - """Calculate projected crown area above a given height. - - This function takes PFT specific parameters (shape parameters) and stem specific - sizes and estimates the projected crown area above a given height $z$. The inputs - can either be scalars describing a single stem or arrays representing a community - of stems. If only a single PFT is being modelled then `m`, `n`, `qm` and `fg` can - be scalars with arrays `H`, `Ac` and `zm` giving the sizes of stems within that - PFT. - - Args: - z: Canopy height - m, n, qm : PFT specific shape parameters - pft, qm, zm: stem data - - """ - - # Calculate q(z) - qz = calculate_relative_canopy_radius_at_z(z, pft.height, m, n) - - # Calculate Ap given z > zm - Ap = pft.crown_area * (qz / qm) ** 2 - # Set Ap = Ac where z <= zm - Ap = np.where(z <= zm, pft.crown_area, Ap) - # Set Ap = 0 where z > H - Ap = np.where(z > pft.height, 0, Ap) - - return Ap -``` - -The code below calculates the projected crown area for each stem and then plots the -vertical profile for individual stems and across the community. - -```{code-cell} ipython3 -:lines_to_next_cell: 2 - -# Calculate the projected area for each stem -Ap_z = calculate_projected_area(z=z[:, None], pft=pft, m=m, n=n, qm=qm, zm=zm) - -# Plot the calculated values for each stem and across the community -fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(10, 5)) - -# Plot the individual stem projected area -ax1.set_ylabel("Height above ground ($z$, m)") -ax1.plot(Ap_z, z) -ax1.set_xlabel("Stem $A_p(z)$ (m2)") - -# Plot the individual stem projected area -ax2.plot(np.nansum(Ap_z, axis=1), z) -ax2.set_xlabel("Total community $A_p(z)$ (m2)") - -plt.tight_layout() -``` - -+++ {"lines_to_next_cell": 2} - -### Canopy closure and canopy gap fraction - -The total cumulative projected area shown above 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. - -The definition of the height of canopy layer closure ($z^*_l$) for a given canopy -layer $l = 1, ..., l_m$ is then: - -$$ -\sum_1^{N_s}{ A_p(z^*_l)} = l A(1 - f_G) -$$ - -This can be found numerically using a root solver as: - -$$ -\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 -def solve_canopy_closure_height( - z: float, - l: int, - A: float, - fG: float, - m: Stems, - n: Stems, - qm: Stems, - pft: Stems, - zm: Stems, -) -> np.ndarray: - """Solver function for canopy closure height. - - This function returns the difference between the total community projected area - at a height $z$ and the total available canopy space for canopy layer $l$, given - the community gap fraction for a given height. It is used with a root solver to - find canopy layer closure heights $z^*_l* for a community. - - Args: - m, n, qm : PFT specific shape parameters - H, Ac, zm: stem specific sizes - A, l: cell area and layer index - fG: community gap fraction - """ - - # Calculate Ap(z) - Ap_z = calculate_projected_area(z=z, pft=pft, m=m, n=n, qm=qm, zm=zm) - - # Return the difference between the projected area and the available space - return Ap_z.sum() - (A * l) * (1 - fG) - - -def calculate_canopy_heights( - A: float, - fG: float, - m: Stems, - n: Stems, - qm: Stems, - pft, - zm: Stems, -): - - # Calculate the number of layers - total_community_ca = pft.crown_area.sum() - n_layers = int(np.ceil(total_community_ca / (A * (1 - fG)))) - - # Data store for z* - z_star = np.zeros(n_layers) - - # Loop over the layers TODO - edge case of completely filled final layer - for lyr in np.arange(n_layers - 1): - z_star[lyr] = root_scalar( - solve_canopy_closure_height, - args=(lyr + 1, A, fG, m, n, qm, pft, zm), - bracket=(0, pft.height.max()), - ).root - - return z_star -``` - -The example below calculates the projected crown area above ground level for the example -stems. These should be identical to the crown area of the stems. - -```{code-cell} ipython3 -# Set the total available canopy space and community gap fraction -canopy_area = 32 -community_gap_fraction = 2 / 32 - -z_star = calculate_canopy_heights( - A=canopy_area, fG=community_gap_fraction, m=m, n=n, qm=qm, pft=pft, zm=zm -) - -print("z_star = ", z_star) -``` - -We can now plot the canopy stems alongside the community $A_p(z)$ profile, and -superimpose the calculated $z^*_l$ values and the cumulative canopy area for each layer -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} ipython3 -community_Ap_z = np.nansum(Ap_z, axis=1) - -fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(10, 5)) - -zcol = "red" - -# Plot the canopy parts -ax1.plot(stem_x + rz_below_zm, z, color="khaki") -ax1.plot(stem_x - rz_below_zm, z, color="khaki") -ax1.plot(stem_x + rz_above_zm, z, color="forestgreen") -ax1.plot(stem_x - rz_above_zm, z, color="forestgreen") - -# Add the maximum radius -ax1.plot(np.vstack((stem_x - rm, stem_x + rm)), np.vstack((zm, zm)), color="firebrick") - -# Plot the stem centre lines -ax1.vlines(stem_x, 0, pft.height, linestyles="-", color="grey") - -ax1.set_ylabel("Height above ground ($z$, m)") -ax1.set_xlabel("Arbitrary stem position") -ax1.hlines(z_star, 0, (stem_x + rm).max(), color=zcol, linewidth=0.5) - -# Plot the projected community crown area by height along with the heights -# at which different canopy layers close -ax2.hlines(z_star, 0, community_Ap_z.max(), color=zcol, linewidth=0.5) -ax2.vlines( - canopy_area * np.arange(1, len(z_star)) * (1 - community_gap_fraction), - 0, - pft.height.max(), - color=zcol, - linewidth=0.5, -) - -ax2.plot(community_Ap_z, z) -ax2.set_xlabel("Projected crown area above height $z$ ($A_p(z)$, m2)") - -# Add z* values on the righthand axis -ax3 = ax2.twinx() - - -def z_star_labels(X): - return [f"$z^*_{l + 1}$ = {z:.2f}" for l, z in enumerate(X)] - - -ax3.set_ylim(ax2.get_ylim()) -ax3.set_yticks(z_star) -ax3.set_yticklabels(z_star_labels(z_star)) - -ax4 = ax2.twiny() - -# Add canopy layer closure areas on top axis -cum_area = np.arange(1, len(z_star)) * canopy_area * (1 - community_gap_fraction) - - -def cum_area_labels(X): - return [f"$A_{l + 1}$ = {z:.1f}" for l, z in enumerate(X)] - - -ax4.set_xlim(ax2.get_xlim()) -ax4.set_xticks(cum_area) -ax4.set_xticklabels(cum_area_labels(cum_area)) - -plt.tight_layout() -``` - -The projected area from individual stems to each canopy layer can then be calculated at -$z^*_l$ and hence the projected area of canopy **within each layer**. - -```{code-cell} ipython3 -# Calculate the canopy area above z_star for each stem -Ap_z_star = calculate_projected_area(z=z_star[:, None], pft=pft, m=m, n=n, qm=qm, zm=zm) - -print(Ap_z_star) -``` - -```{code-cell} ipython3 -:lines_to_next_cell: 2 - -# Calculate the contribution _within_ each layer per stem -Ap_within_layer = np.diff(Ap_z_star, axis=0, prepend=0) - -print(Ap_within_layer) -``` - -+++ {"lines_to_next_cell": 2} - -### Leaf area within canopy layers - -The projected area occupied by leaves at a given height $\tilde{A}_{cp}(z)$ is -needed to calculate light transmission through those layers. This differs from the -projected area $A_p(z)$ because, although a tree occupies an area in the canopy -following the PPA, a **crown gap fraction** ($f_g$) reduces the actual leaf area -at a given height $z$. - -The crown gap fraction does not affect the overall projected canopy area at ground -level or the community gap fraction: the amount of clear sky at ground level is -governed purely by $f_G$. Instead it models how leaf gaps in the upper canopy are -filled by leaf area at lower heights. It captures the vertical distribution of -leaf area within the canopy: a higher $f_g$ will give fewer leaves at the top of -the canopy and more leaves further down within the canopy. - -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} -$$ - -The function below calculates $\tilde{A}_{cp}(z)$. - -```{code-cell} ipython3 -def calculate_leaf_area( - z: float, - fg: float, - pft, - m: Stems, - n: Stems, - qm: Stems, - zm: Stems, -) -> np.ndarray: - """Calculate leaf area above a given height. - - This function takes PFT specific parameters (shape parameters) and stem specific - sizes and estimates the projected crown area above a given height $z$. The inputs - can either be scalars describing a single stem or arrays representing a community - of stems. If only a single PFT is being modelled then `m`, `n`, `qm` and `fg` can - be scalars with arrays `H`, `Ac` and `zm` giving the sizes of stems within that - PFT. - - Args: - z: Canopy height - fg: crown gap fraction - m, n, qm : PFT specific shape parameters - pft, qm, zm: stem data - """ - - # Calculate q(z) - qz = calculate_relative_canopy_radius_at_z(z, pft.height, m, n) - - # Calculate Ac term - Ac_term = pft.crown_area * (qz / qm) ** 2 - # Set Acp either side of zm - Acp = np.where(z <= zm, pft.crown_area - Ac_term * fg, Ac_term * (1 - fg)) - # Set Ap = 0 where z > H - Acp = np.where(z > pft.height, 0, Acp) - - return Acp -``` - -The plot below shows how the vertical leaf area profile for the community changes for -different values of $f_g$. When $f_g = 0$, then $A_cp(z) = A_p(z)$ (red line) because -there are no crown gaps and hence all of the leaf area is within the crown surface. As -$f_g \to 1$, more of the leaf area is displaced deeper into the canopy, leaves in the -lower crown intercepting light coming through holes in the upper canopy. - -```{code-cell} ipython3 -fig, ax1 = plt.subplots(1, 1, figsize=(6, 5)) - -for fg in np.arange(0, 1.01, 0.05): - - if fg == 0: - color = "red" - label = "$f_g = 0$" - lwd = 0.5 - elif fg == 1: - color = "blue" - label = "$f_g = 1$" - lwd = 0.5 - else: - color = "black" - label = None - lwd = 0.25 - - Acp_z = calculate_leaf_area(z=z[:, None], fg=fg, pft=pft, m=m, n=n, qm=qm, zm=zm) - ax1.plot(np.nansum(Acp_z, axis=1), z, color=color, linewidth=lwd, label=label) - -ax1.set_xlabel(r"Projected leaf area above height $z$ ($\tilde{A}_{cp}(z)$, m2)") -ax1.legend(frameon=False) -``` - -We can now calculate the crown area occupied by leaves above the height of each closed -layer $z^*_l$: - -```{code-cell} ipython3 -# Calculate the leaf area above z_star for each stem -crown_gap_fraction = 0.05 -Acp_z_star = calculate_leaf_area( - z=z_star[:, None], fg=crown_gap_fraction, pft=pft, m=m, n=n, qm=qm, zm=zm -) - -print(Acp_z_star) -``` - -And from that, the area occupied by leaves **within each layer**. These values are -similar to the projected crown area within layers (`Ap_within_layer`, above) but -leaf area is displaced into lower layers because $f_g > 0$. - -```{code-cell} ipython3 -# Calculate the contribution _within_ each layer per stem -Acp_within_layer = np.diff(Acp_z_star, axis=0, prepend=0) - -print(Acp_within_layer) -``` - -### 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) - -print("LAI stem = ", LAI_stem) -print("LAI layer = ", LAI_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))) -``` - -## 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. diff --git a/docs/source/users/tmodel/tmodel.md b/docs/source/users/tmodel/tmodel.md deleted file mode 100644 index d7e75b31..00000000 --- a/docs/source/users/tmodel/tmodel.md +++ /dev/null @@ -1,252 +0,0 @@ ---- -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 - 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 T Model - -This module provides a Python implementation of the T-Model {cite:p}`Li:2014bc`, which -provides a physiological model of tree growth given a set of traits on tree growth -scalings and allocation of primary production. - -The module uses three key components of the T Model: - -1. the {class}`~pyrealm.constants.tmodel_const.TModelTraits` class, -2. the {class}`~pyrealm.tmodel.TTree` class, and -3. the {func}`~pyrealm.tmodel.grow_ttree` function. - -## The {class}`~pyrealm.constants.tmodel_const.TModelTraits` class - -The T Model depends on a set of 14 traits that are used to describe the geometry of a -tree and the allocation of carbon within the tree and which are listed in the linked -class description. - -The class can be used to create a default T Model trait set: - -```{code-cell} ipython3 -import numpy as np -from pyrealm import tmodel - -# A tree using the default parameterisation -traits1 = tmodel.TModelTraits() -print(traits1) -``` - -It can also be edited to generate different growth patterns: - -```{code-cell} ipython3 -# A slower growing tree with a higher maximum height -traits2 = tmodel.TModelTraits(a_hd=50, h_max=40) -print(traits2) -``` - -## The {class}`~pyrealm.tmodel.TTree` class - -This class implements the mathematical description of tree growth under the T Model. - -There are three stages to generating predictions: - -1. initialising a model, -2. setting the stem diameters to calculate tree geometry and mass, and -3. calculating growth predictions for a given estimate of gross primary productivity (GPP). - -Note that the {class}`~pyrealm.tmodel.TTree` is not used to model tree growth through -time (see below and {func}`~pyrealm.tmodel.grow_ttree`). It simply calculates the -predictions of the T Model for trees with a given diameter and GPP. Because it accepts -**arrays** of data, it can be used to very quickly visualise the behaviour of the TModel -for a given set of traits and diameters. - -### Initialising a {class}`~pyrealm.tmodel.TTree` object - -A {class}`~pyrealm.tmodel.TTree` object is created by providing an array of initial stem -diameters and an optional set of traits as a -{class}`~pyrealm.constants.tmodel_const.TModelTraits` object. If no traits are provided, -the default {class}`~pyrealm.constants.tmodel_const.TModelTraits` settings are used. - -```{code-cell} ipython3 -# Use a sequence of diameters from sapling to large tree -diameters = np.linspace(0.02, 2, 100) -tree1 = tmodel.TTree(diameters=diameters) # Using default traits -tree2 = tmodel.TTree(diameters=diameters, traits=traits2) -``` - -These inputs are then immediately used to calculate the following properties of the -{class}`~pyrealm.tmodel.TTree` that scale simply with tree diameter: - -* Stem diameter (`diameter`) -* Stem height (`height`) -* Crown fraction (`crown_fraction`) -* Crown area (`crown_area`) -* Mass of stem (`mass_stm`) -* Mass of foliage and fine roots (`mass_fol`) -* Mass of sapwood (`mass_swd`) - -Using an array of diameter values provides an immediate way to visualise the geometric -scaling resulting from a particular set of plant traits: - -```{code-cell} ipython3 -:tags: [hide-input] - -from matplotlib import pyplot - -fig, (ax1, ax2, ax3) = pyplot.subplots(1, 3, figsize=(12, 4)) -ax1.plot(tree1.diameter, tree1.height, label="traits1") -ax1.plot(tree2.diameter, tree2.height, label="traits2") -ax1.set_title("Height") -ax1.set_xlabel("Stem diameter (m)") -ax1.set_ylabel("Height (m)") -ax1.legend() -ax2.plot(tree1.diameter, tree1.crown_area, label="traits1") -ax2.plot(tree2.diameter, tree2.crown_area, label="traits2") -ax2.set_title("Crown Area") -ax2.set_xlabel("Stem diameter (m)") -ax2.set_ylabel("Crown area (m2)") -ax2.legend() -ax3.plot(tree1.diameter, tree1.mass_stm, label="Stem") -ax3.plot(tree1.diameter, tree1.mass_swd, label="Sapwood") -ax3.plot(tree1.diameter, tree1.mass_fol, label="Leaf and fine root") -ax3.set_yscale("log") -ax3.set_title("Mass (traits1)") -ax3.set_xlabel("Stem diameter (m)") -ax3.set_ylabel("Log Mass (kg)") -ax3.legend() -pyplot.show() -``` - -### Calculating growth - -The {meth}`~pyrealm.tmodel.TTree.calculate_growth` method is used to provide estimates -of gross primary productivity (GPP) to a {class}`~pyrealm.tmodel.TTree` instance in -order to apply the T Model predictions of GPP allocation to plant respiration and -growth. - -Running the method populates properties of the {class}`~pyrealm.tmodel.TTree` that -provide estimates of the following growth parameters: - -* Gross primary productivity (`gpp_actual`) -* Net primary productivity (`npp`) -* Sapwood respiration (`resp_swd`) -* Fine root respiration (`resp_frt`) -* Foliage maintenance respiration (`resp_fol`) -* Foliage and fine root turnover (`turnover`) -* Diameter increment (`delta_d`) -* Stem mass increment (`delta_mass_stm`) -* Fine root mass increment (`delta_mass_frt`) - -The code below calculates growth estimates at each diameter under a constant GPP of 7 -TODO - UNITS!. - -```{code-cell} ipython3 -tree1.calculate_growth(np.array([7])) -tree2.calculate_growth(np.array([7])) -``` - -```{code-cell} ipython3 -:tags: [hide-input] - -fig, (ax1, ax2, ax3) = pyplot.subplots(1, 3, figsize=(12, 4)) -ax1.plot(tree1.diameter, tree1.npp, label="traits1") -ax1.plot(tree2.diameter, tree2.npp, label="traits2") -ax1.set_title("NPP") -ax1.set_xlabel("Net primary productivity (m)") -ax1.set_ylabel("Height (m)") -ax1.legend() -ax2.plot(tree1.diameter, tree1.resp_swd, label="Sapwood") -ax2.plot(tree1.diameter, tree1.resp_frt, label="Fine roots") -ax2.plot(tree1.diameter, tree1.resp_fol, label="Foliage") -ax2.set_title("Respiration (traits1)") -ax2.set_xlabel("Stem diameter (m)") -ax2.set_ylabel("Respiration") -ax2.legend() -ax3.plot(tree1.diameter, tree1.delta_d, label="traits1") -ax3.plot(tree2.diameter, tree2.delta_d, label="traits2") -ax3.set_title("Delta diameter ") -ax3.set_xlabel("Stem diameter (m)") -ax3.set_ylabel("Delta diameter (m)") -ax3.legend() -pyplot.show() -``` - -### Resetting stem diameters - -The {meth}`~pyrealm.tmodel.TTree.reset_diameters` can be used to update an existing -{class}`~pyrealm.tmodel.TTree` object with a new set of diameters using the same -{class}`~pyrealm.constants.tmodel_const.TModelTraits` definition. Using -{meth}`~pyrealm.tmodel.TTree.reset_diameters` automatically resets any calculated growth -parameters: they will need to be recalculated for the new diameters. - -```{code-cell} ipython3 -tree1.reset_diameters(np.array([0.0001])) -print(tree1.height) -``` - - diff --git a/pyrealm/demography/community.py b/pyrealm/demography/community.py index 14d1aa65..c7a97949 100644 --- a/pyrealm/demography/community.py +++ b/pyrealm/demography/community.py @@ -20,7 +20,7 @@ >>> import pandas as pd >>> >>> from pyrealm.demography.flora import PlantFunctionalType, Flora ->>> from pyrealm.demography.t_model_functions import ( +>>> from pyrealm.demography.tmodel import ( ... calculate_heights, calculate_crown_areas, calculate_stem_masses, ... calculate_foliage_masses ... ) @@ -139,7 +139,7 @@ from pyrealm.core.utilities import check_input_shapes from pyrealm.demography.core import CohortMethods, PandasExporter from pyrealm.demography.flora import Flora, StemTraits -from pyrealm.demography.t_model_functions import StemAllometry +from pyrealm.demography.tmodel import StemAllometry if sys.version_info[:2] >= (3, 11): import tomllib diff --git a/pyrealm/demography/crown.py b/pyrealm/demography/crown.py index e1aa1083..0ea0e83e 100644 --- a/pyrealm/demography/crown.py +++ b/pyrealm/demography/crown.py @@ -13,7 +13,7 @@ _validate_demography_array_arguments, ) from pyrealm.demography.flora import Flora, StemTraits -from pyrealm.demography.t_model_functions import StemAllometry +from pyrealm.demography.tmodel import StemAllometry def calculate_relative_crown_radius_at_z( diff --git a/pyrealm/demography/t_model_functions.py b/pyrealm/demography/tmodel.py similarity index 99% rename from pyrealm/demography/t_model_functions.py rename to pyrealm/demography/tmodel.py index e8a5d803..176af87f 100644 --- a/pyrealm/demography/t_model_functions.py +++ b/pyrealm/demography/tmodel.py @@ -62,7 +62,7 @@ def calculate_dbh_from_height( This function inverts the normal calculation of stem height (:math:`H`) from diameter at breast height (DBH, :math:`D`) in the T Model (see - :meth:`~pyrealm.demography.t_model_functions.calculate_heights`). This is a helper + :meth:`~pyrealm.demography.tmodel.calculate_heights`). This is a helper function to allow users to convert known stem heights for a plant functional type, with maximum height (:math:`H_{m}`) and initial slope of the height/diameter relationship (:math:`a`) into the expected DBH values. @@ -859,7 +859,7 @@ class StemAllocation(PandasExporter): :class:`~pyrealm.demography.flora.StemTraits`, providing plant functional trait data for a set of stems. stem_allometry: An instance of - :class:`~pyrealm.demography.t_model_functions.StemAllometry` + :class:`~pyrealm.demography.tmodel.StemAllometry` providing the stem size data for which to calculate allocation. at_potential_gpp: An array of potential GPP values at which to predict stem allocation. @@ -885,7 +885,7 @@ class StemAllocation(PandasExporter): :class:`~pyrealm.demography.flora.StemTraits`, providing plant functional trait data for a set of stems.""" stem_allometry: InitVar[StemAllometry] - """An instance of :class:`~pyrealm.demography.t_model_functions.StemAllometry` + """An instance of :class:`~pyrealm.demography.tmodel.StemAllometry` providing the stem size data for which to calculate allocation.""" at_potential_gpp: InitVar[NDArray[np.float64]] """An array of potential gross primary productivity for each stem that should be diff --git a/pyrealm/tmodel.py b/pyrealm/tmodel.py deleted file mode 100644 index 416c53c5..00000000 --- a/pyrealm/tmodel.py +++ /dev/null @@ -1,409 +0,0 @@ -"""The :mod:`~pyrealm.tmodel` module provides an implementation of the T Model of plant -growth given an estimate of gross primary productivity (GPP). - -* The growth form and productivity allocation model of a plant is set using - :class:`~pyrealm.constants.tmodel_const.TModelTraits`. -* The class :class:`~pyrealm.tmodel.TTree` is used to generate an instance of a - plant to be simulated, with methods :meth:`~pyrealm.tmodel.TTree.reset_diameters` - and :meth:`~pyrealm.tmodel.TTree.calculate_growth` to calculate the plant - geometry for a given diameter or predict growth from estimated GPP. -* The function :func:`~pyrealm.tmodel.grow_ttree` predicts plant growth through - time given a time series of GPP. -""" # noqa: D205 - -import numpy as np -from numpy.typing import NDArray - -from pyrealm.constants.tmodel_const import TModelTraits - -# Design Notes: -# -# One functionally easy thing to do the TTree object is to expose the geometry -# methods (e.g. by having TTree.calculate_crown_area()) to allow users to -# subclass the object and substitute their own geometry. Which is neat but the -# geometry and growth are not separate and there is no guarantee that a -# user-supplied geometry behaves as expected. So, easy to implement and a -# natural solution to people wanting to tinker with the geometry but not going -# there now. - - -class TTree: - """Model plant growth using the T model. - - This class provides an implementation of the calculations of tree geometry, mass and - growth described by :cite:t:`Li:2014bc`. All of the properties of the T model are - derived from a set of traits (see - :class:`~pyrealm.constants.tmodel_const.TModelTraits`), stem diameter measurements - and estimates of gross primary productivity. - - See the details of :meth:`~pyrealm.tmodel.TTree.reset_diameters` and - :meth:`~pyrealm.tmodel.TTree.calculate_growth` for details of the properties and - calculations. - - Args: - traits: An object of class :class:`~pyrealm.constants.tmodel_const.TModelTraits` - diameters: A float or np.array of stem diameters. - """ - - def __init__( - self, - diameters: NDArray[np.float64], - traits: TModelTraits = TModelTraits(), - ) -> None: - self.traits: TModelTraits = traits - - # The diameter is used to define all of the geometric scaling - # based on the trait parameters. - - self._diameter: NDArray[np.float64] - self._height: NDArray[np.float64] - self._crown_fraction: NDArray[np.float64] - self._crown_area: NDArray[np.float64] - self._mass_stm: NDArray[np.float64] - self._mass_fol: NDArray[np.float64] - self._mass_swd: NDArray[np.float64] - - self.reset_diameters(diameters) - - # Growth is then applied by providing estimated gpp using the - # calculate_growth() method, which populates the following: - self.growth_calculated: bool = False - self._gpp_raw: NDArray[np.float64] - self._gpp_actual: NDArray[np.float64] - self._npp: NDArray[np.float64] - self._resp_swd: NDArray[np.float64] - self._resp_frt: NDArray[np.float64] - self._resp_fol: NDArray[np.float64] - self._turnover: NDArray[np.float64] - self._d_mass_s: NDArray[np.float64] - self._d_mass_fr: NDArray[np.float64] - self._delta_d: NDArray[np.float64] - self._delta_mass_stm: NDArray[np.float64] - self._delta_mass_frt: NDArray[np.float64] - - def _check_growth_calculated(self, property: str) -> NDArray[np.float64]: - """Helper function to return growth values if calculated. - - This acts as a gatekeeper to make sure that a growth property is not returned - before calculate_growth() has been run on the current diameters. - - Args: - property: The property value to return if available. - """ - if not self.growth_calculated: - raise RuntimeError("Growth estimates not calculated: use calculate_growth") - - return getattr(self, property) - - @property - def diameter(self) -> NDArray[np.float64]: - """Individual diameter (m).""" - return self._diameter - - @property - def height(self) -> NDArray[np.float64]: - """Individual height (m).""" - return self._height - - @property - def crown_fraction(self) -> NDArray[np.float64]: - """Individual crown fraction (unitless).""" - return self._crown_fraction - - @property - def crown_area(self) -> NDArray[np.float64]: - """Individual crown area (m2).""" - return self._crown_area - - @property - def mass_swd(self) -> NDArray[np.float64]: - """Individual softwood mass (kg).""" - return self._mass_swd - - @property - def mass_stm(self) -> NDArray[np.float64]: - """Individual stem mass (kg).""" - return self._mass_stm - - @property - def mass_fol(self) -> NDArray[np.float64]: - """Individual foliage mass (kg).""" - return self._mass_fol - - @property - def gpp_raw(self) -> NDArray[np.float64]: - """Raw gross primary productivity.""" - return self._check_growth_calculated("_gpp_raw") - - @property - def gpp_actual(self) -> NDArray[np.float64]: - """Actual gross primary productivity.""" - return self._check_growth_calculated("_gpp_actual") - - @property - def resp_swd(self) -> NDArray[np.float64]: - """Individual softwood respiration costs.""" - return self._check_growth_calculated("_resp_swd") - - @property - def resp_frt(self) -> NDArray[np.float64]: - """Individual fine root respiration costs.""" - return self._check_growth_calculated("_resp_frt") - - @property - def resp_fol(self) -> NDArray[np.float64]: - """Individual foliar respiration costs.""" - return self._check_growth_calculated("_resp_fol") - - @property - def npp(self) -> NDArray[np.float64]: - """Net primary productivity.""" - return self._check_growth_calculated("_npp") - - @property - def turnover(self) -> NDArray[np.float64]: - """Plant turnover.""" - return self._check_growth_calculated("_turnover") - - @property - def d_mass_s(self) -> NDArray[np.float64]: - """Individual relative change in mass.""" - return self._check_growth_calculated("_d_mass_s") - - @property - def d_mass_fr(self) -> NDArray[np.float64]: - """Individual relative change in fine root mass.""" - return self._check_growth_calculated("_d_mass_fr") - - @property - def delta_d(self) -> NDArray[np.float64]: - """Individual change in diameter.""" - return self._check_growth_calculated("_delta_d") - - @property - def delta_mass_stm(self) -> NDArray[np.float64]: - """Individual total change in stem mass.""" - return self._check_growth_calculated("_delta_mass_stm") - - @property - def delta_mass_frt(self) -> NDArray[np.float64]: - """Individual total change in fine root mass.""" - return self._check_growth_calculated("_delta_mass_frt") - - def reset_diameters(self, values: NDArray[np.float64]) -> None: - """Reset the stem diameters for the T model. - - The set_diameter method can be used to reset the diameter values and then uses - these values to populate the geometric and mass properties that scale with stem - diameter. - - * Height (m, ``height``, :math:`H`): - - .. math:: - - H = H_{max} ( 1 - e^{a D/ H_{max}} - - TODO: complete this description. - """ - - self._diameter = values - - # Height of tree from diameter, Equation (4) of Li ea. - self._height = self.traits.h_max * ( - 1 - np.exp(-self.traits.a_hd * self._diameter / self.traits.h_max) - ) - - # Crown area of tree, Equation (8) of Li ea. - self._crown_area = ( - ((np.pi * self.traits.ca_ratio) / (4 * self.traits.a_hd)) - * self._diameter - * self._height - ) - - # Crown fraction, Equation (11) of Li ea. - self._crown_fraction = self._height / (self.traits.a_hd * self._diameter) - - # Masses - self._mass_stm = ( - (np.pi / 8) * (self._diameter**2) * self._height * self.traits.rho_s - ) - self._mass_fol = self._crown_area * self.traits.lai * (1 / self.traits.sla) - self._mass_swd = ( - self._crown_area - * self.traits.rho_s - * self._height - * (1 - self._crown_fraction / 2) - / self.traits.ca_ratio - ) - - # Flag any calculated growth values as outdated - self.growth_calculated = False - - def calculate_growth(self, gpp: NDArray[np.float64]) -> None: - """Calculate growth predictions given a GPP estimate. - - This method updates the instance with predicted changes in tree geometry, mass - and respiration costs from the initial state given an estimate of gross primary - productivity (GPP). - - Args: - gpp: Primary productivity - """ - - # GPP fixed per m2 of crown - self._gpp_raw = gpp - gpp_unit_cr = self._gpp_raw * ( - 1 - np.exp(-(self.traits.par_ext * self.traits.lai)) - ) - self._gpp_actual = self._crown_area * gpp_unit_cr - - # Respiration costs (Eqn 13 of Li ea) - # - sapwood, fine root and foliage maintenance - self._resp_swd = self._mass_swd * self.traits.resp_s - self._resp_frt = ( - self.traits.zeta * self.traits.sla * self._mass_fol * self.traits.resp_r - ) - self._resp_fol = self._gpp_actual * self.traits.resp_f - - # Net primary productivity - self._npp = self.traits.yld * ( - self._gpp_actual - self._resp_fol - self._resp_frt - self._resp_swd - ) - - # Turnover costs for foliage and fine roots - self._turnover = ( - self._crown_area - * self.traits.lai - * ( - (1 / (self.traits.sla * self.traits.tau_f)) - + (self.traits.zeta / self.traits.tau_r) - ) - ) - - # relative increments - these are used to calculate delta_d and - # then scaled by delta_d to give actual increments - self._d_mass_s = ( - np.pi - / 8 - * self.traits.rho_s - * self._diameter - * ( - self.traits.a_hd - * self._diameter - * (1 - (self._height / self.traits.h_max)) - + 2 * self._height - ) - ) - - self._d_mass_fr = ( - self.traits.lai - * ((np.pi * self.traits.ca_ratio) / (4 * self.traits.a_hd)) - * ( - self.traits.a_hd - * self._diameter - * (1 - self._height / self.traits.h_max) - + self._height - ) - * (1 / self.traits.sla + self.traits.zeta) - ) - - # Actual increments - self._delta_d = (self._npp - self._turnover) / ( - self._d_mass_s + self._d_mass_fr - ) - self._delta_mass_stm = self._d_mass_s * self._delta_d - self._delta_mass_frt = self._d_mass_fr * self._delta_d - - self.growth_calculated = True - - -def grow_ttree( - gpp: NDArray[np.float64], - d_init: NDArray[np.float64], - time_axis: int, - traits: TModelTraits = TModelTraits(), - outvars: tuple[str, ...] = ("diameter", "height", "crown_area", "delta_d"), -) -> dict[str, np.ndarray]: - """Fit a growth time series using the T Model.""" - raise NotImplementedError() - - -# This function fits the T Model incrementally to a set of modelled plants, -# given a time series of GPP estimates. - -# Args: -# gpp: An array of GPP values -# d_init: An array of starting diameters -# traits: Traits to be used -# time_axis: An axis in P0 and d that represents an annual time series -# outvars: A list of Tree properties to store. - -# Returns: -# A dictionary of estimates for the time series, exporting the values -# specified in `outvars` -# """ - -# # The gpp array should contain a time axis to increment over. The d_init input -# # should then be the same shape as P0 _without_ that time axis dimension -# # and the loop will iterate over the years - -# gpp_shape = gpp.shape -# if time_axis < 0 or time_axis >= len(gpp_shape): -# raise RuntimeError(f"time_axis must be >= 0 and <= {len(gpp_shape) - 1}") - -# # Check that the input shapes for a single year match the shape of the -# # initial diameters -# single_year_gpp = np.take(gpp, indices=0, axis=time_axis) -# _ = check_input_shapes(single_year_gpp, d_init) - -# # TODO: - handle 1D GPP time series applied to more than one diameter - -# # Initialise the Tree object -# tree = TTree(d_init, traits) - -# # Check the requested outvars -# if "diameter" not in outvars: -# raise RuntimeWarning("output_vars must include diameter") - -# badvars = [] -# for var in outvars: - -# try: -# _ = getattr(tree, var) -# except AttributeError: -# badvars.append(var) - -# if badvars: -# raise RuntimeError( -# f"Unknown tree properties in outvars: {', '.join(badvars)}" -# ) - -# # Create an array to store the requested variables by adding a dimension -# # to the gpp input with length set to the number of variables -# output = {v: np.zeros(gpp.shape) for v in outvars} - -# # Create an indexing object to insert values into the output. This is -# # a bit obscure: the inputs have a time axis and an arbitrary number -# # of other dimensions and the output adds another dimension for the -# # output variables. So this object creates a slice (:) on _all_ dimensions -# # and the loop then replaces the time axis and variable axis with integers. -# output_index = [slice(None)] * gpp.ndim - -# # Loop over the gpp time axis -# for year in np.arange(gpp_shape[time_axis]): - -# # Calculate the growth based on the current year of gpp -# tree.calculate_growth(np.take(gpp, indices=year, axis=time_axis)) - -# # Store the requested variables into the output arrays -# for each_var in outvars: - -# # Extract variable values from tree into output - set the last index -# # (variable dimension) to the variable index and the time axis to the year -# output_index[time_axis] = year -# output[each_var][tuple(output_index)] = getattr(tree, each_var) - -# # Now update the tree object -# tree.reset_diameters(tree.diameter + getattr(tree, "delta_d")) - -# return output diff --git a/tests/regression/demography/test_t_model_functions_against_rtmodel.py b/tests/regression/demography/test_t_model_functions_against_rtmodel.py index d66550e3..649ecbc8 100644 --- a/tests/regression/demography/test_t_model_functions_against_rtmodel.py +++ b/tests/regression/demography/test_t_model_functions_against_rtmodel.py @@ -114,7 +114,7 @@ def rvalues(): def test_calculate_heights(rvalues): """Test calculation of heights of tree from diameter.""" - from pyrealm.demography.t_model_functions import calculate_heights + from pyrealm.demography.tmodel import calculate_heights for pft, _, data in rvalues: actual_heights = calculate_heights( @@ -126,7 +126,7 @@ def test_calculate_heights(rvalues): def test_calculate_crown_areas(rvalues): """Tests calculation of crown areas of trees.""" - from pyrealm.demography.t_model_functions import calculate_crown_areas + from pyrealm.demography.tmodel import calculate_crown_areas for pft, _, data in rvalues: actual_crown_areas = calculate_crown_areas( @@ -142,7 +142,7 @@ def test_calculate_crown_areas(rvalues): def test_calculate_crown_fractions(rvalues): """Tests calculation of crown fractions of trees.""" - from pyrealm.demography.t_model_functions import calculate_crown_fractions + from pyrealm.demography.tmodel import calculate_crown_fractions for pft, _, data in rvalues: actual_crown_fractions = calculate_crown_fractions( @@ -158,7 +158,7 @@ def test_calculate_crown_fractions(rvalues): def test_calculate_stem_masses(rvalues): """Tests happy path for calculation of stem masses of trees.""" - from pyrealm.demography.t_model_functions import calculate_stem_masses + from pyrealm.demography.tmodel import calculate_stem_masses for pft, _, data in rvalues: actual_stem_masses = calculate_stem_masses( @@ -172,7 +172,7 @@ def test_calculate_stem_masses(rvalues): def test_calculate_foliage_masses(rvalues): """Tests calculation of foliage masses of trees.""" - from pyrealm.demography.t_model_functions import calculate_foliage_masses + from pyrealm.demography.tmodel import calculate_foliage_masses for pft, _, data in rvalues: actual_foliage_masses = calculate_foliage_masses( @@ -184,7 +184,7 @@ def test_calculate_foliage_masses(rvalues): def test_calculate_sapwood_masses(rvalues): """Tests calculation of sapwood masses of trees.""" - from pyrealm.demography.t_model_functions import calculate_sapwood_masses + from pyrealm.demography.tmodel import calculate_sapwood_masses for pft, _, data in rvalues: actual_sapwood_masses = calculate_sapwood_masses( @@ -204,7 +204,7 @@ def test_calculate_whole_crown_gpp(rvalues): foliar respiration. """ - from pyrealm.demography.t_model_functions import calculate_whole_crown_gpp + from pyrealm.demography.tmodel import calculate_whole_crown_gpp for pft, _, data in rvalues: actual_whole_crown_gpp = calculate_whole_crown_gpp( @@ -219,7 +219,7 @@ def test_calculate_whole_crown_gpp(rvalues): def test_calculate_sapwood_respiration(rvalues): """Tests calculation of sapwood respiration of trees.""" - from pyrealm.demography.t_model_functions import calculate_sapwood_respiration + from pyrealm.demography.tmodel import calculate_sapwood_respiration for pft, _, data in rvalues: actual_sapwood_respiration = calculate_sapwood_respiration( @@ -239,7 +239,7 @@ def test_calculate_foliar_respiration(rvalues): circular but is important to validate this difference. """ - from pyrealm.demography.t_model_functions import calculate_foliar_respiration + from pyrealm.demography.tmodel import calculate_foliar_respiration for pft, _, data in rvalues: actual_foliar_respiration = calculate_foliar_respiration( @@ -256,7 +256,7 @@ def test_calculate_foliar_respiration(rvalues): def test_calculate_fine_root_respiration(rvalues): """Tests calculation of fine root respiration of trees.""" - from pyrealm.demography.t_model_functions import calculate_fine_root_respiration + from pyrealm.demography.tmodel import calculate_fine_root_respiration for pft, _, data in rvalues: actual_fine_root_respiration = calculate_fine_root_respiration( @@ -279,7 +279,7 @@ def test_calculate_net_primary_productivity(rvalues): respiration from potential GPP before calculating crown GPP. """ - from pyrealm.demography.t_model_functions import calculate_net_primary_productivity + from pyrealm.demography.tmodel import calculate_net_primary_productivity for pft, _, data in rvalues: actual_npp = calculate_net_primary_productivity( @@ -299,7 +299,7 @@ def test_calculate_net_primary_productivity(rvalues): def test_calculate_foliage_and_fine_root_turnover(rvalues): """Tests calculation of fine root respiration of trees.""" - from pyrealm.demography.t_model_functions import ( + from pyrealm.demography.tmodel import ( calculate_foliage_and_fine_root_turnover, ) @@ -321,7 +321,7 @@ def test_calculate_foliage_and_fine_root_turnover(rvalues): def test_calculate_growth_increments(rvalues): """Tests calculation of fine root respiration of trees.""" - from pyrealm.demography.t_model_functions import ( + from pyrealm.demography.tmodel import ( calculate_growth_increments, ) diff --git a/tests/regression/tmodel/test_tmodel.py b/tests/regression/tmodel/test_tmodel.py deleted file mode 100644 index 2504399a..00000000 --- a/tests/regression/tmodel/test_tmodel.py +++ /dev/null @@ -1,248 +0,0 @@ -"""Test TModel class. - -Tests the init, grow_ttree and other methods of TModel. -""" - -from contextlib import nullcontext as does_not_raise -from importlib import resources - -import numpy as np -import pandas as pd -import pytest - -# Fixtures: inputs and expected values from the original implementation in R - - -@pytest.fixture(scope="module") -def rvalues(): - """Fixture to load test inputs from file. - - This is a time series of growth using the default trait values in R, mapped to the - internal property names used in TTree - """ - from pyrealm.tmodel import TModelTraits - - datapath = ( - resources.files("pyrealm_build_data.t_model") / "rtmodel_output_default.csv" - ) - - data = pd.read_csv(datapath) - - data = data.rename( - columns={ - "dD": "delta_d", - "D": "diameter", - "H": "height", - "fc": "crown_fraction", - "Ac": "crown_area", - "Wf": "mass_fol", - "Ws": "mass_stm", - "Wss": "mass_swd", - "P0": "potential_gpp", - "GPP": "gpp_actual", - "Rm1": "resp_swd", - "Rm2": "resp_frt", - "dWs": "delta_mass_stm", - "dWfr": "delta_mass_frt", - } - ) - - # Fix some scaling differences: - # The R tmodel implementation rescales reported delta_d as a radial increase in - # millimetres, not diameter increase in metres - data["delta_d"] = data["delta_d"] / 500 - - # The R tmodel implementation slices off foliar respiration costs from - # GPP before doing anything - the pyrealm.tmodel implementation keeps - # this cost within the tree calculation - traits = TModelTraits() - data["gpp_actual"] = data["gpp_actual"] / (1 - traits.resp_f) - - return data - - -@pytest.mark.parametrize(argnames="row", argvalues=np.arange(0, 100, 10)) -def test_tmodel_init(rvalues, row): - """Test the geometry calculation of __init__ using the exemplar R data.""" - - from pyrealm.constants.tmodel_const import TModelTraits - from pyrealm.tmodel import TTree - - row = rvalues.iloc[row] - ttree = TTree(traits=TModelTraits, diameters=row["diameter"]) - - for geom_est in ("height", "crown_area", "mass_fol", "mass_stm", "mass_swd"): - assert np.allclose(getattr(ttree, geom_est), row[geom_est]) - - -@pytest.mark.parametrize(argnames="row", argvalues=np.arange(0, 100, 10)) -def test_tmodel_reset_diameters(rvalues, row): - """Test the geometry calculation using reset_diameters using the exemplar R data.""" - - from pyrealm.constants.tmodel_const import TModelTraits - from pyrealm.tmodel import TTree - - ttree = TTree(diameters=0.001, traits=TModelTraits()) - row = rvalues.iloc[row] - ttree.reset_diameters(row["diameter"]) - - for geom_est in ("height", "crown_area", "mass_fol", "mass_stm", "mass_swd"): - assert np.allclose(getattr(ttree, geom_est), row[geom_est]) - - -def test_tmodel_init_array(rvalues): - """Test geometry calculation using an __init__array using the exemplar R data.""" - - from pyrealm.constants.tmodel_const import TModelTraits - from pyrealm.tmodel import TTree - - diams = np.array(rvalues["diameter"]) - ttree = TTree(diameters=diams, traits=TModelTraits) - - for geom_est in ("height", "crown_area", "mass_fol", "mass_stm", "mass_swd"): - vals = rvalues[geom_est] - assert np.allclose(getattr(ttree, geom_est), vals) - - -@pytest.mark.parametrize( - argnames="sequence, raises", - argvalues=[ - ("init_only", pytest.raises(RuntimeError)), - ("init_grow", does_not_raise()), - ("init_set", pytest.raises(RuntimeError)), - ("init_grow_set", pytest.raises(RuntimeError)), - ("init_grow_set_grow", does_not_raise()), - ], -) -def test_tmodel_growth_access(sequence, raises): - """Check masking of growth estimates. - - Tests that the accessors for growth properties are properly masked when not set - or outdated. - """ - from pyrealm.constants.tmodel_const import TModelTraits - from pyrealm.tmodel import TTree - - ttree = TTree(diameters=0.1, traits=TModelTraits()) - - if sequence == "init_only": - pass - elif sequence == "init_grow": - ttree.calculate_growth(7) - elif sequence == "init_set": - ttree.reset_diameters(0.2) - elif sequence == "init_grow_set": - ttree.calculate_growth(7) - ttree.reset_diameters(0.2) - elif sequence == "init_grow_set_grow": - ttree.calculate_growth(7) - ttree.reset_diameters(0.2) - ttree.calculate_growth(7) - - with raises: - _ = ttree.delta_d - - -@pytest.mark.parametrize(argnames="row", argvalues=np.arange(0, 100, 10)) -def test_tmodel_calculate_growth(rvalues, row): - """Test calculate_growth with scalars. - - Runs a test of the tmodel.TTree against output from the R implementation. The values - in the test come from simulating a 100 year run starting from a stem diameter of - 0.1. Each row in the file is the successive growth, but this test just runs some - values from the sequence. - """ - - from pyrealm.constants.tmodel_const import TModelTraits - from pyrealm.tmodel import TTree - - # create a tree with the initial diameter given in the row - row = rvalues.iloc[row] - traits = TModelTraits() - ttree = TTree(diameters=row["diameter"], traits=traits) - ttree.calculate_growth(row["potential_gpp"] / (1 - traits.resp_f)) - - for growth_est in ( - "delta_d", - "gpp_actual", - "resp_swd", - "resp_frt", - "delta_mass_stm", - "delta_mass_frt", - ): - assert np.allclose(getattr(ttree, growth_est), row[growth_est]) - - -def test_tmodel_calculate_growth_array(rvalues): - """Test calculate_growth with an array.""" - - from pyrealm.constants.tmodel_const import TModelTraits - from pyrealm.tmodel import TTree - - # create a tree with the initial diameter given in the row - diams = np.array(rvalues["diameter"]) - traits = TModelTraits() - ttree = TTree(diameters=diams, traits=traits) - ttree.calculate_growth(rvalues["potential_gpp"] / (1 - traits.resp_f)) - - for growth_est in ( - "delta_d", - "gpp_actual", - "resp_swd", - "resp_frt", - "delta_mass_stm", - "delta_mass_frt", - ): - assert np.allclose(getattr(ttree, growth_est), rvalues[growth_est]) - - -# @pytest.mark.parametrize( -# "varname", -# [ -# (0, ("dD", "delta_d")), -# (1, ("D", "diameter")), -# (2, ("H", "height")), -# (3, ("Ac", "crown_area")), -# (4, ("Wf", "mass_fol")), -# (5, ("Ws", "mass_stm")), -# (6, ("Wss", "mass_swd")), -# (7, ("GPP", "gpp_actual")), -# (8, ("Rm1", "resp_swd")), -# (9, ("Rm2", "resp_frt")), -# (10, ("dWs", "delta_mass_stm")), -# (11, ("dWfr", "delta_mass_frt")), -# ], -# ) -# def test_grow_ttree(rvalues, pyvalues, varname): -# """Runs a test of the tmodel.grow_ttree against the same R output. -# In this case, the iteration process through time is tested. If the -# previous test is successful then the values being fed forward in the -# iteration _should_ be all fine, but this checks that the time iteration -# is being run correctly in addition to the implementation of the model being -# correct - -# Args: -# values: Expected outputs from R implementation - -# """ - -# idx, (rv, pyv) = varname -# traits, pyvalues = pyvalues - -# # Get all the R values across timesteps into an array -# r_var = np.array([rw[rv] for rw in rvalues]) -# # Get the matching py output -# py_var = pyvalues[..., idx] - -# # Some implementation differences -# if pyv == "delta_d": -# # The R tmodel implementation rescales reported delta_d as -# # a radial increase in millimetres. -# py_var *= 500 -# elif pyv == "gpp_actual": -# # The R tmodel implementation slices off foliar respiration costs from -# # GPP before doing anything - the pyrealm.tmodel implementation keeps -# # this cost within the tree calculation -# py_var *= 1 - traits.resp_f - -# assert np.allclose(r_var, py_var) diff --git a/tests/unit/demography/test_community.py b/tests/unit/demography/test_community.py index bc6620a7..b4a24635 100644 --- a/tests/unit/demography/test_community.py +++ b/tests/unit/demography/test_community.py @@ -15,7 +15,7 @@ def fixture_expected(fixture_flora): Community instance maintains row order and calculation as expected. """ - from pyrealm.demography.t_model_functions import calculate_heights + from pyrealm.demography.tmodel import calculate_heights dbh = np.array([0.2, 0.5]) diff --git a/tests/unit/demography/test_flora.py b/tests/unit/demography/test_flora.py index e0ea14b1..6bc95b3d 100644 --- a/tests/unit/demography/test_flora.py +++ b/tests/unit/demography/test_flora.py @@ -422,7 +422,7 @@ def test_StemTraits(fixture_flora): def test_StemTraits_CohortMethods(fixture_flora): """Test the StemTraits inherited cohort methods.""" - from pyrealm.demography.t_model_functions import StemTraits + from pyrealm.demography.tmodel import StemTraits # Construct some input data with duplicate PFTs by doubling the fixture_flora data flora_df = fixture_flora.to_pandas() diff --git a/tests/unit/demography/test_t_model_functions.py b/tests/unit/demography/test_t_model_functions.py index 535165c0..44bc49d4 100644 --- a/tests/unit/demography/test_t_model_functions.py +++ b/tests/unit/demography/test_t_model_functions.py @@ -1,4 +1,4 @@ -"""test the functions in t_model_functions.py.""" +"""test the functions in tmodel.py.""" from contextlib import nullcontext as does_not_raise from unittest.mock import patch @@ -17,7 +17,7 @@ def test_calculate_crown_r_0_values(crown_areas, expected_r0): """Test happy path for calculating r_0.""" - from pyrealm.demography.t_model_functions import calculate_crown_r0 + from pyrealm.demography.tmodel import calculate_crown_r0 q_m = np.array([2.9038988210485766, 2.3953681843215673]) actual_r0_values = calculate_crown_r0(q_m=q_m, crown_area=crown_areas) @@ -211,7 +211,7 @@ def test_calculate_heights( exp_shape, ): """Tests calculation of heights of tree from diameter.""" - from pyrealm.demography.t_model_functions import calculate_heights + from pyrealm.demography.tmodel import calculate_heights with outcome as excep: result = calculate_heights( @@ -239,7 +239,7 @@ def test_calculate_dbh_from_height( ): """Tests inverted calculation of dbh from height.""" - from pyrealm.demography.t_model_functions import calculate_dbh_from_height + from pyrealm.demography.tmodel import calculate_dbh_from_height with outcome as excep: result = calculate_dbh_from_height( @@ -267,7 +267,7 @@ def test_calculate_crown_areas( ): """Tests calculation of crown areas of trees.""" - from pyrealm.demography.t_model_functions import calculate_crown_areas + from pyrealm.demography.tmodel import calculate_crown_areas with outcome as excep: result = calculate_crown_areas( @@ -296,7 +296,7 @@ def test_calculate_crown_fractions( ): """Tests calculation of crown fraction of trees.""" - from pyrealm.demography.t_model_functions import calculate_crown_fractions + from pyrealm.demography.tmodel import calculate_crown_fractions with outcome as excep: result = calculate_crown_fractions( @@ -324,7 +324,7 @@ def test_calculate_stem_masses( ): """Tests calculation of stem masses of trees.""" - from pyrealm.demography.t_model_functions import calculate_stem_masses + from pyrealm.demography.tmodel import calculate_stem_masses with outcome as excep: result = calculate_stem_masses( @@ -352,7 +352,7 @@ def test_calculate_foliage_masses( ): """Tests calculation of stem masses of trees.""" - from pyrealm.demography.t_model_functions import calculate_foliage_masses + from pyrealm.demography.tmodel import calculate_foliage_masses with outcome as excep: result = calculate_foliage_masses( @@ -380,7 +380,7 @@ def test_calculate_sapwood_masses( ): """Tests calculation of stem masses of trees.""" - from pyrealm.demography.t_model_functions import calculate_sapwood_masses + from pyrealm.demography.tmodel import calculate_sapwood_masses with outcome as excep: result = calculate_sapwood_masses( @@ -410,7 +410,7 @@ def test_calculate_whole_crown_gpp( ): """Tests calculation of whole crown GPP.""" - from pyrealm.demography.t_model_functions import calculate_whole_crown_gpp + from pyrealm.demography.tmodel import calculate_whole_crown_gpp with outcome as excep: result = calculate_whole_crown_gpp( @@ -439,7 +439,7 @@ def test_calculate_sapwood_respiration( ): """Tests calculation of sapwood respiration.""" - from pyrealm.demography.t_model_functions import calculate_sapwood_respiration + from pyrealm.demography.tmodel import calculate_sapwood_respiration with outcome as excep: result = calculate_sapwood_respiration( @@ -472,7 +472,7 @@ def test_calculate_foliar_respiration( broadcasting """ - from pyrealm.demography.t_model_functions import calculate_foliar_respiration + from pyrealm.demography.tmodel import calculate_foliar_respiration with outcome as excep: result = calculate_foliar_respiration( @@ -503,7 +503,7 @@ def test_calculate_fine_root_respiration( ): """Tests calculation of fine root respiration.""" - from pyrealm.demography.t_model_functions import calculate_fine_root_respiration + from pyrealm.demography.tmodel import calculate_fine_root_respiration with outcome as excep: result = calculate_fine_root_respiration( @@ -532,7 +532,7 @@ def test_calculate_net_primary_productivity( ): """Tests calculation of net primary productivity.""" - from pyrealm.demography.t_model_functions import ( + from pyrealm.demography.tmodel import ( calculate_net_primary_productivity, ) @@ -566,7 +566,7 @@ def test_calculate_foliage_and_fine_root_turnover( ): """Tests calculation of foliage and fine root turnover.""" - from pyrealm.demography.t_model_functions import ( + from pyrealm.demography.tmodel import ( calculate_foliage_and_fine_root_turnover, ) @@ -598,7 +598,7 @@ def test_calculate_growth_increments( ): """Tests calculation of growth increments.""" - from pyrealm.demography.t_model_functions import ( + from pyrealm.demography.tmodel import ( calculate_growth_increments, ) @@ -639,7 +639,7 @@ def test_calculate_dbh_from_height_edge_cases(): * If H = h_max, dbh is infinite. """ - from pyrealm.demography.t_model_functions import calculate_dbh_from_height + from pyrealm.demography.tmodel import calculate_dbh_from_height pft_h_max_values = np.array([20, 30]) pft_a_hd_values = np.array([116.0, 116.0]) @@ -693,12 +693,12 @@ def test_StemAllometry_validation(rtmodel_flora, at_dbh, outcome, excep_msg): """ from pyrealm.demography.core import _validate_demography_array_arguments - from pyrealm.demography.t_model_functions import StemAllometry + from pyrealm.demography.tmodel import StemAllometry with ( outcome as excep, patch( - "pyrealm.demography.t_model_functions._validate_demography_array_arguments", + "pyrealm.demography.tmodel._validate_demography_array_arguments", wraps=_validate_demography_array_arguments, ) as val_func_patch, ): @@ -712,7 +712,7 @@ def test_StemAllometry_validation(rtmodel_flora, at_dbh, outcome, excep_msg): def test_StemAllometry(rtmodel_flora, rtmodel_data): """Test the StemAllometry class and inherited methods.""" - from pyrealm.demography.t_model_functions import StemAllometry + from pyrealm.demography.tmodel import StemAllometry stem_allometry = StemAllometry( stem_traits=rtmodel_flora, at_dbh=rtmodel_data["dbh"][:, [0]] @@ -740,7 +740,7 @@ def test_StemAllometry(rtmodel_flora, rtmodel_data): def test_StemAllometry_CohortMethods(rtmodel_flora, rtmodel_data): """Test the StemAllometry inherited cohort methods.""" - from pyrealm.demography.t_model_functions import StemAllometry + from pyrealm.demography.tmodel import StemAllometry stem_allometry = StemAllometry( stem_traits=rtmodel_flora, at_dbh=rtmodel_data["dbh"][:, [0]] @@ -772,7 +772,7 @@ def test_StemAllometry_CohortMethods(rtmodel_flora, rtmodel_data): def test_StemAllocation(rtmodel_flora, rtmodel_data): """Test the StemAllometry class.""" - from pyrealm.demography.t_model_functions import StemAllocation, StemAllometry + from pyrealm.demography.tmodel import StemAllocation, StemAllometry stem_allometry = StemAllometry( stem_traits=rtmodel_flora, at_dbh=rtmodel_data["dbh"][:, [0]] @@ -840,7 +840,7 @@ def test_StemAllocation_validation(rtmodel_flora, pot_gpp, outcome, excep_msg): """ from pyrealm.demography.core import _validate_demography_array_arguments - from pyrealm.demography.t_model_functions import StemAllocation, StemAllometry + from pyrealm.demography.tmodel import StemAllocation, StemAllometry # Calculate a constant allometry allom = StemAllometry(stem_traits=rtmodel_flora, at_dbh=np.ones((4, 3))) @@ -848,7 +848,7 @@ def test_StemAllocation_validation(rtmodel_flora, pot_gpp, outcome, excep_msg): with ( outcome as excep, patch( - "pyrealm.demography.t_model_functions._validate_demography_array_arguments", + "pyrealm.demography.tmodel._validate_demography_array_arguments", wraps=_validate_demography_array_arguments, ) as val_func_patch, ):