The canopy model

Run this notebook

Warning

This area of pyrealm is in active development and this notebook currently contains notes and initial demonstration code.

Hide code cell source

 
import warnings

from matplotlib import pyplot as plt
import matplotlib as mpl
from matplotlib.patches import Polygon
import numpy as np
import pandas as pd

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.tmodel import StemAllometry
from pyrealm.core.experimental import ExperimentalFeatureWarning

warnings.filterwarnings(
    "ignore",
    category=ExperimentalFeatureWarning,
)

np.set_printoptions(precision=2)

The canopy module in pyrealm is used to calculate a vertically structured model of leaf distribution for a plant community. The purpose of the canopy module is two-fold:

  1. to calculate the vertical distribution of crown and leaf area, including the partitioning of canopy into discrete vertical layers, and

  2. to estimate light capture through canopy layers.

This page describes the model of the vertical canopy structure and the light capture model is presented separately.

Vertically structured canopies

The simplest “big leaf” model of a canopy, represents the entire crown crown (\(A_c\)) of an individual tree as a single flat disc at the top of the stem. However, in pyrealm the canopy model uses the crown shape described in the crown model to model how leaf area accumulates vertically through a community.

Community definition

To recap, a plant community consists of a number of individual stems that are growing together in a location. The community structure groups individuals together into cohorts that are defined as:

However, for the purposes of the canopy model, it is much simpler to consider the community as a collection of individual stems, some of which just happen to share identical properties.

Community projected crown area

Each individual stem in a community has an projected crown area that describes how the crown area of that stem accumulates with height (\(z\)) from the top of the tree to the ground. The form of that curve varies with the size and PFT of the stem. If we have a community with \(N_s\) individuals, then each stem \(i = 1, ..., N_s\) has a projected crown area function \(A_{p}(z)_i\).

If we take the sum across all individuals of this function at a height \(z\), we get the total community projected crown area across all individual stems at that same height:

\[ C_p(z) = \sum_{j=1}^{N_s}{A_{p}(z)_i} \]

To demonstrate this, the code below creates a plant community that contains 12 individual stems, grouped into 3 cohorts with different stem sizes and crown shapes.

# Define PFTs
short_pft = PlantFunctionalType(name="short", h_max=15, m=2, n=4, f_g=0, ca_ratio=380)
tall_pft = PlantFunctionalType(
    name="tall", h_max=30, m=2, n=3, par_ext=0.6, f_g=0, ca_ratio=500
)

# Create the flora
flora = Flora([short_pft, tall_pft])

# Define a simply community with three cohorts
community = Community(
    flora=flora,
    cell_area=32,
    cell_id=1,
    cohorts=Cohorts(
        dbh_values=np.array([0.1, 0.20, 0.5]),
        n_individuals=np.array([7, 3, 2]),
        pft_names=np.array(["short", "short", "tall"]),
    ),
)

The total crown area of the individuals in each of the three cohorts is shown below:

community.stem_allometry.crown_area

array([[ 2.08,  6.07, 43.43]])

The total crown area across the community is then the sum of those crown areas across all individuals:

total_crown_area = np.sum(
    community.stem_allometry.crown_area * community.cohorts.n_individuals, keepdims=True
)
total_crown_area

array([[119.64]])

So, the total community crown area is around 120 m2. The plot below shows how \(C_p(z)\) changes with height for this community from zero at the top of the canopy to around 120 m2 at ground level. The horizontal dashed lines show the stem heights of the individuals, where the crown area of each individual starts contributing to the community wide projected crown area.

Hide code cell source

 
# Calculate the crown profiles of the individuals across vertical heights
hghts = np.linspace(community.stem_allometry.stem_height.max(), 0, num=101)[:, None]
crown_profiles = CrownProfile(community.stem_traits, community.stem_allometry, z=hghts)

# Calculate the cumulative crown area across individuals at each height
crown_area = np.nansum(
    crown_profiles.projected_crown_area * community.cohorts.n_individuals,
    axis=1,
)

# Extract the crown profiles as XY arrays for plotting
profiles = get_crown_xy(
    crown_profile=crown_profiles,
    stem_allometry=community.stem_allometry,
    attr="crown_radius",
    as_xy=True,
)


