Annual fAPAR_max and LAI_max

pyrealm/phenology/fapar_limitation.py

@@ -0,0 +1,208 @@
+"""Class to compute the fAPAR_max and annual peak Leaf Area Index (LAI)."""
+
+import numpy as np
+from numpy.ma.core import zeros
+from numpy.typing import NDArray
+from typing_extensions import Self
+
+from pyrealm.pmodel import PModel, SubdailyScaler
+
+
+def get_annual(
+    x: NDArray,
+    datetimes: NDArray[np.datetime64],
+    growing_season: NDArray[np.bool],
+    method: str,
+) -> NDArray:
+    """Computes an array of the annual total or mean of an entity x given datetimes.
+
+    Args:
+        x: Array of values to be converted to annual values. Should be either daily (
+        same datetimes as growing_season) or subdaily (same datetimes as datetimes
+        array)
+        datetimes: Datetimes of the measurements as np.datetime62 arrays.
+        growing_season: Bool array of days, indicating whether they are ain growing
+        season or not.
+        method: Either "total" (sum all values of the year) or "mean" (take the mean
+        of all values of the year)
+    """
+
+    # Extract years from datetimes
+    all_years = [np.datetime64(i, "Y") for i in datetimes]
+
+    # Create scaler object to handle conversion between scales
+    scaler = SubdailyScaler(datetimes)
+    scaler.set_nearest(np.timedelta64(12, "h"))
+
+    if len(x) == len(growing_season):
+        daily_x = x
+    elif len(x) == len(datetimes):
+        # Convert values to daily to match with growing_season
+        daily_x = scaler.get_daily_means(x)
+    else:
+        raise ValueError("Input array does not fit datetimes nor growing_season array")
+
+    # Get rid of that extra dimension of size 1
+    years_by_day = np.squeeze(scaler.get_window_values(np.asarray(all_years)))
+    # Which years are present?
+    years = np.unique(all_years)
+
+    # Compute annual totals or means, only taking into account the days which are in
+    # growing season.
+
+    annual_x = zeros(len(years))
+
+    if method == "total":
+        for i in range(len(years)):
+            annual_x[i] = sum(
+                daily_x[growing_season & (years_by_day == years[i])].astype(np.int64)
+            )
+    elif method == "mean":
+        for i in range(len(years)):
+            annual_x[i] = np.mean(
+                daily_x[growing_season & (years_by_day == years[i])].astype(np.int64)
+            )
+    else:
+        raise ValueError("No valid method given for annual values")
+
+    return annual_x
+
+
+def compute_annual_total_precip(
+    precip: NDArray, datetimes: NDArray[np.datetime64], growing_season: NDArray[np.bool]
+) -> NDArray:
+    """Returns the sum of annual precipitation."""
+
+    annual_precip = get_annual(precip, datetimes, growing_season, "total")
+
+    return convert_precipitation_to_molar(annual_precip)
+
+
+def convert_precipitation_to_molar(precip_mm: NDArray) -> NDArray:
+    """Convert precipitation from mm/m2 to mol/m2.
+
+    - 1 mm/m2 = 1000000 mm3 = 1000 mL = 1 L
+    - Molar mass of water = 18g
+    - Assuming density = 1 (but really temp varying), molar volume = 18mL
+    - So 1 mm/m2 = 1000 / 18 ~ 55.55 mols/m2
+    """
+
+    water_mm_to_mol = 1000 / 18
+    precip_molar = precip_mm * water_mm_to_mol
+
+    return precip_molar
+
+
+class FaparLimitation:
+    r"""FaparLimitation class to compute fAPAR_max and LAI.
+
+    This class takes the annual total potential GPP and precipitation, the annual
+    mean CA, Chi and VPD, as well as the aridity index and some constants to compute
+    the annual peak fractional absorbed photosynthetically active radiation (
+    fAPAR_max) and annual peak Leaf Area Index (LAI).
+    """
+
+    def __init__(
+        self,
+        annual_total_potential_gpp: NDArray[np.float64],
+        annual_mean_ca: NDArray[np.float64],
+        annual_mean_chi: NDArray[np.float64],
+        annual_mean_vpd: NDArray[np.float64],
+        annual_total_precip: NDArray[np.float64],
+        aridity_index: NDArray[np.float64],
+        z: float = 12.227,
+        k: float = 0.5,
+    ) -> None:
+        r"""Annual peak fractional absorbed photosynthetically active radiation (fAPAR).
+
+        Computes fAPAR_max as minimum of a water-limited and an energy-limited
+        quantity.
+
+        Args:
+            annual_total_potential_gpp: Aka A_0, the annual sum of potential GPP.
+              Potential GPP would be achieved if fAPAR = 1. [mol m^{-2} year^{-1}]
+            annual_mean_ca: Ambient CO2 partial pressure. [Pa]
+            annual_mean_chi: Annual mean ratio of leaf-internal CO2 partial pressure to
+               c_a during the growing season (>0℃). [Pa]
+            annual_mean_vpd: Aka D, annual mean vapour pressure deficit (VPD) during the
+              growing season (>0℃). [Pa]
+            annual_total_precip: Aka P, annual total precipitation. [mol m^{-2}
+             year^{-1}]
+            aridity_index: Aka AI, climatological estimate of local aridity index.
+            z: accounts for the costs of building and maintaining leaves and the total
+              below-ground allocation required to support the nutrient demand of those
+              leaves. [mol m^{-2} year^{-1}]
+            k: Light extinction coefficient.
+        """
+
+        #  f_0 is the ratio of annual total transpiration of annual total
+        #  precipitation, which is an empirical function of the climatic Aridity Index
+        #  (AI).
+        b = 0.604169
+        f_0 = 0.65 * np.exp(-b * np.log(aridity_index / 1.9))
+
+        fapar_energylim = 1.0 - z / (k * annual_total_potential_gpp)
+        fapar_waterlim = (
+            f_0
+            * annual_total_precip
+            * annual_mean_ca
+            * (1 - annual_mean_chi)
+            / (1.6 * annual_mean_vpd * annual_total_potential_gpp)
+        )
+
+        self.faparmax = -9999 * np.ones(np.shape(fapar_waterlim))
+        self.energylim = -9999 * np.ones(np.shape(fapar_waterlim))
+        self.annual_precip_molar = annual_total_precip
+
+        for i in range(len(fapar_waterlim)):
+            if fapar_waterlim[i] < fapar_energylim[i]:
+                self.faparmax[i] = fapar_waterlim[i]
+                self.energylim[i] = False
+            else:
+                self.faparmax[i] = fapar_energylim[i]
+                self.energylim[i] = True
+
+        self.laimax = -(1 / k) * np.log(1.0 - self.faparmax)
+
+    @classmethod
+    def from_pmodel(
+        cls,
+        pmodel: PModel,
+        growing_season: NDArray[np.bool],
+        datetimes: NDArray[np.datetime64],
+        precip: NDArray[np.float64],
+        aridity_index: NDArray[np.float64],
+    ) -> Self:
+        r"""Get FaparLimitation from PModel input.
+
+        Computes the input for fAPAR_max from the P Model and additional inputs.
+
+        Args:
+            pmodel: pyrealm.pmodel.PModel
+            growing_season: Bool array indicating which times are within growing
+              season by some definition and implementation.
+            datetimes: Array of datetimes to consider.
+            precip: Precipitation for given datetimes.
+            aridity_index: Climatological estimate of local aridity index.
+        """
+
+        annual_total_potential_gpp = get_annual(
+            pmodel.gpp, datetimes, growing_season, "total"
+        )
+        annual_mean_ca = get_annual(pmodel.env.ca, datetimes, growing_season, "mean")
+        annual_mean_chi = get_annual(
+            pmodel.optchi.chi.round(5), datetimes, growing_season, "mean"
+        )
+        annual_mean_vpd = get_annual(pmodel.env.vpd, datetimes, growing_season, "mean")
+        annual_total_precip = compute_annual_total_precip(
+            precip, datetimes, growing_season
+        )
+
+        return cls(
+            annual_total_potential_gpp,
+            annual_mean_ca,
+            annual_mean_chi,
+            annual_mean_vpd,
+            annual_total_precip,
+            aridity_index,
+        )

pyrealm_build_data/LAI_in_pyrealm/.gitignore

@@ -0,0 +1,2 @@
+outputs/*
+reports/*

pyrealm_build_data/LAI_in_pyrealm/DE-GRI_site_data.json

@@ -0,0 +1,10 @@
+{
+  "Dataset": "Flux2015",
+  "Site": "DE-Gri",
+  "Lat": 50.95,
+  "Lon": 13.5126,
+  "Vegetation": "GRA",
+  "Elev": 385,
+  "AI": 1.094517,
+  "f0": 0.540873739473232
+}