diff --git a/docs/source/users/tmodel/canopy.md b/docs/source/users/demography/canopy.md similarity index 99% rename from docs/source/users/tmodel/canopy.md rename to docs/source/users/demography/canopy.md index b0911b87..a285909b 100644 --- a/docs/source/users/tmodel/canopy.md +++ b/docs/source/users/demography/canopy.md @@ -6,13 +6,20 @@ jupytext: format_name: myst format_version: 0.13 kernelspec: - display_name: python3 + display_name: Python 3 (ipykernel) language: python name: python3 --- # Canopy model +:::{admonition} Warning + +This area of `pyrealm` is in active development and this notebook currently contains +notes and initial demonstration code. + +::: + This notebook walks through the steps in generating the canopy model as (hopefully) used in Plant-FATE. diff --git a/docs/source/users/demography/crown.md b/docs/source/users/demography/crown.md index a107675e..1eba4486 100644 --- a/docs/source/users/demography/crown.md +++ b/docs/source/users/demography/crown.md @@ -43,31 +43,114 @@ from pyrealm.demography.crown import ( ``` The {mod}`pyrealm.demography` module uses three-dimensional model of crown shape to -define the vertical distribution of leaf area. This is driven by four parameters within -the {class}`~pyrealm.demography.flora.PlantFunctionalType`: +define the vertical distribution of leaf area, following the implementation of crown +shape in the Plant-FATE model {cite}`joshi:2022a`. -* The `m` and `n` parameters that define the vertical shape of the crown profile. -* The `ca_ratio` parameters that defines the size of the crown relative to the stem size. -* The `f_g` parameter that define the crown gap fraction and sets the vertical - distribution of leaves within the crown. - -The code below provides a demonstration of the impacts of each parameter and shows the -use of `pyrealm` tools to visualise the crown model for individual stems. - -## Crown shape parameters +## Crown traits -The `m` and `n` parameters define the vertical profile of the crown but are also define -two further parameters: +The crown model for a plant functional type (PFT) is driven by four traits within +the {class}`~pyrealm.demography.flora.PlantFunctionalType` class: -* `q_m`, which is used to scale the size of the crown radius to match the expected crown - area, given the crown shape. -* `z_max_prop`, which sets the height on the stem at which the maximum crown radius is found. +* 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. -The code below calculates `q_m` and `z_max_prop` for combinations of `m` and `n` and -then plots the resulting values. +### 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]) ``` @@ -90,36 +173,53 @@ ax2.set_ylabel("n") ax2.set_aspect("equal") ``` -## Plotting crown profiles +## Crown calculations in `pyrealm` -The examples below show the calculation of crown profiles for a set of plant functional -types (PFTs) with differing crown trait values. We first need to create a -{class}`~pyrealm.demography.flora.Flora` object that defines those PFTs. +The examples below show the calculation of crown variables in `pyrealm`. -```{code-cell} -flora = Flora( - [ - PlantFunctionalType(name="short", h_max=20, m=1.5, n=1.5, f_g=0, ca_ratio=20), - PlantFunctionalType( - name="medium", h_max=20, m=1.5, n=4, f_g=0.05, ca_ratio=500 - ), - PlantFunctionalType(name="tall", h_max=20, m=4, n=1.5, f_g=0.2, ca_ratio=2000), - ] -) +### 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} -pd.DataFrame({k: getattr(flora, k) for k in flora.trait_attrs}) +# 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"]] ``` -We then also need to specify the size of a stem of each PFT, along with the resulting -allometric predictions from the T Model. +The next section of code generates the `StemAllometry` to use for the profiles. +The T Model uses DBH to define stem size - here the the code is being used to +back-calculate the required DBH values to give three stems with similar heights +near the maximum height for each PFT. ```{code-cell} -# Generate the expected stem allometries at a single height for each PFT -stem_height = np.array([18, 18, 18]) +# 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 ) @@ -127,79 +227,64 @@ stem_dbh ``` ```{code-cell} +# Calculate the stem allometries allometry = StemAllometry(stem_traits=flora, at_dbh=stem_dbh) ``` -We can use {mod}`pandas` to visualise those allometric predictions. +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}) ``` -We can now use the {class}`~pyrealm.demography.crown.CrownProfile` class to -calculate the crown profile for each stem for a set of vertical heights. The heights -need to be defined as a column array, that is with a shape `(N, 1)`, in order to show -that we want the crown variables to be calculated at each height for each PFT. +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(0, 18.0, num=181)[:, None] +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 is a dictionary containing values for four crown profile -variables. +The `crown_profiles` object then provides the four crown profile attributes describe +above calculated at each height $z$: -* The relative crown radius: this value is driven purely by `m`, `n` and the stem height - and defines the vertical profile of the crown. -* The crown radius: this is simply a rescaling of the relative radius so that the - maximum crown radius matches the expected crown area predicted by the allometry of the - T Model. -* The projected crown area: this value shows how crown area accumulates from the top of - the crown to the ground. -* The projected leaf area: this value shows how _leaf_ area accumulates from the top of - the crown to the ground. The difference from the crown area is driven by the crown gap - fraction for a given PFT. +* The relative crown radius $q(z)$ +* The crown radius $r(z)$ +* The projected crown area +* The projected leaf area ```{code-cell} crown_profiles ``` -### Crown radius values - -The relative crown radius values are arbitrary - they simply define the shape of the -vertical profile. For the example PFTs in the `Flora` object, the maximum relative -radius values on each stem are: - -```{code-cell} -max_relative_crown_radius = crown_profiles.relative_crown_radius.max(axis=0) -print(max_relative_crown_radius) -``` - -However the scaled maximum radius values match the expected crown area in the allometry -table above - -```{code-cell} -max_crown_radius = crown_profiles.crown_radius.max(axis=0) -print(max_crown_radius) -print(max_crown_radius**2 * np.pi) -``` +### 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 varies with height, -* the solid line shows the actual crown radius profile with height, and -* the dotted line shows the height at which the maximum crown radius is found. +* 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)$ does define values where +$z <0$ or $z > 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 @@ -229,32 +314,51 @@ for pft_idx, offset, colour in zip((0, 1, 2), (0, 5, 12), ("r", "g", "b")): ax.set_aspect(aspect=1) ``` -### Projected crown and leaf areas +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 use the crown profile to generate projected area plots, but it is much easier to -compare the effects when comparing stems with similar crown areas. The code below -generates these new predictions for a new set of PFTs. +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} -flora2 = Flora( - [ - PlantFunctionalType(name="short", h_max=20, m=1.5, n=1.5, f_g=0, ca_ratio=380), - PlantFunctionalType( - name="medium", h_max=20, m=1.5, n=4, f_g=0.05, ca_ratio=400 - ), - PlantFunctionalType(name="tall", h_max=20, m=4, n=1.5, f_g=0.2, ca_ratio=420), - ] +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 ) -allometry2 = StemAllometry(stem_traits=flora2, at_dbh=stem_dbh) - +# 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 profile across those heights for each PFT -crown_profiles2 = CrownProfile(stem_traits=flora2, stem_allometry=allometry2, z=z) +# 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 shows how projected crown area (solid lines) and leaf area (dashed lines) -change with height along the stem. +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 @@ -271,13 +375,17 @@ fig, ax = plt.subplots(ncols=1) for pft_idx, offset, colour in zip((0, 1, 2), (0, 5, 10), ("r", "g", "b")): - ax.plot(crown_profiles2.projected_crown_area[:, pft_idx], z, color=colour) + ax.plot(area_crown_profiles.projected_crown_area[:, pft_idx], area_z, color=colour) ax.plot( - crown_profiles2.projected_leaf_area[:, pft_idx], - z, + 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)") ``` + +```{code-cell} + +``` diff --git a/docs/source/users/demography/flora.md b/docs/source/users/demography/flora.md index cb3605c8..dfe4f4e2 100644 --- a/docs/source/users/demography/flora.md +++ b/docs/source/users/demography/flora.md @@ -95,7 +95,7 @@ the PlantFATE model {cite}`joshi:2022a`. The {class}`~pyrealm.demography.flora.PlantFunctionalType` class is used define a PFT with a given name, along with the trait values associated with the PFT. By default, -values for each trait are taken from Table 1 of {cite}`Li:2014bc`, but these can be +values for each trait are taken from Table 1 of {cite:t}`Li:2014bc`, but these can be adjusted for different PFTs. The code below contains three examples that just differ in their maximum height.