# Helper function to add vertical or horizontal lines to plots
def add_hvlines(at, axes, vertical=True):

    # Line styling
    style = dict(color="black", linewidth=0.5, linestyle="--")

    # Add lines
    if vertical:
        axes.vlines(at, ymin=0, ymax=1, transform=axes.get_xaxis_transform(), **style)
    else:
        axes.hlines(at, xmin=0, xmax=1, transform=axes.get_yaxis_transform(), **style)


def add_second_axis(at, axes, fmt, xaxis=True):

    # Get labels
    labels = [fmt.format(val=val, idx=idx + 1) for idx, val in enumerate(at)]

    if xaxis:
        secax = axes.secondary_xaxis("top")
        secax.set_xticks(at, labels=labels)
    else:
        secax = axes.secondary_yaxis("right")
        secax.set_yticks(at, labels=labels)


# Create a side by side plot of the canopies and the cumulative community crown area
fig, (ax1, ax2) = plt.subplots(
    ncols=2, sharey=True, width_ratios=[1, 2], figsize=(8, 4)
)

for idx, crown in enumerate(profiles):

    # Get spaced but slightly randomized stem locations
    n_stems = community.cohorts.n_individuals[idx]
    stem_locations = np.linspace(0, 10, num=n_stems) + np.random.normal(size=n_stems)

    # Plot the crown model for each stem
    for stem_loc in stem_locations:
        ax1.add_patch(Polygon(crown + np.array([stem_loc, 0]), color="#00550055"))

ax1.autoscale_view()
add_hvlines(community.stem_allometry.stem_height, ax1, False)
ax1.set_ylabel("Height (m)")

ax2.plot(crown_area, hghts)
add_hvlines(community.stem_allometry.stem_height, ax2, False)
add_hvlines([0, total_crown_area], ax2, True)
_ = ax2.set_xlabel("Community projected crown area ($C_{p}(z), m^2$)")

plt.tight_layout()
../../_images/a8f6877b5d7145122ccfb1cff869d65dc0eadedde3c0b44e8829c8dd5ac989d8.png

If the community is growing in an area greater than 120 m2, then the individuals can avoid overlapping any of their crown area. However, this is not possible if the community is growing in a smaller area, and some crown area must overlie other parts of the community, creating a vertical structure.

The canopy module implements the perfect plasticity approximation (PPA) model (Purves et al., 2008) to generate this structure. The PPA model assumes that all the individuals within the community are able to plastically arrange their crown at each height within the broader canopy of the community to fill available space. When the available horizontal space (\(A\)) is filled by the cumulative crown area across all individuals at a height \(z\), that canopy at that height forms a layer that is considered “full” (or “closed”). The remaining crown area below this layer then begins to accumulate into a new layer. This process repeats downward until the entire canopy structure is formed down to the ground.

To fit this model, we need to find the heights \(z^*_l\) at which the projected community crown area occupies multiples of the available area \(A\). This gives the heights at which each successive canopy layer closes for canopy layers \(l = 1, ..., l_m\). This is given as the equality:

\[ C_p(z^*_l) = l A \]

An extra detail is that the canopy model includes a community-level canopy gap fraction (\(f_G\)) that captures the overall proportion of the canopy area that is left unfilled by canopy. This gap fraction, capturing processes such as crown shyness, describes the proportion of open sky visible from the forest floor. This term reduces the amount of \(A\) that can be occupied by crown area before a layer closes:

\[ C_p(z^*_l) = l A (1- f_G) \]

The values of \(z^*_l\) cannot be found analytically, so we use numerical methods to find values of \(z^*_l\) for \(l = 1,..., l_m\) that satisfy:

\[ C_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 \]

The community above has an available space \(A=32\;m^2\). If we fit the PPA model to this community we get the following heights:

# Fit the canopy model
canopy_ppa = Canopy(community=community, fit_ppa=True)
canopy_ppa.heights

array([[23.45],
       [21.79],
       [10.91],
       [ 0.  ]])

We plot those layer heights below. The vertical dashed lines in the right hand plot show the values \(lA\) at which the cumulative space across layers forms closed layers: that is \(32 \times 1 = 32\), \(32 \times 2 = 64\), etc. . The horizontal dashed lines then show those calculated vertical heights. These are the canopy heights at which the vertical lines intersect the cumulative community projected crown area: where the total crown area across the individuals in the community fills each available layer.

