Regression tests for FaparLimitation class

pyrealm_build_data/LAI_in_pyrealm/fapar_limitation_example.py

@@ -6,6 +6,7 @@
 import numpy as np
 import pandas as pd
 import xarray as xr
+from faparlim_data_io import write_faparlim_input, write_pmodel_faparlim_input
 from numpy.typing import NDArray
 from pandas import DataFrame
 from scipy.special import lambertw
@@ -242,6 +243,22 @@ def load_fluxnet_data(
     join="left",
 )
 
+write_faparlim_input(
+    ann_total_A0_subdaily_penalised, ann_total_P_fnet, aridity_index, annual_values
+)
+
+write_pmodel_faparlim_input(
+    de_gri_hh_xr["TA_F"].to_numpy(),
+    de_gri_hh_xr["VPD_F"].to_numpy(),
+    de_gri_hh_xr["CO2_F_MDS"].to_numpy(),
+    de_gri_hh_xr["PA_F"].to_numpy(),
+    de_gri_daily_values["growing_day"].to_numpy(),
+    de_gri_hh_pd.axes[0].to_numpy(),
+    de_gri_daily_values["pre"].data,
+    aridity_index.data,
+    de_gri_hh_xr["PPFD"].to_numpy(),
+)
+
 # Calculate fapar_max and LAI_max
 
 if use_pmodel:

pyrealm_build_data/LAI_in_pyrealm/faparlim_data_io.py

@@ -0,0 +1,141 @@
+"""Methods for input and output of test data for FaparLimitation class."""
+
+import numpy as np
+import xarray as xr
+from netCDF4 import Dataset
+from xarray import DataArray
+
+
+def write_faparlim_input(
+    ann_total_A0_subdaily_penalised: xr.Dataset,
+    ann_total_P_fnet: xr.Dataset,
+    aridity_index: DataArray,
+    annual_values: xr.Dataset,
+) -> None:
+    """Method for writing out the data needed for the FaparLimitation constructor."""
+
+    with Dataset("faparlim_input.nc", "w") as f:
+        f.createDimension("year", len(ann_total_A0_subdaily_penalised.year.data))
+        year = f.createVariable("year", "i", ("year",))
+        year[:] = ann_total_A0_subdaily_penalised.year.data
+
+        annual_total_A0_subdaily = f.createVariable(
+            "annual_total_A0_subdaily", float, ("year",)
+        )
+        annual_total_A0_subdaily[:] = ann_total_A0_subdaily_penalised.data
+
+        annual_total_P = f.createVariable("annual_total_P", float, ("year",))
+        annual_total_P[:] = ann_total_P_fnet.data
+
+        f.createVariable("aridity_idx", float, ())
+        f.variables["aridity_idx"][...] = aridity_index.data
+
+    annual_values.to_netcdf("faparlim_input.nc", mode="a")
+
+
+def write_pmodel_faparlim_input(
+    tc: DataArray,
+    vpd: DataArray,
+    co2: DataArray,
+    patm: DataArray,
+    growing_season: DataArray,
+    datetimes: DataArray,
+    precipitation: DataArray,
+    aridity_index: DataArray,
+    ppfd: DataArray,
+) -> None:
+    """Writing out the data needed for the from_pmodel FaparLimitation class method."""
+
+    with Dataset("faparlim_pmodel_input.nc", "w") as f:
+        f.createDimension("datatimes_length", 29)
+
+        f.createDimension("datetimes", len(datetimes))
+        datetimes_var = f.createVariable("datetimes", str, ("datetimes",))
+        datetimes_var[:] = datetimes[:].astype(str)
+
+        f.createDimension("days", len(growing_season))
+        growing_season_var = f.createVariable("growing_season", "i", ("days",))
+        growing_season_var[:] = growing_season
+
+        tc_var = f.createVariable("tc", float, ("datetimes",))
+        tc_var[:] = tc
+
+        vpd_var = f.createVariable("vpd", float, ("datetimes",))
+        vpd_var[:] = vpd
+
+        co2_var = f.createVariable("co2", float, ("datetimes",))
+        co2_var[:] = co2
+
+        patm_var = f.createVariable("patm", float, ("datetimes",))
+        patm_var[:] = patm
+
+        precip_var = f.createVariable("precipitation", float, ("days",))
+        precip_var[:] = precipitation
+
+        f.createVariable("aridity_idx", float, ())
+        f.variables["aridity_idx"][...] = aridity_index.data
+
+        ppfd_var = f.createVariable("ppfd", float, ("datetimes",))
+        ppfd_var[:] = ppfd
+
+
+def read_faparlim_input(
+    filename: str,
+) -> tuple[xr.Dataset, xr.Dataset, DataArray, DataArray, DataArray, DataArray]:
+    """Method for reading in the data needed for the FaparLimitation constructor."""
+
+    with Dataset(filename, "r") as f:
+        annual_total_A0_subdaily = f.variables["annual_total_A0_subdaily"][...].data
+        annual_total_P = f.variables["annual_total_P"][...].data
+        aridity_index = f.variables["aridity_idx"][...].data
+        annual_mean_ca = f.variables["annual_mean_ca_in_GS"][...].data
+        annual_mean_chi = f.variables["annual_mean_chi_in_GS"][...].data
+        annual_mean_vpd = f.variables["annual_mean_VPD_in_GS"][...].data
+
+    return (
+        annual_total_A0_subdaily,
+        annual_total_P,
+        aridity_index,
+        annual_mean_ca,
+        annual_mean_chi,
+        annual_mean_vpd,
+    )
+
+
+def read_pmodel_faparlim_input(
+    filename: str,
+) -> tuple[
+    DataArray,
+    DataArray,
+    DataArray,
+    DataArray,
+    DataArray,
+    DataArray,
+    DataArray,
+    DataArray,
+    DataArray,
+]:
+    """Reading in the data needed for the from_pmodel FaparLimitation class method."""
+
+    with Dataset(filename, "r") as f:
+        tc = f.variables["tc"][...].data
+        vpd = f.variables["vpd"][...].data
+        co2 = f.variables["co2"][...].data
+        patm = f.variables["patm"][...].data
+        growing_season = f.variables["growing_season"][...].data
+        datetimes = f.variables["datetimes"][...].astype(np.datetime64)
+        precipitation = f.variables["precipitation"][...].data
+        aridity_index = f.variables["aridity_idx"][...].data
+        ppfd = f.variables["ppfd"][...].data
+
+    return (
+        tc,
+        vpd,
+        co2,
+        patm,
+        growing_season,
+        datetimes,
+        precipitation,
+        aridity_index,
+        ppfd,
+    )