Hide code cell source

 
# Create a plot of the canopy structure
fig, (ax1, ax2) = plt.subplots(
    ncols=2, sharey=True, width_ratios=[1, 2], figsize=(8, 4)
)

# Create a canopy plot
for idx, crown in enumerate(profiles):

    # Get spaced but slightly randomized stem locations
    n_stems = community.cohorts.n_individuals[idx]
    stem_locations = np.linspace(0, 10, num=n_stems) + np.random.normal(size=n_stems)

    # Plot the crown model for each stem
    for stem_loc in stem_locations:
        ax1.add_patch(Polygon(crown + np.array([stem_loc, 0]), color="#00550055"))

add_hvlines(canopy_ppa.heights, ax1, False)
ax1.set_ylabel("Height (m)")

# Plot community projected crown area
ax2.plot(crown_area, hghts)
ax2.set_xlabel("Cumulative crown area (m2)")

# Add lines to show PPA closure areas and heights
area_values = np.arange(1, 5) * 32
add_hvlines(canopy_ppa.heights, ax2, False)
add_hvlines(area_values, ax2, True)

# Add secondary axes to give values
add_second_axis(canopy_ppa.heights.squeeze(), ax2, "$z^*_{idx} = ${val:0.2f}", False)
add_second_axis(area_values, ax2, "$A_{idx} = ${val}", True)

plt.tight_layout()
../../_images/8fc13c207796d636c1bf218f996f213034c0c7ee6687b8585e89bec7e1e72067.png

We can look in detail at the how much crown area is in each layer from each individual. The canopy model object contains a crown profile which records the projected crown area at each height. Note how the bottom row - the total projected (i.e. cumulative) crown for the individuals in each cohort - matches the values above from the stem allometry calculations for the community.

canopy_ppa.crown_profile.projected_crown_area

array([[ 0.  ,  0.  , 16.  ],
       [ 0.  ,  0.  , 32.  ],
       [ 0.  ,  3.04, 43.43],
       [ 2.08,  6.07, 43.43]])

Those values are the projected crown areas of each individual within the cohorts: the accumulating crown area from the top of the canopy down to the ground. We can take the differences between vertical layers to give the actual amount of crown area in each layer for each individual within a cohort.

individual_crown_in_layer = np.diff(
    canopy_ppa.crown_profile.projected_crown_area, prepend=0, axis=0
)
individual_crown_in_layer

array([[ 0.  ,  0.  , 16.  ],
       [ 0.  ,  0.  , 16.  ],
       [ 0.  ,  3.04, 11.43],
       [ 2.08,  3.03,  0.  ]])

To bring those values back to the community model: since each value represents a cohort, we can multiply those values by the number of individuals and then find the row sums to show the total amount of crown area in each layer across all individuals and cohorts. As expected, the crown area in each layer matches the 32 m2 of available space, except for the last layer that is not completely filled.

np.sum(individual_crown_in_layer * community.cohorts.n_individuals, axis=1)

array([32.  , 32.  , 32.  , 23.64])

Crown and canopy gap fractions

The canopy and crown models are extended by providing two gap fractions:

  • The crown gap fraction (\(f_g\)) is described in detail in the crown model and captures how an individual tree crown may contain holes that displace leaf area further down into the canopy. Crown gaps do not push all the way down to the ground - they only allow light to penetrate more deeply into the crown.

  • The canopy gap fraction (\(f_G\)) is described above and captures how the canopy across the whole community may leave space unfilled in the canopy. Light passing through canopy gaps passes unimpeded from the top of the canopy down to the ground.

The code below alters the community and canopy model used above to include both crown and canopy gap fractions.

# Define gappy PFTs
gappy_short_pft = PlantFunctionalType(
    name="short",
    h_max=15,
    m=1.5,
    n=1.5,
    f_g=0.1,
    ca_ratio=380,
)
gappy_tall_pft = PlantFunctionalType(
    name="tall", h_max=30, m=1.5, n=2, par_ext=0.6, f_g=0.1, ca_ratio=500
)

# Create the flora
gappy_flora = Flora([gappy_short_pft, gappy_tall_pft])

# Define community with three cohorts
gappy_community = Community(
    flora=gappy_flora,
    cell_area=32,
    cell_id=1,
    cohorts=Cohorts(
        dbh_values=np.array([0.1, 0.20, 0.5]),
        n_individuals=np.array([7, 3, 2]),
        pft_names=np.array(["short", "short", "tall"]),
    ),
)

# Calculate the canopy profile across vertical heights
gappy_canopy_ppa = Canopy(
    community=gappy_community, fit_ppa=True, canopy_gap_fraction=1 / 8
)

/home/docs/checkouts/readthedocs.org/user_builds/pyrealm/checkouts/latest/pyrealm/demography/crown.py:72: RuntimeWarning: invalid value encountered in power
  q_z = m * n * z_over_height ** (n - 1) * (1 - z_over_height**n) ** (m - 1)

Hide code cell source

 
# Calculate the crown profiles for the gappy community
gappy_crown_profiles = CrownProfile(
    gappy_community.stem_traits, gappy_community.stem_allometry, z=hghts
)

# Calculate the cumulative crown area across individuals for the gappy community
gappy_crown_area = np.nansum(
    gappy_crown_profiles.projected_crown_area * gappy_community.cohorts.n_individuals,
    axis=1,
)

# Calculate leaf areas for each community
leaf_area = np.nansum(
    crown_profiles.projected_leaf_area * community.cohorts.n_individuals,
    axis=1,
)

gappy_leaf_area = np.nansum(
    gappy_crown_profiles.projected_leaf_area * gappy_community.cohorts.n_individuals,
    axis=1,
)

The plots below show the cumulative projected leaf and crown areas across the individuals in the two communities along with the PPA layer closure heights. The first thing to note is that the projected crown area for the communities is identical: the cohorts have the same crown shapes, stem heights and number of individuals. For the non-gappy community, the projected leaf area and projected crown area are also identical. However:

  • The gappy community has a canopy gap fraction of 1/8, reducing the available total crown area within a single layer to 28 m2. The layer closure heights under the PPA model are displaced upwards and an extra closed canopy layer is needed to fit the total community crown area of ~120 m2 into the available area.

  • The crown gap fractions in the gappy community are not zero and so the projected leaf area is displaced downwards within the canopy. Note that this does not affect the location of the canopy layer closure heights, but it does affect the location of leaf area for the purposes of light capture.

Hide code cell source

 
fig, (ax1, ax2) = plt.subplots(ncols=2, sharey=True, figsize=(10, 5))

# Plot projected community crown and leaf area for non gappy community
ax1.plot(crown_area, hghts)
ax1.plot(leaf_area, hghts, linestyle="--", color="red")

# Add PPA closure height and area lines and axes
add_hvlines(canopy_ppa.heights, ax1, False)
add_hvlines(area_values, ax1, True)
add_second_axis(canopy_ppa.heights.squeeze(), ax1, "$z^*_{idx} = ${val:0.2f}", False)
add_second_axis(area_values, ax1, "$A_{idx} = ${val}", True)

text_args = dict(x=0.95, y=0.95, ha="right", va="top", backgroundcolor="white")
ax1.text(s="Non-gappy community", transform=ax1.transAxes, **text_args)
ax1.set_xlabel("Projected community area (m2)")
ax1.set_ylabel("Height (m)")

# Gappy community
gappy_area_values = np.arange(1, 5) * 28

ax2.plot(gappy_crown_area, hghts, label="Crown area")
ax2.plot(gappy_leaf_area, hghts, linestyle="--", color="red", label="Leaf area")

# Add PPA closure height and area lines and axes
add_hvlines(gappy_canopy_ppa.heights.squeeze(), ax2, False)
add_hvlines(gappy_area_values, ax2, True)
add_second_axis(
    gappy_canopy_ppa.heights.squeeze(), ax2, "$z^*_{idx} = ${val:0.2f}", False
)
add_second_axis(gappy_area_values, ax2, "$A_{idx} = ${val}", True)

ax2.text(s="Gappy community", transform=ax2.transAxes, **text_args)
ax2.set_xlabel("Projected community area (m2)")
ax2.legend(framealpha=1.0)

plt.tight_layout()
../../_images/90d6d31505391077aa4586e5925c2224a3b2a005c3976d2f89e3181f448a75d3.png