tests/regression/phenology/test_phenology_faparlimitation.py

@@ -0,0 +1,170 @@
+"""Test the FaparLimitation class."""
+
+from importlib import resources
+
+import numpy as np
+import pytest
+
+
+@pytest.fixture()
+def annual_data():
+    """Load the input data from netcdf file."""
+
+    from pyrealm_build_data.LAI_in_pyrealm import faparlim_data_io
+
+    datafile = (
+        resources.files("pyrealm_build_data") / "LAI_in_pyrealm/faparlim_input.nc"
+    )
+
+    return faparlim_data_io.read_faparlim_input(datafile)
+
+
+@pytest.mark.parametrize(
+    argnames="exp_faparmax, exp_laimax",
+    argvalues=[
+        (
+            np.array(
+                [
+                    0.00741829,
+                    0.00657795,
+                    0.00663436,
+                    0.00792586,
+                    0.00744838,
+                    0.00790967,
+                    0.00790302,
+                    0.00552947,
+                    0.00616514,
+                    0.00852866,
+                    0.00539698,
+                ]
+            ),
+            np.array(
+                [
+                    0.01489189,
+                    0.01319937,
+                    0.01331292,
+                    0.01591488,
+                    0.01495251,
+                    0.01588223,
+                    0.01586883,
+                    0.01108962,
+                    0.01236845,
+                    0.01713047,
+                    0.0108232,
+                ]
+            ),
+        )
+    ],
+)
+def test_faparlimitation(annual_data, exp_faparmax, exp_laimax):
+    """Regression test for FaparLimitation constructor."""
+
+    from pyrealm.phenology.fapar_limitation import FaparLimitation
+
+    (
+        annual_total_A0_subdaily,
+        annual_total_P,
+        aridity_index,
+        annual_mean_ca,
+        annual_mean_chi,
+        annual_mean_vpd,
+    ) = annual_data
+
+    faparlim = FaparLimitation(
+        annual_total_A0_subdaily,
+        annual_mean_ca,
+        annual_mean_chi,
+        annual_mean_vpd,
+        annual_total_P,
+        aridity_index,
+    )
+
+    assert np.allclose(exp_faparmax, faparlim.fapar_max)
+    assert np.allclose(exp_laimax, faparlim.lai_max)
+
+
+@pytest.fixture()
+def pmodel_data():
+    """Load the input data for the from_pmodel class function from netcdf file."""
+
+    from pyrealm_build_data.LAI_in_pyrealm import faparlim_data_io
+
+    datafile = (
+        resources.files("pyrealm_build_data")
+        / "LAI_in_pyrealm/faparlim_pmodel_input.nc"
+    )
+
+    return faparlim_data_io.read_pmodel_faparlim_input(datafile)
+
+
+@pytest.mark.parametrize(
+    argnames="exp_faparmax, exp_laimax",
+    argvalues=[
+        (
+            np.array(
+                [
+                    0.00890042,
+                    0.00748109,
+                    0.00538652,
+                    0.00979542,
+                    0.00671602,
+                    0.00627187,
+                    0.00778598,
+                    0.00556823,
+                    0.00648213,
+                    0.00921882,
+                    0.00805771,
+                ]
+            ),
+            np.array(
+                [
+                    0.01788053,
+                    0.01501843,
+                    0.01080216,
+                    0.01968741,
+                    0.01347735,
+                    0.01258324,
+                    0.01563289,
+                    0.01116758,
+                    0.01300646,
+                    0.01852314,
+                    0.0161807,
+                ]
+            ),
+        )
+    ],
+)
+def test_faparlimitation_frompmodel(pmodel_data, exp_faparmax, exp_laimax):
+    """Regression test for from_pmodel FaparLimitation class method."""
+
+    from pyrealm.phenology.fapar_limitation import FaparLimitation
+    from pyrealm.pmodel import PModel, PModelEnvironment
+
+    (
+        tc,
+        vpd,
+        co2,
+        patm,
+        growing_season,
+        datetimes,
+        precipitation,
+        aridity_index,
+        ppfd,
+    ) = pmodel_data
+
+    env = PModelEnvironment(tc=tc, vpd=vpd, co2=co2, patm=patm)
+
+    pmodel = PModel(
+        env=env,
+        reference_kphio=1 / 8,
+        method_kphio="temperature",
+    )
+
+    pmodel.estimate_productivity(fapar=np.ones_like(env.ca), ppfd=ppfd)
+
+    faparlim = FaparLimitation.from_pmodel(
+        pmodel, growing_season, datetimes, precipitation, aridity_index
+    )
+
+    assert np.allclose(exp_faparmax, faparlim.fapar_max)
+    assert np.allclose(exp_laimax, faparlim.lai_max)