diff --git a/docs/source/refs.bib b/docs/source/refs.bib index 7a4abd0d..d356ed95 100644 --- a/docs/source/refs.bib +++ b/docs/source/refs.bib @@ -947,3 +947,30 @@ @misc{zotero-2443a howpublished = {http://localhost:8888/lab/tree/source/users/demography/crown.md\#crown-radius-values}, file = {/Users/dorme/Zotero/storage/LHHVRPS9/crown.html} } + +@article{Pury&Farquhar:1997, + title = {Simple scaling of photosynthesis from leaves to canopies without the errors of big-leaf models}, + volume = {20}, + issn = {}, + url = {https://onlinelibrary.wiley.com/doi/10.1111/j.1365-3040.1997.00094.x}, + doi = {}, + abstract = {In big-leaf models of canopy photosynthesis, the Rubisco activity per unit ground area is taken as the sum of activities per unit leaf area within the canopy, and electron transport capacity is similarly summed. Such models overestimate rates of photosynthesis and require empirical curvature factors in the response to irradiance. We show that, with any distribution of leaf nitrogen within the canopy (including optimal), the required curvature factors are not constant but vary with canopy leaf area index and leaf nitrogen content. We further show that the underlying reason is the difference between the time-averaged and instantaneous distributions of absorbed irradiance, caused by penetration of sunflecks and the range of leaf angles in canopies. + These errors are avoided in models that treat the canopy in terms of a number of layers – the multi-layer models. We present an alternative to the multi-layer model: by separately integrating the sunlit and shaded leaf fractions of the canopy, a single layered sun/shade model is obtained, which is as accurate and simpler. The model is a scaled version of a leaf model as distinct from an integrative approach.}, + language = {}, + number = {}, + urldate = {}, + journal = {}, + author = {}, + month = {}, + year = {}, + note = {}, + pages = {}, + file = {} + } + +@article{Ryu_et_al:2011, + title = {Integration of MODIS land and atmosphere products with a coupled-process model to estimate gross primary productivity and evapotranspiration from 1 km to global scales}, + volume = {}, + + +} diff --git a/profiling/benchmark-fails.csv b/profiling/benchmark-fails.csv index d427274d..c8d93c62 100644 --- a/profiling/benchmark-fails.csv +++ b/profiling/benchmark-fails.csv @@ -1,7 +1,7 @@ ncalls,tottime,tottime_percall,cumtime,cumtime_percall,timestamp,filename,lineno,function,process_id,label,ignore_result,ignore_justification,is_incoming,tottime_percall_target,cumtime_percall_target,tottime_target,cumtime_target,relative_tottime_percall,relative_cumtime_percall,relative_tottime,relative_cumtime -2,0.13,0.065,0.13,0.065,2024-05-07T16:26:25,pyrealm/pmodel/functions.py,644,calc_co2_to_ca,pyrealm/pmodel/functions.py:calc_co2_to_ca,7c3fe42,False,,False,0.065,0.065,0.13,0.13,1.0,1.0,1.0,1.0 -1,0.016,0.016,0.016,0.016,2024-05-07T16:26:25,pyrealm/pmodel/isotopes.py,102,calc_c4_discrimination,pyrealm/pmodel/isotopes.py:CalcCarbonIsotopes.calc_c4_discrimination,7c3fe42,False,,False,0.016,0.016,0.016,0.016,1.0,1.0,1.0,1.0 -1,0.058,0.058,0.359,0.359,2024-05-07T16:26:25,pyrealm/core/hygro.py,259,calc_psychrometric_constant,pyrealm/core/hygro.py:calc_psychrometric_constant,7c3fe42,False,,False,0.058,0.359,0.058,0.359,1.0,1.0,1.0,1.0 -2,0.169,0.085,0.169,0.085,2024-06-05T11:06:21,pyrealm/pmodel/functions.py,644,calc_co2_to_ca,pyrealm/pmodel/functions.py:calc_co2_to_ca,9760db6,False,,True,0.065,0.065,0.13,0.13,1.3076923076923077,1.3076923076923077,1.3,1.3 -1,0.017,0.017,0.017,0.017,2024-06-05T11:06:21,pyrealm/pmodel/isotopes.py,102,calc_c4_discrimination,pyrealm/pmodel/isotopes.py:CalcCarbonIsotopes.calc_c4_discrimination,9760db6,False,,True,0.016,0.016,0.016,0.016,1.0625,1.0625,1.0625,1.0625 -1,0.061,0.061,0.37,0.37,2024-06-05T11:06:21,pyrealm/core/hygro.py,259,calc_psychrometric_constant,pyrealm/core/hygro.py:calc_psychrometric_constant,9760db6,False,,True,0.058,0.359,0.058,0.359,1.0517241379310345,1.0306406685236769,1.0517241379310345,1.0306406685236769 +4,0.101,0.025,0.647,0.162,2024-05-07T16:26:25,pyrealm/core/water.py,67,calc_density_h2o_fisher,pyrealm/core/water.py:calc_density_h2o_fisher,7c3fe42,False,,False,0.025,0.162,0.101,0.647,1.0,1.0,1.0,1.0 +5,0.005,0.001,1.721,0.344,2024-05-07T16:26:25,pyrealm/core/water.py,126,calc_density_h2o,pyrealm/core/water.py:calc_density_h2o,7c3fe42,False,,False,0.001,0.344,0.005,1.721,1.0,1.0,1.0,1.0 +17,1.463,0.086,1.574,0.093,2024-05-07T16:26:25,pyrealm/core/utilities.py,339,evaluate_horner_polynomial,pyrealm/core/utilities.py:evaluate_horner_polynomial,7c3fe42,False,,False,0.086,0.093,1.463,1.574,1.0,1.0,1.0,1.0 +4,0.09,0.022,0.937,0.234,2024-06-05T11:31:06,pyrealm/core/water.py,67,calc_density_h2o_fisher,pyrealm/core/water.py:calc_density_h2o_fisher,62d7a74,False,,True,0.025,0.162,0.101,0.647,0.8799999999999999,1.4444444444444444,0.891089108910891,1.4482225656877898 +5,0.005,0.001,1.936,0.387,2024-06-05T11:31:06,pyrealm/core/water.py,126,calc_density_h2o,pyrealm/core/water.py:calc_density_h2o,62d7a74,False,,True,0.001,0.344,0.005,1.721,1.0,1.1250000000000002,1.0,1.124927367809413 +17,1.7,0.1,1.805,0.106,2024-06-05T11:31:06,pyrealm/core/utilities.py,339,evaluate_horner_polynomial,pyrealm/core/utilities.py:evaluate_horner_polynomial,62d7a74,False,,True,0.086,0.093,1.463,1.574,1.1627906976744187,1.139784946236559,1.161995898838004,1.1467598475222363 diff --git a/profiling/call-graph.svg b/profiling/call-graph.svg index 411bb0ca..526cddb8 100644 --- a/profiling/call-graph.svg +++ b/profiling/call-graph.svg @@ -4,20 +4,20 @@ - + %3 - + 73 - -utilities:113:summarize_attrs -6.64% -(0.71%) -6× + +utilities:113:summarize_attrs +6.39% +(0.67%) +6× @@ -25,139 +25,139 @@ 885 - -nanfunctions:236:nanmin -1.11% -(0.00%) -47× + +nanfunctions:236:nanmin +0.89% +(0.00%) +47× 73->885 - - -0.65% -38× + + +0.63% +38× 886 - -nanfunctions:952:nanmean -4.51% -(0.13%) -38× + +nanfunctions:952:nanmean +4.35% +(0.12%) +38× 73->886 - - -4.51% -38× + + +4.35% +38× 891 - -nanfunctions:369:nanmax -0.63% -(0.00%) -40× + +nanfunctions:369:nanmax +0.59% +(0.00%) +40× 73->891 - - -0.63% -38× + + +0.59% +38× - - -979 - - -~:0:<method 'reduce' of 'numpy.ufunc' objects> -2.96% -(2.96%) -3271× + + +971 + + +~:0:<method 'reduce' of 'numpy.ufunc' objects> +2.69% +(2.69%) +3271× - - -885->979 - - -1.11% -45× + + +885->971 + + +0.89% +45× - + 884 - - -nanfunctions:68:_replace_nan -3.24% -(1.32%) -40× + + +nanfunctions:68:_replace_nan +3.10% +(1.31%) +40× - + 886->884 - - -3.24% -38× + + +3.10% +38× - + 1009 - - -fromnumeric:2177:sum -1.13% -(0.00%) -76× + + +fromnumeric:2177:sum +1.12% +(0.00%) +76× - + 886->1009 - - -1.13% -76× + + +1.12% +76× - - -891->979 - - -0.62% -40× + + +891->971 + + +0.59% +40× 91 - -functions:14:calc_ftemp_arrh -1.87% -(1.87%) -13× + +functions:14:calc_ftemp_arrh +1.81% +(1.81%) +13× @@ -165,90 +165,90 @@ 92 - -functions:119:calc_ftemp_inst_vcmax -3.83% -(3.01%) -3× + +functions:119:calc_ftemp_inst_vcmax +3.70% +(2.89%) +3× 92->91 - - -0.82% -3× + + +0.81% +3× 93 - -functions:256:calc_gammastar -0.75% -(0.59%) -2× + +functions:256:calc_gammastar +0.80% +(0.64%) +2× 93->91 - - -0.16% -2× + + +0.16% +2× 94 - -functions:354:calc_kmm -1.22% -(0.65%) -2× + +functions:354:calc_kmm +1.34% +(0.78%) +2× 94->91 - - -0.57% -4× + + +0.55% +4× 95 - -subdaily:218:__init__ -5.09% -(1.08%) -1× + +subdaily:218:__init__ +4.92% +(1.03%) +1× 95->91 - - -0.32% -4× + + +0.30% +4× 97 - + pmodel:358:estimate_productivity -9.50% -(3.17%) +9.12% +(3.09%) 3× @@ -256,29 +256,29 @@ 95->97 - - -0.15% + + +0.14% 1× 98 - -functions:197:calc_ftemp_kphio -1.56% -(1.42%) -4× + +functions:197:calc_ftemp_kphio +1.64% +(1.52%) +4× 95->98 - - -0.59% + + +0.52% 1× @@ -287,8 +287,8 @@ pmodel:186:__init__ -7.81% -(1.06%) +7.91% +(0.98%) 3× @@ -296,396 +296,396 @@ 95->99 - + -0.13% +0.12% 1× 100 - -pmodel_environment:76:__init__ -22.84% -(0.02%) -2× + +pmodel_environment:76:__init__ +23.66% +(0.02%) +2× 95->100 - - -0.55% -1× + + +0.55% +1× 127 - -optimal_chi:95:__init__ -3.22% -(0.00%) -4× + +optimal_chi:95:__init__ +3.30% +(0.00%) +4× 95->127 - - -1.19% + + +1.16% 1× 130 - -optimal_chi:253:estimate_chi -3.24% -(3.24%) -4× + +optimal_chi:253:estimate_chi +3.26% +(3.26%) +4× 95->130 - - -0.88% -1× + + +0.89% +1× 97->92 - - -3.83% -3× + + +3.70% +3× 96 - -functions:77:calc_ftemp_inst_rd -1.12% -(1.12%) -3× + +functions:77:calc_ftemp_inst_rd +1.01% +(1.01%) +3× 97->96 - - -1.12% -3× + + +1.01% +3× 996 - -numeric:2170:allclose -1.37% -(0.02%) -3× + +numeric:2170:allclose +1.33% +(0.02%) +3× 97->996 - - -1.37% -3× + + +1.33% +3× 1001 - -fromnumeric:2100:clip -1.88% -(0.02%) -4398× + +fromnumeric:2100:clip +1.88% +(0.02%) +4398× 98->1001 - - -0.14% -4× + + +0.12% +4× 99->98 - - -0.97% -3× + + +1.12% +3× 99->127 - - -2.03% -3× + + +2.14% +3× 124 - -jmax_limitation:75:__init__ -3.75% -(0.00%) -3× + +jmax_limitation:75:__init__ +3.68% +(0.00%) +3× 99->124 - - -3.75% -3× + + +3.68% +3× 100->885 - - -0.29% -4× + + +0.16% +4× 100->93 - - -0.75% -2× + + +0.80% +2× 100->94 - - -1.22% -2× + + +1.34% +2× 101 - -functions:310:calc_ns_star -17.49% -(0.04%) -2× + +functions:310:calc_ns_star +18.66% +(0.04%) +2× 100->101 - - -17.49% -2× + + +18.66% +2× 102 - -functions:644:calc_co2_to_ca -0.72% -(0.72%) -2× + +functions:644:calc_co2_to_ca +0.59% +(0.59%) +2× 100->102 - - -0.72% -2× + + +0.59% +2× - + -147 +148 - -utilities:213:bounds_checker -3.02% -(2.98%) -13× + +utilities:213:bounds_checker +2.75% +(2.72%) +13× - + -100->147 - - -2.35% -8× +100->148 + + +2.11% +8× 127->130 - - -2.36% -3× + + +2.37% +3× 132 - -optimal_chi:385:estimate_chi -0.86% -(0.75%) -1× + +optimal_chi:385:estimate_chi +0.92% +(0.80%) +1× 127->132 - - -0.86% -1× + + +0.92% +1× - - -961 - - -numeric:2249:isclose -1.35% -(0.95%) -9× + + +953 + + +numeric:2249:isclose +1.31% +(0.94%) +9× - - -996->961 - - -1.35% -3× + + +996->953 + + +1.31% +3× - - -908 - - -fromnumeric:53:_wrapfunc -1.86% -(0.02%) -4526× + + +897 + + +fromnumeric:53:_wrapfunc +1.86% +(0.02%) +4526× - - -1001->908 - - -1.85% -4398× + + +1001->897 + + +1.85% +4398× 125 - -jmax_limitation:127:wang17 -3.75% -(3.75%) -3× + +jmax_limitation:127:wang17 +3.68% +(3.68%) +3× 124->125 - - -3.75% -3× + + +3.68% +3× 105 - -water:181:calc_viscosity_h2o -17.45% -(14.55%) -4× + +water:181:calc_viscosity_h2o +18.62% +(14.46%) +4× 101->105 - - -17.45% -4× + + +18.62% +4× 104 - -water:126:calc_density_h2o -7.39% -(0.02%) -5× + +water:126:calc_density_h2o +8.40% +(0.02%) +5× 105->104 - - -2.90% -4× + + +4.15% +4× 103 - -water:67:calc_density_h2o_fisher -2.75% -(0.44%) -4× + +water:67:calc_density_h2o_fisher +4.07% +(0.39%) +4× @@ -791,195 +791,179 @@ -148 +149 - -utilities:339:evaluate_horner_polynomial -6.72% -(6.23%) -17× + +utilities:339:evaluate_horner_polynomial +7.83% +(7.37%) +17× - + -103->148 - - -2.31% -12× - - - -968 - - -numeric:67:zeros_like -0.51% -(0.51%) -24× - - - - - -148->968 - - -0.49% -17× +103->149 + + +3.68% +12× 104->885 - - -0.18% -5× + + +0.11% +5× 104->103 - - -2.75% -4× + + +4.07% +4× -150 +151 - -water:12:calc_density_h2o_chen -4.44% -(0.90%) -1× + +water:12:calc_density_h2o_chen +4.20% +(0.85%) +1× - + -104->150 - - -4.44% -1× +104->151 + + +4.20% +1× - - -150->148 - - -3.55% -4× + + +151->149 + + +3.35% +4× 106 - -evap:74:__post_init__ -8.21% -(1.19%) -1× + +evap:74:__post_init__ +7.78% +(1.15%) +1× 106->104 - - -4.49% -1× + + +4.25% +1× - + -1866 +1865 - -hygro:188:calc_saturation_vapour_pressure_slope -0.64% -(0.64%) -1× + +hygro:188:calc_saturation_vapour_pressure_slope +0.62% +(0.62%) +1× + + + + + +106->1865 + + +0.62% +1× + + + +1866 + + +hygro:210:calc_enthalpy_vaporisation +0.56% +(0.56%) +2× + 106->1866 - - -0.64% -1× + + +0.28% +1× + 1867 +1867 - -hygro:210:calc_enthalpy_vaporisation -0.60% -(0.60%) -2× + +hygro:259:calc_psychrometric_constant +1.56% +(0.25%) +1× + 106->1867 - - -0.30% -1× - - - -1868 - - -hygro:259:calc_psychrometric_constant -1.58% -(0.26%) -1× - - - - - -106->1868 - - -1.58% -1× - - - -1868->1867 - - -0.30% -1× - - - -149 - - -hygro:227:calc_specific_heat -1.01% -(0.03%) -1× + + +1.48% +1× + + + +1867->1866 + + +0.28% +1× + + + +150 + + +hygro:227:calc_specific_heat +0.96% +(0.03%) +1× - - -1868->149 - - -1.01% -1× + + +1867->150 + + +0.96% +1× @@ -987,8 +971,8 @@ isotopes:42:__init__ -1.11% -(0.79%) +1.05% +(0.74%) 2× @@ -997,37 +981,37 @@ 108 - -test_profiling_pmodel:7:test_profiling_pmodel -26.87% -(0.17%) -1× + +test_profiling_pmodel:7:test_profiling_pmodel +26.24% +(0.16%) +1× 108->97 - - -9.35% -2× + + +8.98% +2× 108->99 - - -7.69% -2× + + +7.79% +2× 108->107 - - -1.11% -2× + + +1.05% +2× @@ -1035,7 +1019,7 @@ isotopes:226:summarize -2.03% +1.95% (0.00%) 2× @@ -1044,10 +1028,10 @@ 108->113 - - -2.03% -2× + + +1.95% +2× @@ -1055,7 +1039,7 @@ pmodel:456:summarize -2.78% +2.66% (0.00%) 2× @@ -1064,10 +1048,10 @@ 108->123 - - -2.78% -2× + + +2.66% +2× @@ -1075,8 +1059,8 @@ competition:182:__init__ -1.47% -(0.38%) +1.46% +(0.37%) 1× @@ -1084,10 +1068,10 @@ 108->136 - - -1.47% -1× + + +1.46% +1× @@ -1095,7 +1079,7 @@ competition:289:summarize -1.83% +1.78% (0.00%) 2× @@ -1104,174 +1088,162 @@ 108->140 - - -1.83% -2× + + +1.78% +2× 113->73 - - -2.03% -2× + + +1.95% +2× 123->73 - - -2.78% -2× + + +2.66% +2× 137 - -competition:62:calculate_tree_proportion -0.62% -(0.54%) -1× + +competition:62:calculate_tree_proportion +0.63% +(0.55%) +1× 136->137 - - -0.62% -1× + + +0.63% +1× 140->73 - - -1.83% -2× + + +1.78% +2× 126 - -conftest:11:pmodel_profile_data -22.76% -(0.00%) -2× + +conftest:11:pmodel_profile_data +23.55% +(0.00%) +2× 126->100 - - -22.29% -1× + + +23.11% +1× - + -142 - - -evap:109:estimate_aet -12.74% -(12.32%) -1463× +141 + + +python:187:pytest_pyfunc_call +75.86% +(0.19%) +3× - + -142->1001 - - -0.34% -1463× +141->108 + + +26.24% +1× - + -143 - - -solar:71:__post_init__ -19.51% -(19.02%) -1× +1659 + + +test_profiling_subdaily:7:test_profiling_subdaily +4.98% +(0.01%) +1× - + -143->1001 - - -0.49% -2× +141->1659 + + +4.98% +1× - + -144 - - -splash:58:__init__ -29.88% -(0.00%) -1× +1660 + + +test_profiling_splash:6:test_profile_splash +44.46% +(0.00%) +1× - + -144->147 - - -0.66% -5× - - - -1863 - - -pressure:10:calc_patm -1.50% -(1.50%) -1× - - +141->1660 + + +44.46% +1× - - -144->1863 - - -1.50% -1× + + +1659->95 + + +4.92% +1× - + 145 - - -splash:224:estimate_daily_water_balance -14.22% -(0.69%) -1463× + + +splash:58:__init__ +29.56% +(0.00%) +1× - - -145->1001 - - -0.69% -2926× + + +1660->145 + + +29.56% +1× @@ -1305,13 +1277,13 @@ 0.13% 1× - - -149->148 - - -0.86% -1× + + +1660->890 + + +7.32% +1× @@ -1325,677 +1297,694 @@ - - -188->108 - - -26.87% -1× + + +1660->991 + + +7.58% +1× - - -1663 - - -test_profiling_subdaily:7:test_profiling_subdaily -5.15% -(0.01%) -1× + + +143 + + +evap:109:estimate_aet +13.06% +(12.63%) +1463× - - -188->1663 - - -5.15% -1× + + +143->1001 + + +0.35% +1463× - - -1664 - - -test_profiling_splash:6:test_profile_splash -44.45% -(0.00%) -1× + + +144 + + +solar:71:__post_init__ +19.64% +(19.14%) +1× - - -188->1664 - - -44.45% -1× - - - -1663->95 - - -5.09% -1× + + +144->1001 + + +0.49% +2× - - -1664->144 - - -29.88% -1× + + +145->148 + + +0.64% +5× - - -890 - - -splash:110:estimate_initial_soil_moisture -7.15% -(0.01%) -1× + + +1864 + + +pressure:10:calc_patm +1.50% +(1.50%) +1× - - -1664->890 - - -7.15% -1× + + +145->1864 + + +1.50% +1× - - -991 - - -splash:287:calculate_soil_moisture -7.41% -(0.11%) -1× + + +146 + + +splash:224:estimate_daily_water_balance +14.55% +(0.68%) +1463× - - -1664->991 - - -7.41% -1× + + +146->1001 + + +0.70% +2926× - - -231 - - -~:0:<built-in method builtins.next> -23.26% -(0.00%) -337× + + +146->143 + + +13.06% +1463× + + + +150->1001 + + +0.12% +1× + + + +150->149 + + +0.80% +1× + + + +255 + + +~:0:<built-in method builtins.next> +24.04% +(0.00%) +337× - - -231->126 - - -22.76% -2× + + +255->126 + + +23.55% +2× - + 648 - - -fixtures:592:_get_active_fixturedef -23.01% -(0.00%) -18× + + +fixtures:592:_get_active_fixturedef +23.84% +(0.00%) +18× - - -1204 - - -fixtures:620:_compute_fixture_value -23.01% -(0.00%) -6× + + +1203 + + +fixtures:620:_compute_fixture_value +23.84% +(0.00%) +6× - - -648->1204 - - -23.01% -6× + + +648->1203 + + +23.84% +6× - - -1121 - - -fixtures:1041:execute -23.01% -(0.00%) -6× + + +1120 + + +fixtures:1041:execute +23.84% +(0.00%) +6× - - -1204->1121 - - -23.01% -6× + + +1203->1120 + + +23.84% +6× - - -939 - - -~:0:<built-in method numpy.array> -1.92% -(1.92%) -96× + + +928 + + +~:0:<built-in method numpy.array> +1.79% +(1.79%) +96× - - -884->939 - - -1.91% -40× + + +884->928 + + +1.79% +40× - - -905 - - -fromnumeric:71:_wrapreduction -1.26% -(0.03%) -3043× + + +894 + + +fromnumeric:71:_wrapreduction +1.24% +(0.03%) +3043× - - -1009->905 - - -1.13% -76× + + +1009->894 + + +1.11% +76× - - -890->145 - - -7.14% -732× + + +890->146 + + +7.31% +732× - - -905->979 - - -1.22% -3027× + + +894->971 + + +1.20% +3027× - - -918 - - -~:0:<method 'clip' of 'numpy.ndarray' objects> -1.83% -(0.01%) -4401× + + +907 + + +~:0:<method 'clip' of 'numpy.ndarray' objects> +1.83% +(0.01%) +4401× - - -908->918 - - -1.83% -4398× + + +897->907 + + +1.83% +4398× - + 1026 - - -_methods:90:_clip -1.82% -(1.82%) -4401× + + +_methods:90:_clip +1.82% +(1.82%) +4401× - - -918->1026 - - -1.82% -4401× + + +907->1026 + + +1.82% +4401× - - -991->145 - - -7.08% -731× + + +991->146 + + +7.24% +731× - + 1059 - - -_callers:52:_multicall -99.98% -(0.05%) -57× + + +_callers:52:_multicall +99.98% +(0.05%) +57× + + + + + +1059->141 + + +75.86% +3× + + + +1138 + + +runner:160:pytest_runtest_call +75.91% +(0.00%) +3× - + -1059->188 - - -76.65% -3× +1059->1138 + + +75.91% +3× - + -1150 - - -runner:160:pytest_runtest_call -76.69% -(0.00%) -3× +1200 + + +fixtures:1108:pytest_fixture_setup +23.84% +(0.00%) +3× - + -1059->1150 - - -76.69% -3× +1059->1200 + + +23.84% +3× - + -1201 - - -fixtures:1108:pytest_fixture_setup -23.01% -(0.00%) -3× +1213 + + +runner:155:pytest_runtest_setup +23.84% +(0.00%) +3× - + -1059->1201 - - -23.01% -3× +1059->1213 + + +23.84% +3× - + -1214 - - -runner:155:pytest_runtest_setup -23.02% -(0.00%) -3× +1119 + + +python:1790:runtest +75.91% +(0.00%) +3× - - -1059->1214 - - -23.02% -3× + + +1138->1119 + + +75.91% +3× - - -1120 - - -python:1790:runtest -76.69% -(0.00%) -3× + + +1207 + + +fixtures:886:call_fixture_func +23.84% +(0.00%) +3× - - -1150->1120 - - -76.69% -3× + + +1200->1207 + + +23.84% +3× - - -1208 - - -fixtures:886:call_fixture_func -23.01% -(0.00%) -3× + + +1140 + + +runner:478:setup +23.84% +(0.00%) +3× - - -1201->1208 - - -23.01% -3× + + +1213->1140 + + +23.84% +3× - - -1152 - - -runner:478:setup -23.02% -(0.00%) -3× + + +1395 + + +_hooks:486:__call__ +99.98% +(0.00%) +57× - - -1214->1152 - - -23.02% -3× + + +1119->1395 + + +75.91% +3× - - -1396 - - -_hooks:486:__call__ -99.98% -(0.00%) -57× + + +1394 + + +_manager:110:_hookexec +99.98% +(0.00%) +57× - + + +1395->1394 + + +99.98% +33× + + -1120->1396 - - -76.69% -3× +1120->1395 + + +23.84% +3× - - -1395 - - -_manager:110:_hookexec -99.98% -(0.00%) -57× + + +1121 + + +runner:111:pytest_runtest_protocol +100.00% +(0.00%) +3× - - -1396->1395 - - -99.98% -33× - - - -1121->1396 - - -23.01% -3× - - + -1122 +1212 - -runner:111:pytest_runtest_protocol -100.00% -(0.00%) -3× + +runner:119:runtestprotocol +99.99% +(0.01%) +3× - + + +1121->1212 + + +99.99% +3× + + -1213 - - -runner:119:runtestprotocol -99.99% -(0.01%) -3× +1122 +1122 + + +runner:219:call_and_report +99.98% +(0.00%) +9× - - -1122->1213 - - -99.99% -3× + + +1212->1122 + + +99.98% +9× 1123 - -runner:219:call_and_report -99.98% -(0.00%) -9× + +runner:247:call_runtest_hook +99.97% +(0.00%) +9× - - -1213->1123 - - -99.98% -9× + + +1122->1123 + + +99.97% +9× - + -1124 +1218 - -runner:247:call_runtest_hook -99.97% -(0.00%) -9× + +runner:318:from_call +99.97% +(0.00%) +9× - + -1123->1124 - - -99.97% -9× +1123->1218 + + +99.97% +9× - - -1219 - - -runner:318:from_call -99.97% -(0.00%) -9× + + +1217 + + +runner:262:<lambda> +99.97% +(0.00%) +9× - - -1124->1219 - - -99.97% -9× + + +1218->1217 + + +99.97% +9× - - -1218 - - -runner:262:<lambda> -99.97% -(0.00%) -9× + + +1139 + + +python:1794:setup +23.84% +(0.00%) +3× - - -1219->1218 - - -99.97% -9× - - + -1151 - - -python:1794:setup -23.02% -(0.00%) -3× +1194 + + +fixtures:561:_fillfixtures +23.84% +(0.00%) +3× - - -1195 - - -fixtures:561:_fillfixtures -23.02% -(0.00%) -3× - - - - - -1151->1195 - - -23.02% -3× + + +1139->1194 + + +23.84% +3× - - -1200 - - -fixtures:568:getfixturevalue -23.01% -(0.00%) -9× + + +1199 + + +fixtures:568:getfixturevalue +23.84% +(0.00%) +9× - - -1195->1200 - - -23.01% -9× - - + -1152->1151 - - -23.02% -3× +1194->1199 + + +23.84% +9× - - -1200->648 - - -23.01% -9× + + +1140->1139 + + +23.84% +3× - - -1208->231 - - -23.01% -2× + + +1199->648 + + +23.84% +9× - - -1218->1396 - - -99.97% -9× + + +1207->255 + + +23.84% +2× - - -1395->1059 - - -99.98% -33× + + +1217->1395 + + +99.97% +9× + + + +1394->1059 + + +99.98% +33× diff --git a/pyrealm/constants/two_leaf_canopy.py b/pyrealm/constants/two_leaf_canopy.py new file mode 100644 index 00000000..4d1d9684 --- /dev/null +++ b/pyrealm/constants/two_leaf_canopy.py @@ -0,0 +1,28 @@ +"""To do.""" + +from dataclasses import dataclass + +import numpy as np + +from pyrealm.constants import ConstantsClass + + +@dataclass(frozen=True) +class TwoLeafConst(ConstantsClass): + """Pyrealm two leaf canopy model constants dataclass.""" + + # two leaf canopy model constants + k_PA0 = 101325 + """Reference standard pressure""" + k_fa: float = 0.426 # needs citation + """scattering coefficient of PAR in the atmosphere, dimensionless""" + k_sigma: float = 0.15 # needs citation + """leaf scattering coefficient of PAR (relections and transmissivity), + dimensionless""" + k_rho_cd: float = 0.036 # needs citation + """canopy reflection coefficient for diffuse PAR, dimensionless""" + k_kd_prime: float = 0.719 # needs citation + """diffuse and scattered diffuse PAR extinction coefficient, dimensionless""" + + k_sol_obs_angle: float = 1 / (2 * np.pi) # 1 degree in rads, needs a citation + """ solar obscurity angle, radians""" diff --git a/pyrealm/pmodel/two_leaf_irradience.py b/pyrealm/pmodel/two_leaf_irradience.py new file mode 100644 index 00000000..b4612044 --- /dev/null +++ b/pyrealm/pmodel/two_leaf_irradience.py @@ -0,0 +1,1041 @@ +"""To do.""" + +import numpy as np +from numpy.typing import NDArray + +from pyrealm.constants.core_const import CoreConst +from pyrealm.constants.two_leaf_canopy import TwoLeafConst +from pyrealm.pmodel import PModel +from pyrealm.pmodel.optimal_chi import OptimalChiABC +from pyrealm.pmodel.subdaily import SubdailyPModel + + +class TwoLeafIrradience: + """Running the two leaf, two stream model within Pyrealm. + + This class implements the methodology of Pury and Farquhar (1997) + :cite:p:`Pury&Farquhar:1997` two leaf, two stream model. This model is chosen to + provide a better representation than the big leaf model and to align closely to + the workings of the ``BESS`` model :cite:alp:`Ryu_et_al:2011`. + + Args: + beta_angle (NDArray): Array of beta angles (radians). + ppfd (NDArray): Array of photosynthetic photon flux density values. + leaf_area_index (NDArray): Array of leaf area index values. + patm (NDArray): Array of atmospheric pressure values. + """ + + def __init__( + self, + beta_angle: NDArray, + ppfd: NDArray, + leaf_area_index: NDArray, + patm: NDArray, + constants: TwoLeafConst = TwoLeafConst(), + ): + self.beta_angle = beta_angle + self.ppfd = ppfd + self.leaf_area_index = leaf_area_index + self.patm = patm + + self.pa0: float = constants.k_PA0 + self.k_fa: float = constants.k_fa + self.k_sigma: float = constants.k_sigma + self.k_rho_cd: float = constants.k_rho_cd + self.k_kd_prime: float = constants.k_kd_prime + self.k_sol_obs_angle: float = constants.k_sol_obs_angle + + self.kb: NDArray + self.kb_prime: NDArray + self.fd: NDArray + self.rho_h: float + self.rho_cb: NDArray + self.I_d: NDArray + self.T_b: NDArray + self.I_c: NDArray + self.Isun_beam: NDArray + self.Isun_diffuse: NDArray + self.Isun_scattered: NDArray + self.I_csun: NDArray + self.I_cshade: NDArray + + self.arrays = [self.beta_angle, self.ppfd, self.leaf_area_index, self.patm] + self.shapes_agree: bool = self._check_input_consistency() + self.no_NaNs: bool = self._check_for_NaN() + self.no_negatives: bool = self._check_for_negative_values() + self.pass_checks = all([self.shapes_agree, self.no_NaNs, self.no_negatives]) + + def _check_input_consistency(self) -> bool: + """Check if input arrays have consistent shapes. + + Returns: + bool: True if all input arrays have the same shape, False otherwise. + """ + arrays = [self.beta_angle, self.ppfd, self.leaf_area_index, self.patm] + try: + shapes = [array.shape for array in arrays] + if not all(shape == shapes[0] for shape in shapes): + raise ValueError("Input arrays have inconsistent shapes.") + return True + except Exception as e: + print(f"Error in input consistency check: {e}") + return False + + def _check_for_NaN(self) -> bool: + """Tests for any NaN in input arrays.""" + arrays = [self.beta_angle, self.ppfd, self.leaf_area_index, self.patm] + try: + for array in arrays: + if np.isnan(array).any(): + raise ValueError("Nan data in input array") + return True + except Exception as e: + print(f"Error in input NaN check: {e}") + return False + + def _check_for_negative_values(self) -> bool: + """Tests for negative values in arrays.""" + try: + for array in [self.ppfd, self.leaf_area_index, self.patm]: + if not (array >= 0).all(): + raise ValueError("Input arrays contain negative values.") + return True + except Exception as e: + print(f"Error in input negative number check: {e}") + return False + + def calc_absorbed_irradience(self) -> None: + """Calculate absorbed irradiance for sunlit and shaded leaves.""" + + self.kb = beam_extinction_coeff(self.beta_angle, self.k_sol_obs_angle) + self.kb_prime = scattered_beam_extinction_coeff( + self.beta_angle, self.k_sol_obs_angle + ) + self.fd = fraction_of_diffuse_rad( + self.patm, self.pa0, self.beta_angle, self.k_fa + ) + self.rho_h = beam_irradience_h_leaves(self.k_sigma) + self.rho_cb = beam_irrad_unif_leaf_angle_dist(self.rho_h, self.kb) + self.I_d = diffuse_radiation(self.fd, self.ppfd) + self.I_b = beam_irradience(self.ppfd, self.fd) + self.I_c = canopy_irradience( + self.rho_cb, + self.I_b, + self.kb_prime, + self.I_d, + self.leaf_area_index, + self.k_rho_cd, + ) + self.Isun_beam = sunlit_beam_irrad( + self.I_b, self.k_sigma, self.kb, self.leaf_area_index + ) + self.Isun_diffuse = sunlit_diffuse_irrad( + self.I_d, self.k_rho_cd, self.k_kd_prime, self.kb, self.leaf_area_index + ) + self.Isun_scattered = sunlit_scattered_irrad( + self.I_b, + self.rho_cb, + self.kb_prime, + self.kb, + self.leaf_area_index, + self.k_sigma, + ) + self.I_csun = sunlit_absorbed_irrad( + self.Isun_beam, self.Isun_diffuse, self.Isun_scattered + ) + self.I_cshade = shaded_absorbed_irrad( + self.beta_angle, self.k_sol_obs_angle, self.I_c, self.I_csun + ) + + +class TwoLeafAssimilation: + """A class to estimate gross primary production (``GPP``) using a two-leaf approach. + + This class integrates a photosynthesis model (:class:`~pyrealm.pmodel.pmodel.PModel` + or :class:`~pyrealm.pmodel.subdaily.SubdailyPModel`) and irradiance data + (:class:`~pyrealm.pmodel.two_leaf_irradience.TwoLeafIrradience`) to compute various + canopy photosynthetic properties and ``GPP`` estimates. + + Attributes: + pmodel (PModel | SubdailyPModel): The photosynthesis model used for + assimilation. + irrad (TwoLeafIrradience): Irradiance data required for ``GPP`` calculations. + leaf_area_index (NDArray): Array of leaf area index values. + vcmax_pmod (NDArray): Maximum carboxylation rate at current conditions. + vcmax25_pmod (NDArray): Maximum carboxylation rate at 25°C. + optchi_obj (OptimalChiABC): Object containing optimization of :math:`chi`. + core_const (CoreConst): Core constants retrieved from the PModel. + gpp (NDArray): Estimated gross primary production. + """ + + def __init__( + self, + pmodel: PModel | SubdailyPModel, + irrad: TwoLeafIrradience, + leaf_area_index: NDArray, + ): + """Initialize the ``TwoLeafAssimilation`` class. + + Args: + pmodel (PModel | SubdailyPModel): The photosynthesis model. + irrad (TwoLeafIrradience): The irradiance data. + leaf_area_index (NDArray): Array of leaf area index values. + """ + self.pmodel = pmodel + self.irrad = irrad + self.leaf_area_index = leaf_area_index + + self.vcmax_pmod: NDArray + self.vcmax25_pmod: NDArray + self.optchi_obj: OptimalChiABC + self.core_const = self._get_core_const() + self.gpp: NDArray + + self._get_vcmax_pmod() + self._get_vcmax_25_pmod() + self._get_optchi_obj() + + def _get_vcmax_pmod(self) -> None: + """Retrieve the maximum carboxylation rate from the photosynthesis model. + + Sets: + vcmax_pmod (NDArray): Maximum carboxylation rate based on the model type. + """ + self.vcmax_pmod = ( + self.pmodel.vcmax + if isinstance(self.pmodel, PModel) + else self.pmodel.subdaily_vcmax + ) + + def _get_vcmax_25_pmod(self) -> None: + """Retrieve the max carboxylation rate at 25°C from the photosynthesis model. + + Sets: + vcmax25_pmod (NDArray): Maximum carboxylation rate at 25°C based on the + model type. + """ + self.vcmax25_pmod = ( + self.pmodel.vcmax25 + if isinstance(self.pmodel, PModel) + else self.pmodel.subdaily_vcmax25 + ) + + def _get_optchi_obj(self) -> None: + """Retrieve the optimal chi object from the photosynthesis model. + + Sets: + optchi_obj (:class:`~pyrealm.pmodel.optimal_chi.OptimalChiABC'): Optimal chi + object based on the model type. + """ + self.optchi_obj = ( + self.pmodel.optchi + if isinstance(self.pmodel, PModel) + else self.pmodel.optimal_chi + ) + + def _get_core_const(self) -> CoreConst: + """Retrieve the core constants from the photosynthesis model. + + Returns: + CoreConst: The core constants from the + :class:`~pyrealm.pmodel.pmodel.PModel` or + :class:`~pyrealm.pmodel.subdaily.SubdailyPModel`. + """ + if isinstance(self.pmodel, PModel): + return self.pmodel.core_const + else: + return self.pmodel.env.core_const + + def gpp_estimator(self) -> None: + """Estimate the gross primary production (``GPP``) using the two-leaf model. + + This method calculates various parameters related to photosynthesis, including + carboxylation rates, assimilation rates, and electron transport rates for sunlit + and shaded leaves, ultimately leading to the estimation of ``GPP``. + + Sets: + gpp_estimate (NDArray): Estimated gross primary production. + """ + self.kv_Lloyd = canopy_extinction_coefficient(self.vcmax_pmod) + + self.Vmax25_canopy = Vmax25_canopy( + self.leaf_area_index, self.vcmax25_pmod, self.kv_Lloyd + ) + self.Vmax25_sun = Vmax25_sun( + self.leaf_area_index, self.vcmax25_pmod, self.kv_Lloyd, self.irrad.kb + ) + self.Vmax25_shade = Vmax25_shade(self.Vmax25_canopy, self.Vmax25_sun) + + self.Vmax_sun = carboxylation_scaling_to_T(self.Vmax25_sun, self.pmodel.env.tc) + self.Vmax_shade = carboxylation_scaling_to_T( + self.Vmax25_shade, self.pmodel.env.tc + ) + + self.Av_sun = photosynthetic_estimate(self.Vmax_sun, self.optchi_obj.mc) + self.Av_shade = photosynthetic_estimate(self.Vmax_shade, self.optchi_obj.mc) + + self.Jmax25_sun = Jmax25(self.Vmax25_sun) + self.Jmax25_shade = Jmax25(self.Vmax25_shade) + + self.Jmax_sun = Jmax25_temp_correction(self.Jmax25_sun, self.pmodel.env.tc) + self.Jmax_shade = Jmax25_temp_correction(self.Jmax25_shade, self.pmodel.env.tc) + + self.J_sun = electron_transport_rate(self.Jmax_sun, self.irrad.I_csun) + self.J_shade = electron_transport_rate(self.Jmax_shade, self.irrad.I_cshade) + + self.Aj_sun = assimilation_rate(self.optchi_obj.mj, self.J_sun) + self.Aj_shade = assimilation_rate(self.optchi_obj.mj, self.J_shade) + + self.Acanopy_sun = assimilation_canopy( + self.Aj_sun, self.Av_sun, self.irrad.beta_angle, self.irrad.k_sol_obs_angle + ) + self.Acanopy_shade = assimilation_canopy( + self.Aj_shade, + self.Av_shade, + self.irrad.beta_angle, + self.irrad.k_sol_obs_angle, + ) + + self.gpp_estimate = gross_primary_product( + self.core_const.k_c_molmass, self.Acanopy_sun, self.Acanopy_shade + ) + + +def beam_extinction_coeff( + beta_angle: NDArray, + k_sol_obs_angle: float, + clip_angle: float = 30, + kb_numerator: float = 0.5, +) -> NDArray: + r"""Calculate the beam extinction coefficient :math:`kb`. + + The beam extinction coefficient :math:`kb` represents the attenuation of direct + sunlight through the canopy. It is influenced by the solar elevation angle. + + .. math:: + + \text{kb} = + \begin{cases} + \frac{\text{kb\_numerator}}{\sin(\beta\_angle)} & \text{if } \beta\_angle > + \text{k\_sol\_obs\_angle} \\ + \text{clip\_angle} & \text{otherwise} + \end{cases} + + Args: + beta_angle (NDArray): Array of ``solar elevation`` angles. + k_sol_obs_angle (float): Solar angle threshold for calculating the extinction + coefficient. + clip_angle (float, optional): Angle used when solar elevation is below the + threshold. Defaults to 30. + kb_numerator (float, optional): Numerator used in the calculation of the + extinction coefficient. Defaults to 0.5. + + Returns: + NDArray: Array of beam extinction coefficients. + """ + kb = np.where( + beta_angle > k_sol_obs_angle, kb_numerator / np.sin(beta_angle), clip_angle + ) + return kb + + +def scattered_beam_extinction_coeff( + beta_angle: NDArray, k_sol_obs_angle: float +) -> NDArray: + r"""Calculate the scattered beam extinction coefficient :math:`kb_prime`. + + The scattered beam extinction coefficient accounts for the attenuation of + scattered sunlight in the canopy. + + .. math:: + \text{kb\_prime} = \text{beam\_extinction\_coeff}(\beta\_angle, + \text{k\_sol\_obs\_angle}, 27, 0.46) + + Args: + beta_angle (NDArray): Array of solar elevation angles. + k_sol_obs_angle (float): Solar angle threshold for calculating the extinction + coefficient. + + Returns: + NDArray: Array of scattered beam extinction coefficients. + """ + + kb_prime = beam_extinction_coeff(beta_angle, k_sol_obs_angle, 27, 0.46) + + return kb_prime + + +def fraction_of_diffuse_rad( + patm: NDArray, pa0: float, beta_angle: NDArray, k_fa: float +) -> NDArray: + r"""Calculate the fraction of diffuse radiation ``fd``. + + The fraction of diffuse radiation represents the proportion of sunlight that is + scattered before reaching the canopy. + + .. math:: + + m = \frac{\text{patm}}{\text{pa0}} \cdot \frac{1}{\sin(\beta\_angle)} + fd = \frac{1 - 0.72^m}{1 + 0.72^m \cdot \left(\frac{1}{k\_fa} - 1\right)} + + Args: + patm (NDArray): Array of atmospheric pressure values. + pa0 (float): Reference atmospheric pressure. + beta_angle (NDArray): Array of solar elevation angles. + k_fa (float): Constant for calculating the fraction of diffuse radiation. + + Returns: + NDArray: Array of fractions of diffuse radiation. + """ + + m = (patm / pa0) / np.sin(beta_angle) + fd = (1 - 0.72**m) / (1 + (0.72**m * (1 / k_fa - 1))) + + return fd + + +def beam_irradience_h_leaves(k_sigma: float) -> float: + r"""Calculate the beam irradiance :math:`rho_h` for horizontal leaves. + + The beam irradiance for horizontal leaves considers the leaf orientation and + the direct sunlight received. + + .. math:: + + \rho_h = \frac{1 - \sqrt{1 - k\_sigma}}{1 + \sqrt{1 - k\_sigma}} + + Args: + k_sigma (float): Parameter representing the fraction of light intercepted by + horizontal leaves. + + Returns: + NDArray: Array of beam irradiances for horizontal leaves. + """ + + rho_h = (1 - np.sqrt(1 - k_sigma)) / (1 + np.sqrt(1 - k_sigma)) + + return rho_h + + +def beam_irrad_unif_leaf_angle_dist(rho_h: float, kb: NDArray) -> NDArray: + r"""Calculate the beam irradiance with a uniform leaf angle distribution. + + The beam irradiance with a uniform leaf angle distribution :math:`rho_cb` considers + different leaf orientations within the canopy. + + .. math:: + + \rho_{cb} = 1 - \exp\left(-\frac{2 \rho_h kb}{1 + kb}\right) + + Args: + rho_h (NDArray): Array of beam irradiances for horizontal leaves. + kb (NDArray): Array of beam extinction coefficients. + + Returns: + NDArray: Array of beam irradiances with uniform leaf angle distribution. + """ + + rho_cb = 1 - np.exp(-2 * rho_h * kb / (1 + kb)) + + return rho_cb + + +def diffuse_radiation(fd: NDArray, ppfd: NDArray) -> NDArray: + r"""Calculate the diffuse radiation :math`I\d`. + + The diffuse radiation is the portion of sunlight that is scattered in the + atmosphere before reaching the canopy. + + .. math:: + + I_d = \max(0, \text{ppfd} \cdot \text{fd}) + + Args: + fd (NDArray): Array of fractions of diffuse radiation. + ppfd (NDArray): Array of photosynthetic photon flux density values. + + Returns: + NDArray: Array of diffuse radiation values. + """ + + I_d = np.clip(ppfd * fd, a_min=0, a_max=np.inf) + + return I_d + + +def beam_irradience(ppfd: NDArray, fd: NDArray) -> NDArray: + r"""Calculate the beam irradiance :math:`I_b`. + + The beam irradiance is the direct component of sunlight that reaches the canopy + without being scattered. + + .. math:: + + I_b = \text{ppfd} \cdot (1 - \text{fd}) + + Args: + ppfd (NDArray): Array of photosynthetic photon flux density values. + fd (NDArray): Array of fractions of diffuse radiation. + + Returns: + NDArray: Array of beam irradiance values. + """ + + I_b = ppfd * (1 - fd) + + return I_b + + +def scattered_beam_irradience( + I_b: NDArray, + kb: NDArray, + kb_prime: NDArray, + rho_cb: NDArray, + leaf_area_index: NDArray, + k_sigma: NDArray, +) -> NDArray: + r"""Calculate the scattered beam irradiance :math:`I_bs`. + + The scattered beam irradiance is the portion of direct sunlight that is + scattered within the canopy. + + .. math:: + + I_{bs} = I_b \cdot (1 - \rho_{cb}) \cdot kb\_prime \cdot \exp(-kb\_prime \cdot + leaf\_area\_index) - (1 - k\_sigma) \cdot kb \cdot \exp(-kb \cdot + leaf\_area\_index) + + Args: + I_b (NDArray): Array of beam irradiance values. + kb (NDArray): Array of beam extinction coefficients. + kb_prime (NDArray): Array of scattered beam extinction coefficients. + rho_cb (NDArray): Array of beam irradiances with uniform leaf angle + distribution. + leaf_area_index (NDArray): Array of leaf area index values. + k_sigma (NDArray): Array of sigma values. + + Returns: + NDArray: Array of scattered beam irradiance values. + """ + + I_bs = I_b * (1 - rho_cb) * kb_prime * np.exp(-kb_prime * leaf_area_index) - ( + 1 - k_sigma + ) * kb * np.exp(-kb * leaf_area_index) + + return I_bs + + +def canopy_irradience( + rho_cb: NDArray, + I_b: NDArray, + kb_prime: NDArray, + I_d: NDArray, + leaf_area_index: NDArray, + k_rho_cd: float, +) -> NDArray: + r"""Calculate the canopy irradiance :math:`I_c`. + + The canopy irradiance is the total irradiance within the canopy, including both + direct and diffuse radiation components. + + .. math:: + + I_c = (1 - \rho_{cb}) \cdot I_b \cdot + (1 - \exp(-kb\_prime \cdot leaf\_area\_index)) + + (1 - k\_rho\_cd) \cdot I_d \cdot (1 - \exp(-kb\_prime \cdot leaf\_area\_index)) + + Args: + rho_cb (NDArray): Array of beam irradiances with uniform leaf angle + distribution. + I_b (NDArray): Array of beam irradiance values. + kb_prime (NDArray): Array of scattered beam extinction coefficients. + I_d (NDArray): Array of diffuse radiation values. + leaf_area_index (NDArray): Array of leaf area index values. + k_rho_cd (float): Constant for calculating the diffuse radiation component. + + Returns: + NDArray: Array of canopy irradiance values. + """ + I_c = (1 - rho_cb) * I_b * (1 - np.exp(-kb_prime * leaf_area_index)) + ( + 1 - k_rho_cd + ) * I_d * (1 - np.exp(-kb_prime * leaf_area_index)) + + return I_c + + +def sunlit_beam_irrad( + I_b: NDArray, k_sigma: float, kb: NDArray, leaf_area_index: NDArray +) -> NDArray: + r"""Calculate the sunlit beam irradiance :math:`Isun_beam`. + + The sunlit beam irradiance is the direct sunlight received by the sunlit portion + of the canopy. + + .. math:: + + I_{sun\_beam} = I_b \cdot (1 - k\_sigma) \cdot (1 - \exp(-kb \cdot + leaf\_area\_index)) + + Args: + I_b (NDArray): Array of beam irradiance values. + k_sigma (float): Constant for calculating the sunlit beam irradiance. + kb (NDArray): Array of beam extinction coefficients. + leaf_area_index (NDArray): Array of leaf area index values. + + Returns: + NDArray: Array of sunlit beam irradiance values. + """ + Isun_beam = I_b * (1 - k_sigma) * (1 - np.exp(-kb * leaf_area_index)) + + return Isun_beam + + +def sunlit_diffuse_irrad( + I_d: NDArray, + k_rho_cd: float, + k_kd_prime: float, + kb: NDArray, + leaf_area_index: NDArray, +) -> NDArray: + r"""Calculate the sunlit diffuse irradiance :math:`Isun_diffuse`. + + The sunlit diffuse irradiance is the diffuse radiation received by the sunlit + portion of the canopy. + + .. math:: + + I_{sun\_diffuse} = I_d \cdot (1 - k\_rho\_cd) \cdot + (1 - \exp(-(k\_kd\_prime + kb) \cdot + leaf\_area\_index)) \cdot \frac{k\_kd\_prime}{k\_kd\_prime + kb} + + Args: + I_d (NDArray): Array of diffuse radiation values. + k_rho_cd (float): Constant rho_cd value. + k_kd_prime (float): Constant for calculating the sunlit diffuse irradiance. + kb (NDArray): Array of beam extinction coefficients. + leaf_area_index (NDArray): Array of leaf area index values. + + Returns: + NDArray: Array of sunlit diffuse irradiance values. + """ + Isun_diffuse = ( + I_d + * (1 - k_rho_cd) + * (1 - np.exp(-(k_kd_prime + kb) * leaf_area_index)) + * k_kd_prime + / (k_kd_prime + kb) + ) + + return Isun_diffuse + + +def sunlit_scattered_irrad( + I_b: NDArray, + rho_cb: NDArray, + kb_prime: NDArray, + kb: NDArray, + leaf_area_index: NDArray, + k_sigma: float, +) -> NDArray: + r"""Calculate the sunlit scattered irradiance :math:`Isun_scattered`. + + The sunlit scattered irradiance is the scattered radiation received by the sunlit + portion of the canopy. + + .. math:: + + I_{sun\_scattered} = I_b \cdot ((1 - \rho_{cb}) \cdot + (1 - \exp(-(kb\_prime + kb) \cdot + leaf\_area\_index)) \cdot \frac{kb\_prime}{kb\_prime + kb} - (1 - k\_sigma) + \cdot (1 - \exp(-2 \cdot kb \cdot leaf\_area\_index)) / 2) + + Args: + I_b (NDArray): Array of beam irradiance values. + rho_cb (NDArray): Array of beam irradiances with uniform leaf angle + distribution. + kb_prime (NDArray): Array of scattered beam extinction coefficients. + kb (NDArray): Array of beam extinction coefficients. + leaf_area_index (NDArray): Array of leaf area index values. + k_sigma (float): Constant for calculating the sunlit scattered irradiance. + + Returns: + NDArray: Array of sunlit scattered irradiance values. + """ + Isun_scattered = I_b * ( + (1 - rho_cb) + * (1 - np.exp(-(kb_prime + kb) * leaf_area_index)) + * kb_prime + / (kb_prime + kb) + - (1 - k_sigma) * (1 - np.exp(-2 * kb * leaf_area_index)) / 2 + ) + + return Isun_scattered + + +def sunlit_absorbed_irrad( + Isun_beam: NDArray, Isun_diffuse: NDArray, Isun_scattered: NDArray +) -> NDArray: + r"""Calculate the sunlit absorbed irradiance :math:`I_csun`. + + The sunlit absorbed irradiance is the total irradiance absorbed by the sunlit + portion of the canopy, combining beam, diffuse, and scattered irradiance. + + .. math:: + + I_{csun} = I_{sun\_beam} + I_{sun\_diffuse} + I_{sun\_scattered} + + Args: + Isun_beam (NDArray): Array of sunlit beam irradiance values. + Isun_diffuse (NDArray): Array of sunlit diffuse irradiance values. + Isun_scattered (NDArray): Array of sunlit scattered irradiance values. + + Returns: + NDArray: Array of sunlit absorbed irradiance values. + """ + I_csun = Isun_beam + Isun_scattered + Isun_diffuse + + return I_csun + + +def shaded_absorbed_irrad( + beta_angle: NDArray, k_sol_obs_angle: float, I_c: NDArray, I_csun: NDArray +) -> NDArray: + r"""Calculate the irradiance absorbed by the shaded fraction of the canopy. + + The irradiance absorbed by the shaded fraction of the canopy :math:`I_cshade`is + calculated by subtracting the sunlit absorbed irradiance from the total canopy + irradiance. + + .. math:: + + I_{cshade} = \max(0, I_c - I_{csun}) + + Args: + beta_angle (NDArray): Array of solar elevation angles. + k_sol_obs_angle (float): Solar angle threshold. + I_c (NDArray): Array of canopy irradiance values. + I_csun (NDArray): Array of sunlit absorbed irradiance values. + + Returns: + NDArray: Array of irradiance absorbed by the shaded fraction of the canopy. + """ + I_cshade = np.where(beta_angle > k_sol_obs_angle, I_c - I_csun, 0) + + return I_cshade + + +def canopy_extinction_coefficient(vcmax_pmod: NDArray) -> NDArray: + r"""Calculate :math:`k_v` parameter. + + This function calculates the extinction coefficient, :math:`k_v`, which represents + how the photosynthetic capacity (Vcmax) decreases with depth in the plant canopy. + The exponential model used here is derived from empirical data and represents how + light attenuation affects photosynthetic capacity vertically within the canopy. + + Equation is sourced from Figure 10 of Lloyd et al. (2010). + + .. math:: + \text{kv\_Lloyd} = \exp(0.00963 \cdot \text{vcmax\_pmod} - 2.43) + + Note: + ``Vcmax`` is used here rather than ``Vcmax_25``. + + Args: + vcmax_pmod (NDArray): The ``Vcmax`` parameter for the pmodel. + + Returns: + NDArray: The calculated :math:`kv_Lloyd` values. + """ + kv_Lloyd = np.exp(0.00963 * vcmax_pmod - 2.43) + return kv_Lloyd + + +def Vmax25_canopy( + leaf_area_index: NDArray, vcmax25_pmod: NDArray, kv: NDArray +) -> NDArray: + r"""Calculate carboxylation rate, :math:`V_max25` in the canopy at 25C. + + This function calculates the maximum carboxylation rate of the canopy at a + reference temperature of 25°C. It integrates the vertical gradient of + photosynthetic capacity across the leaf area index (``LAI``), considering the + extinction coefficient for light. + + .. math:: + \text{Vmax25\_canopy} = \text{LAI} \cdot \text{vcmax25\_pmod} \cdot + \left(\frac{1 - \exp(-\text{kv})}{\text{kv}}\right) + + + Args: + leaf_area_index (NDArray): The ``leaf area index``. + vcmax25_pmod (NDArray): The ``Vcmax25`` parameter for the pmodel. + kv (NDArray): The kv parameter. + + Returns: + NDArray: The calculated Vmax25 canopy values. + """ + Vmax25_canopy = leaf_area_index * vcmax25_pmod * ((1 - np.exp(-kv)) / kv) + + return Vmax25_canopy + + +def Vmax25_sun( + leaf_area_index: NDArray, vcmax25_pmod: NDArray, kv: NDArray, kb: NDArray +) -> NDArray: + r"""Calculate carboxylation in sunlit areas, :math:`Vmax25_sun` at 25C. + + This function calculates the maximum carboxylation rate for the sunlit portions of + the canopy at 25°C. It considers both the extinction coefficient and the direct + sunlight penetration parameter :math:`k_b`. + + .. math:: + + \text{Vmax25\_sun} = \text{leaf\_area\_index} \cdot \text{vcmax25\_pmod} \cdot + \left( \frac{1 - \exp(-\text{kv} - \text{kb} \cdot \text{leaf\_area\_index})} + {\text{kv} + \text{kb} \cdot \text{leaf\_area\_index}} \right) + + Args: + leaf_area_index (NDArray): The leaf area index. + vcmax25_pmod (NDArray): The Vcmax25 parameter for the pmodel. + kv (NDArray): The kv parameter. + kb (NDArray): The irradiation kb parameter. + + Returns: + NDArray: The calculated Vmax25 sun values. + """ + Vmax25_sun = ( + leaf_area_index + * vcmax25_pmod + * ((1 - np.exp(-kv - kb * leaf_area_index)) / (kv + kb * leaf_area_index)) + ) + + return Vmax25_sun + + +def Vmax25_shade(Vmax25_canopy: NDArray, Vmax_25_sun: NDArray) -> NDArray: + r"""Calculate carboxylation in shaded areas, :math:`Vmax25_shade` at 25C. + + This function calculates the maximum carboxylation rate for the shaded portions of + the canopy by subtracting the sunlit carboxylation rate from the total canopy + carboxylation rate. + + .. math:: + + \text{Vmax25\_shade} = \text{Vmax25\_canopy} - \text{Vmax25\_sun} + + Args: + Vmax25_canopy (NDArray): The ``Vmax25`` parameter for the canopy. + Vmax_25_sun (NDArray): The ``Vmax25`` parameter for the sunlit areas. + + Returns: + NDArray: The calculated Vmax25 shade values. + """ + Vmax25_shade = Vmax25_canopy - Vmax_25_sun + + return Vmax25_shade + + +def carboxylation_scaling_to_T(Vmax25: NDArray, tc: NDArray) -> NDArray: + r"""Convert carboxylation rates at 25C to rate at ambient temperature (C). + + This function adjusts the carboxylation rate from the reference temperature of 25°C + to the ambient temperature using an Arrhenius-type function, which describes the + temperature dependence of enzymatic reactions. + + .. math:: + \text{Vmax\_sun} = \text{Vmax25} \cdot \exp \left(\frac{64800 \cdot (\text{tc} + - 25)}{298 \cdot 8.314 \cdot (\text{tc} + 273)}\right) + + Args: + Vmax25 (NDArray): The ``Vmax25`` parameter. + tc (NDArray): The temperature in Celsius. + + Returns: + NDArray: The carboxylation rates adjusted to ambient temperature. + """ + + Vmax = Vmax25 * np.exp(64800 * (tc - 25) / (298 * 8.314 * (tc + 273))) + + return Vmax + + +def photosynthetic_estimate(Vmax: NDArray, mc: NDArray) -> NDArray: + r"""Calculate photosynthetic rate estimate, :math:`Av`. + + This function estimates photosynthetic rates by multiplying the carboxylation + capacity by mc, the limitation term for Rubisco-limited assimilation. + + .. math:: + + A_v = V_{max} \cdot m_c + + Args: + Vmax (NDArray): The ``Vmax`` parameter. + mc (NDArray): limitation term for Rubisco-limited assimilation + + Returns: + NDArray: The calculated photosynthetic estimates. + """ + + Av = Vmax * mc + + return Av + + +def Jmax25(Vmax25: NDArray) -> NDArray: + r"""Calculate the maximum rate of electron transport :math:`Jmax25`. + + This function calculates the maximum rate of electron transport :math:`Jmax25` + at 25°C, which is related to the capacity for light-driven electron transport in + photosynthesis. + + Uses Eqn 31, after Wullschleger. + + .. math:: + + J_{max25} = 29.1 + 1.64 \cdot V_{max25} + + + Args: + Vmax25 (NDArray): The ``Vmax25`` parameter. + + Returns: + NDArray: The calculated ``Jmax25`` values. + """ + + Jmax25 = 29.1 + 1.64 * Vmax25 + + return Jmax25 + + +def Jmax25_temp_correction(Jmax25: NDArray, tc: NDArray) -> NDArray: + r"""Corrects Jmax value to ambient temperature. + + This function adjusts the maximum electron transport rate ``Jmax25`` for temperature + using a temperature correction formula, similar to the Arrhenius equation. + + Correction derived from Mengoli 2021 Eqn 3b. + + .. math:: + + J_{max} = J_{max25} \cdot \exp\left(\frac{43990}{8.314} + \left( \frac{1}{298} - \frac{1}{\text{tc} + 273} \right)\right) + + + Args: + Jmax25 (NDArray): The ``Jmax25`` parameter. + tc (NDArray): The temperature in Celsius. + + Returns: + NDArray: The temperature-corrected Jmax values. + """ + + Jmax = Jmax25 * np.exp((43990 / 8.314) * (1 / 298 - 1 / (tc + 273))) + + return Jmax + + +def electron_transport_rate(Jmax: NDArray, I_c: NDArray) -> NDArray: + r"""Calculate electron transport rate :math:`J`. + + This function calculates the electron transport rate :math:`J`,considering the + irradiance :math:`I_c` and the maximum electron transport rate :math:`Jmax`. + + .. math:: + + J = J_{max} \cdot I_c \cdot \frac{(1 - 0.15)}{(I_c + 2.2 \cdot J_{max})} + + Args: + Jmax (NDArray): maximum rate of electron transport. + I_c (NDArray): The irradiance parameter. + + Returns: + NDArray: The calculated J values. + """ + + J = Jmax * I_c * (1 - 0.15) / (I_c + 2.2 * Jmax) + + return J + + +def assimilation_rate(mj: NDArray, J: NDArray) -> NDArray: + r"""Calculate assimilation rate :math:`A`. + + This function calculates the assimilation rate driven by electron transport, + :math:`Aj`, using the parameter :math:`mj` and the electron transport rate + :math:`J`. + + .. math:: + + A = m_j \cdot \frac{J}{4} + + Args: + mj (NDArray): The ``mj`` parameter. + J (NDArray): The ``J`` parameter. + + Returns: + NDArray: The calculated Aj values. + """ + A = mj * J / 4 + + return A + + +def assimilation_canopy( + Aj: NDArray, Av: NDArray, beta_angle: NDArray, solar_obscurity_angle: float +) -> NDArray: + r"""Calculate assimilation in canopy, :math:`Acanopy`. + + This function calculates the total canopy assimilation by taking the minimum of + :math:`Aj` and :math:`Av` for each partition and clipping the data when the sun is + below the solar obscurity angle. + + .. math:: + + \text{ew\_minima} = \min(A_j, A_v) \\ + A_{\text{canopy}} = \begin{cases} + 0 & \text{if } \beta_{\text{angle}} < \text{solar\_obscurity\_angle} \\ + \text{ew\_minima} & \text{otherwise} + \end{cases} + + Args: + Aj (NDArray): The ``Aj`` parameter. + Av (NDArray): The ``Av`` parameter. + beta_angle (NDArray): The solar elevation :math:`beta` angle. + solar_obscurity_angle (float): The solar obscurity angle. + + Returns: + NDArray: The calculated assimilation values for the canopy. + """ + + ew_minima = np.minimum(Aj, Av) + Acanopy = np.where(beta_angle < solar_obscurity_angle, 0, ew_minima) + + return Acanopy + + +def gross_primary_product( + k_c_molmass: float, Acanopy_sun: NDArray, Acanopy_shade: NDArray +) -> NDArray: + r"""Calculate gross primary productivity (``GPP``). + + This function calculates the ``GPP`` by combining the assimilation rates of sunlit + (:math:`Acanopy_sun`) and shaded (:math:`Acanopy_shade`) portions of the canopy, + scaled by a molar mass constant, :math:`k_c_molmass`. + + .. math:: + + \text{gpp} = k_{\text{c\_molmass}} \cdot A_{\text{canopy\_sun}} + + A_{\text{canopy\_shade}} + + Args: + k_c_molmass (float): The constant for molar mass. + Acanopy_sun (NDArray): The assimilation values for the sunlit canopy. + Acanopy_shade (NDArray): The assimilation values for the shaded canopy. + + Returns: + NDArray: The calculated GPP values. + """ + + gpp = k_c_molmass * Acanopy_sun + Acanopy_shade + + return gpp diff --git a/pyrealm_build_data/two_leaf/Untitled.ipynb b/pyrealm_build_data/two_leaf/Untitled.ipynb new file mode 100644 index 00000000..e8487fdc --- /dev/null +++ b/pyrealm_build_data/two_leaf/Untitled.ipynb @@ -0,0 +1,828 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "id": "214800c6-1e69-4d26-af47-7d227b62397c", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "from pvlib.location import Location\n", + "\n", + "from pyrealm.pmodel import PModel, PModelEnvironment, SubdailyPModel" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "b7f0cdec-7722-40a5-a7dd-d2e4aad2cd42", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
timeppfdvpdtcco2patmgpp_fluxnetLAIfaparGPP_JAMES
02014-01-01 00:00:002.052.23.33418.40795300.00.01.9426980.523071NaN
12014-01-01 00:30:002.00.02.71418.83695300.00.01.9426980.523071NaN
22014-01-01 01:00:002.00.02.91418.72995300.00.01.9426980.523071NaN
32014-01-01 01:30:002.00.03.12418.47295300.00.01.9426980.523071NaN
42014-01-01 02:00:002.00.03.21418.39595300.00.01.9426980.523071NaN
.................................
175152014-12-31 21:30:000.00.0-1.24394.19397200.00.01.9538210.444753NaN
175162014-12-31 22:00:000.00.0-1.46396.10997200.00.01.9538210.444753NaN
175172014-12-31 22:30:000.00.0-1.70395.77697200.00.01.9538210.444753NaN
175182014-12-31 23:00:000.00.0-1.93394.60197200.00.01.9538210.444753NaN
175192014-12-31 23:30:000.00.0-2.16393.45897300.00.01.9538210.444753NaN
\n", + "

17520 rows × 10 columns

\n", + "
" + ], + "text/plain": [ + " time ppfd vpd tc co2 patm gpp_fluxnet \\\n", + "0 2014-01-01 00:00:00 2.0 52.2 3.33 418.407 95300.0 0.0 \n", + "1 2014-01-01 00:30:00 2.0 0.0 2.71 418.836 95300.0 0.0 \n", + "2 2014-01-01 01:00:00 2.0 0.0 2.91 418.729 95300.0 0.0 \n", + "3 2014-01-01 01:30:00 2.0 0.0 3.12 418.472 95300.0 0.0 \n", + "4 2014-01-01 02:00:00 2.0 0.0 3.21 418.395 95300.0 0.0 \n", + "... ... ... ... ... ... ... ... \n", + "17515 2014-12-31 21:30:00 0.0 0.0 -1.24 394.193 97200.0 0.0 \n", + "17516 2014-12-31 22:00:00 0.0 0.0 -1.46 396.109 97200.0 0.0 \n", + "17517 2014-12-31 22:30:00 0.0 0.0 -1.70 395.776 97200.0 0.0 \n", + "17518 2014-12-31 23:00:00 0.0 0.0 -1.93 394.601 97200.0 0.0 \n", + "17519 2014-12-31 23:30:00 0.0 0.0 -2.16 393.458 97300.0 0.0 \n", + "\n", + " LAI fapar GPP_JAMES \n", + "0 1.942698 0.523071 NaN \n", + "1 1.942698 0.523071 NaN \n", + "2 1.942698 0.523071 NaN \n", + "3 1.942698 0.523071 NaN \n", + "4 1.942698 0.523071 NaN \n", + "... ... ... ... \n", + "17515 1.953821 0.444753 NaN \n", + "17516 1.953821 0.444753 NaN \n", + "17517 1.953821 0.444753 NaN \n", + "17518 1.953821 0.444753 NaN \n", + "17519 1.953821 0.444753 NaN \n", + "\n", + "[17520 rows x 10 columns]" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "forcing_data = pd.read_csv(\"merged_BE-Vie_data.csv\")\n", + "forcing_data[\"time\"] = pd.to_datetime(forcing_data[\"time\"])\n", + "forcing_data" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "1952b1fd-e3ed-424a-bf73-b2eb427e7e53", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
apparent_zenithzenithapparent_elevationelevationazimuthequation_of_time
time
2014-01-01 00:00:00152.420122152.420122-62.420122-62.42012210.321119-3.304259
2014-01-01 00:30:00150.978178150.978178-60.978178-60.97817824.585260-3.314157
2014-01-01 01:00:00148.506258148.506258-58.506258-58.50625837.398486-3.324052
2014-01-01 01:30:00145.237642145.237642-55.237642-55.23764248.541226-3.333946
2014-01-01 02:00:00141.393309141.393309-51.393309-51.39330958.150628-3.343837
.....................
2014-12-31 21:30:00142.950364142.950364-52.950364-52.950364305.337753-3.135913
2014-12-31 22:00:00146.596355146.596355-56.596355-56.596355315.510680-3.145789
2014-12-31 22:30:00149.585960149.585960-59.585960-59.585960327.298609-3.155662
2014-12-31 23:00:00151.687297151.687297-61.687297-61.687297340.729300-3.165533
2014-12-31 23:30:00152.679438152.679438-62.679438-62.679438355.391879-3.175401
\n", + "

17520 rows × 6 columns

\n", + "
" + ], + "text/plain": [ + " apparent_zenith zenith apparent_elevation \\\n", + "time \n", + "2014-01-01 00:00:00 152.420122 152.420122 -62.420122 \n", + "2014-01-01 00:30:00 150.978178 150.978178 -60.978178 \n", + "2014-01-01 01:00:00 148.506258 148.506258 -58.506258 \n", + "2014-01-01 01:30:00 145.237642 145.237642 -55.237642 \n", + "2014-01-01 02:00:00 141.393309 141.393309 -51.393309 \n", + "... ... ... ... \n", + "2014-12-31 21:30:00 142.950364 142.950364 -52.950364 \n", + "2014-12-31 22:00:00 146.596355 146.596355 -56.596355 \n", + "2014-12-31 22:30:00 149.585960 149.585960 -59.585960 \n", + "2014-12-31 23:00:00 151.687297 151.687297 -61.687297 \n", + "2014-12-31 23:30:00 152.679438 152.679438 -62.679438 \n", + "\n", + " elevation azimuth equation_of_time \n", + "time \n", + "2014-01-01 00:00:00 -62.420122 10.321119 -3.304259 \n", + "2014-01-01 00:30:00 -60.978178 24.585260 -3.314157 \n", + "2014-01-01 01:00:00 -58.506258 37.398486 -3.324052 \n", + "2014-01-01 01:30:00 -55.237642 48.541226 -3.333946 \n", + "2014-01-01 02:00:00 -51.393309 58.150628 -3.343837 \n", + "... ... ... ... \n", + "2014-12-31 21:30:00 -52.950364 305.337753 -3.135913 \n", + "2014-12-31 22:00:00 -56.596355 315.510680 -3.145789 \n", + "2014-12-31 22:30:00 -59.585960 327.298609 -3.155662 \n", + "2014-12-31 23:00:00 -61.687297 340.729300 -3.165533 \n", + "2014-12-31 23:30:00 -62.679438 355.391879 -3.175401 \n", + "\n", + "[17520 rows x 6 columns]" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "pv_location = Location(\n", + " latitude=50.30493, longitude=5.99812, tz=\"Etc/GMT+1\", altitude=493, name=\"BE-Vie\"\n", + ")\n", + "solar_position = pv_location.get_solarposition(forcing_data[\"time\"])\n", + "solar_position" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "ce5d67d1-0fb5-4325-a06c-d34444597e59", + "metadata": {}, + "outputs": [], + "source": [ + "# Add solar elevation as radians to the forcing dataframe\n", + "forcing_data[\"solar_zenith\"] = np.deg2rad(solar_position[\"zenith\"].to_numpy())\n", + "forcing_data[\"solar_elevation\"] = np.deg2rad(solar_position[\"elevation\"].to_numpy())\n", + "data = forcing_data.loc[(forcing_data[\"time\"] == \"2014-08-01 12:30:00\")]" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "381e43b7-83af-413e-b861-a0f8f8a82243", + "metadata": {}, + "outputs": [], + "source": [ + "class TwoLeafIrradiance:\n", + " def __init__(\n", + " self,\n", + " beta_angle,\n", + " PPFD,\n", + " LAI,\n", + " PATM,\n", + " solar_obscurity_angle=np.deg2rad(1),\n", + " PA0=101325,\n", + " fa=0.426,\n", + " sigma=0.15,\n", + " rho_cd=0.036,\n", + " kd_prime=0.719,\n", + " ):\n", + " \"\"\"Partition solar irradiance\n", + "\n", + " The references in parentheses e.g. (A23) refer to the equation sequence in\n", + " deP&F 1997 and recreated here.\n", + " \"\"\"\n", + "\n", + " # Store key inputs\n", + " self.solar_obscurity_angle = solar_obscurity_angle\n", + " self.beta_angle = beta_angle\n", + "\n", + " # beam extinction coefficient, clipped below solar_obscurity_angle\n", + " self.kb = np.where(\n", + " beta_angle > solar_obscurity_angle, 0.5 / np.sin(self.beta_angle), 30\n", + " )\n", + "\n", + " # beam and scattered-beam extinction coefficient\n", + " kb_prime = np.where(\n", + " beta_angle > solar_obscurity_angle, 0.46 / np.sin(self.beta_angle), 27\n", + " )\n", + "\n", + " m = (PATM / PA0) / np.sin(beta_angle)\n", + "\n", + " fd = (1 - 0.72**m) / (1 + ((0.72**m) * ((1 / fa) - 1)))\n", + "\n", + " # beam irradiance horizontal leaves (A20)\n", + " rho_h = (1 - np.sqrt(1 - sigma)) / (1 + np.sqrt(1 - sigma))\n", + "\n", + " # beam irradiance uniform leaf-angle distribution (A19)\n", + " rho_cb = 1 - np.exp(-2 * rho_h * self.kb / (1 + self.kb))\n", + "\n", + " # diffuse radiation, truncating at zero\n", + " I_d = np.clip(PPFD * fd, a_min=0, a_max=np.inf)\n", + "\n", + " # beam irradiance\n", + " I_b = PPFD * (1 - fd)\n", + "\n", + " # scattered beam irradiance (A8)\n", + " I_bs = I_b * (1 - rho_cb) * kb_prime * np.exp(-kb_prime * LAI) - (\n", + " 1 - sigma\n", + " ) * self.kb * np.exp(-self.kb * LAI)\n", + "\n", + " # Irradiance absorbed by the whole canopy (Eqn 13)\n", + " self.I_c = (1 - rho_cb) * I_b * (1 - np.exp(-kb_prime * LAI)) + (\n", + " 1 - rho_cd\n", + " ) * I_d * (1 - np.exp(-kb_prime * LAI))\n", + "\n", + " # for the sunlit leaves, the components are (Eqn 20):\n", + "\n", + " self.Isun_beam = I_b * (1 - sigma) * (1 - np.exp(-self.kb * LAI))\n", + " self.Isun_diffuse = (\n", + " I_d\n", + " * (1 - rho_cd)\n", + " * (1 - np.exp(-(kd_prime + self.kb) * LAI))\n", + " * kd_prime\n", + " / (kd_prime + self.kb)\n", + " )\n", + " self.Isun_scattered = I_b * (\n", + " (\n", + " (1 - rho_cb)\n", + " * (1 - np.exp(-(kb_prime + self.kb) * LAI))\n", + " * kb_prime\n", + " / (kb_prime + self.kb)\n", + " )\n", + " - ((1 - sigma) * (1 - np.exp(-2 * self.kb * LAI)) / 2)\n", + " )\n", + "\n", + " # Irradiance absorbed by the sunlit fraction of the canopy (Eqn 20):\n", + " self.I_csun = self.Isun_beam + self.Isun_diffuse + self.Isun_scattered\n", + "\n", + " # Irradiance absorbed by the shaded fraction of the canopy (Eqn 21):\n", + " # and including a clause to exclude the hours of obscurity\n", + " self.I_cshade = np.where(\n", + " beta_angle > solar_obscurity_angle, self.I_c - self.I_csun, 0\n", + " )\n", + "\n", + " self.rho_h = rho_h\n", + " self.rho_cb = rho_cb\n", + " self.I_d = I_d\n", + " self.I_b = I_b\n", + " self.I_bs = I_bs\n", + " self.kb_prime = kb_prime" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "bd01339c-0ff5-468a-8a6f-4806abf0d6a2", + "metadata": {}, + "outputs": [], + "source": [ + "two_leaf_irrad = TwoLeafIrradiance(\n", + " beta_angle=data[\"solar_elevation\"].to_numpy(),\n", + " PPFD=data[\"ppfd\"].to_numpy(),\n", + " LAI=data[\"LAI\"].to_numpy(),\n", + " PATM=data[\"patm\"].to_numpy(),\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "e92e3b61-94db-406a-b7c3-54590f7ca902", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "kb = [0.60119491]\n", + "kb_prime = [0.55309931]\n", + "rho_h = 0.04060739027615024\n", + "rho_cb = [0.03003319]\n", + "I_d = [136.29450038]\n", + "I_b = [714.70549962]\n", + "I_bs = [87.564122]\n", + "I_c = [636.08707023]\n", + "Isun_beam = [485.32862357]\n", + "Isun_diffuse = [69.44255865]\n", + "Isun_scattered = [25.43894747]\n", + "beta_angle = [0.98212117]\n" + ] + } + ], + "source": [ + "print(f\"kb = {two_leaf_irrad.kb}\")\n", + "print(f\"kb_prime = {two_leaf_irrad.kb_prime}\")\n", + "print(f\"rho_h = {two_leaf_irrad.rho_h}\")\n", + "print(f\"rho_cb = {two_leaf_irrad.rho_cb}\")\n", + "print(f\"I_d = {two_leaf_irrad.I_d}\")\n", + "print(f\"I_b = {two_leaf_irrad.I_b}\")\n", + "print(f\"I_bs = {two_leaf_irrad.I_bs}\")\n", + "print(f\"I_c = {two_leaf_irrad.I_c}\")\n", + "print(f\"Isun_beam = {two_leaf_irrad.Isun_beam}\")\n", + "print(f\"Isun_diffuse = {two_leaf_irrad.Isun_diffuse}\")\n", + "print(f\"Isun_scattered = {two_leaf_irrad.Isun_scattered}\")\n", + "print(f\"beta_angle = {two_leaf_irrad.beta_angle}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "8cf793c0-9c2a-4bd4-a655-1c437de434c7", + "metadata": {}, + "outputs": [], + "source": [ + "pmod_env = PModelEnvironment(\n", + " tc=data[\"tc\"].to_numpy(),\n", + " patm=data[\"patm\"].to_numpy(),\n", + " co2=data[\"co2\"].to_numpy(),\n", + " vpd=data[\"vpd\"].to_numpy(),\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "id": "a49030e8-bc5b-42ad-b344-6473665556c7", + "metadata": {}, + "outputs": [], + "source": [ + "standard_pmod = PModel(pmod_env, kphio=1 / 8)\n", + "standard_pmod.estimate_productivity(\n", + " fapar=data[\"fapar\"].to_numpy(), ppfd=data[\"ppfd\"].to_numpy()\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "id": "cde288c0-a6bd-43b4-923a-4a9a113a3ee0", + "metadata": {}, + "outputs": [], + "source": [ + "class TwoLeafAssimilation:\n", + " def __init__(\n", + " self, pmodel: PModel | SubdailyPModel, irrad: TwoLeafIrradiance, LAI: np.array\n", + " ):\n", + " # A load of needless inconsistencies in the object attribute names - to be resolved!\n", + " vcmax_pmod = (\n", + " pmodel.vcmax if isinstance(pmodel, PModel) else pmodel.subdaily_vcmax\n", + " )\n", + " vcmax25_pmod = (\n", + " pmodel.vcmax25 if isinstance(pmodel, PModel) else pmodel.subdaily_vcmax25\n", + " )\n", + " optchi_obj = pmodel.optchi if isinstance(pmodel, PModel) else pmodel.optimal_chi\n", + " core_const = (\n", + " pmodel.core_const if isinstance(pmodel, PModel) else pmodel.env.core_const\n", + " )\n", + "\n", + " # Generate extinction coefficients to express the vertical gradient in photosynthetic\n", + " # capacity after the equation provided in Figure 10 of Lloyd et al. (2010):\n", + "\n", + " kv_Lloyd = np.exp(\n", + " 0.00963 * vcmax_pmod - 2.43\n", + " ) # KB: Here I opt for vcmax rather than Vcmax25\n", + "\n", + " # Calculate carboxylation in the two partitions at standard conditions\n", + " self.Vmax25_canopy = LAI * vcmax25_pmod * ((1 - np.exp(-kv_Lloyd)) / kv_Lloyd)\n", + " self.Vmax25_sun = (\n", + " LAI\n", + " * vcmax25_pmod\n", + " * ((1 - np.exp(-kv_Lloyd - irrad.kb * LAI)) / (kv_Lloyd + irrad.kb * LAI))\n", + " )\n", + " self.Vmax25_shade = self.Vmax25_canopy - self.Vmax25_sun\n", + "\n", + " # Convert carboxylation rates to ambient temperature using an Arrhenius function\n", + " self.Vmax_sun = self.Vmax25_sun * np.exp(\n", + " 64800 * (pmodel.env.tc - 25) / (298 * 8.314 * (pmodel.env.tc + 273))\n", + " )\n", + " self.Vmax_shade = self.Vmax25_shade * np.exp(\n", + " 64800 * (pmodel.env.tc - 25) / (298 * 8.314 * (pmodel.env.tc + 273))\n", + " )\n", + "\n", + " # Now the photosynthetic estimates as V_cmax * mc\n", + " self.Av_sun = self.Vmax_sun * optchi_obj.mc\n", + " self.Av_shade = self.Vmax_shade * optchi_obj.mc\n", + "\n", + " ## Jmax estimates for sun and shade;\n", + " self.Jmax25_sun = 29.1 + 1.64 * self.Vmax25_sun # Eqn 31, after Wullschleger\n", + " self.Jmax25_shade = 29.1 + 1.64 * self.Vmax25_shade\n", + "\n", + " # Temperature correction (Mengoli 2021 Eqn 3b); relevant temperatures given in Kelvin\n", + " self.Jmax_sun = self.Jmax25_sun * np.exp(\n", + " (43990 / 8.314) * (1 / 298 - 1 / (pmodel.env.tc + 273))\n", + " )\n", + " self.Jmax_shade = self.Jmax25_shade * np.exp(\n", + " (43990 / 8.314) * (1 / 298 - 1 / (pmodel.env.tc + 273))\n", + " )\n", + "\n", + " # Calculate J and Aj for each partition\n", + " self.J_sun = (\n", + " self.Jmax_sun\n", + " * irrad.I_csun\n", + " * (1 - 0.15)\n", + " / (irrad.I_csun + 2.2 * self.Jmax_sun)\n", + " )\n", + " self.J_shade = (\n", + " self.Jmax_shade\n", + " * irrad.I_cshade\n", + " * (1 - 0.15)\n", + " / (irrad.I_cshade + 2.2 * self.Jmax_shade)\n", + " )\n", + "\n", + " self.Aj_sun = (self.J_sun / 4) * optchi_obj.mj\n", + " self.Aj_shade = (self.J_shade / 4) * optchi_obj.mj\n", + "\n", + " # Calculate the assimilation in each partition as the minimum of Aj and Av,\n", + " # in both cases clipping when the sun is below the angle of obscurity\n", + " Acanopy_sun = np.minimum(self.Aj_sun, self.Av_sun)\n", + " self.Acanopy_sun = np.where(\n", + " irrad.beta_angle < irrad.solar_obscurity_angle, 0, Acanopy_sun\n", + " )\n", + " Acanopy_shade = np.minimum(self.Aj_shade, self.Av_shade)\n", + " self.Acanopy_shade = np.where(\n", + " irrad.beta_angle < irrad.solar_obscurity_angle, 0, Acanopy_shade\n", + " )\n", + "\n", + " self.gpp = core_const.k_c_molmass * self.Acanopy_sun + self.Acanopy_shade" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "id": "205a40de-b49f-43ef-a32a-02319fdeb6b5", + "metadata": {}, + "outputs": [], + "source": [ + "two_leaf_assim_standard = TwoLeafAssimilation(\n", + " irrad=two_leaf_irrad, pmodel=standard_pmod, LAI=data[\"LAI\"].to_numpy()\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "id": "424196e7-f7bc-4d23-bcfb-8a245fa22745", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Vmax25_canopy = [220.32138308]\n", + "Vmax25_sun = [112.19495779]\n", + "Vmax25_shade = [108.12642529]\n", + "Vmax_sun = [88.33176716]\n", + "Vmax_shade = [85.12858699]\n", + "Av_sun = [23.7468972]\n", + "Av_shade = [22.88576204]\n", + "Jmax25_sun = [213.09973077]\n", + "Jmax25_shade = [206.42733748]\n", + "J_sun = [91.28498756]\n", + "J_shade = [18.85937627]\n", + "Jmax_sun = [181.16701827]\n", + "Jmax_shade = [175.49447428]\n", + "Aj_sun = [15.25224592]\n", + "Aj_shade = [3.15109694]\n", + "Acanopy_sun = [15.25224592]\n", + "Acanopy_shade = [3.15109694]\n", + "gpp = [186.34124706]\n" + ] + } + ], + "source": [ + "print(f\"Vmax25_canopy = {two_leaf_assim_standard.Vmax25_canopy}\")\n", + "print(f\"Vmax25_sun = {two_leaf_assim_standard.Vmax25_sun}\")\n", + "print(f\"Vmax25_shade = {two_leaf_assim_standard.Vmax25_shade}\")\n", + "print(f\"Vmax_sun = {two_leaf_assim_standard.Vmax_sun}\")\n", + "print(f\"Vmax_shade = {two_leaf_assim_standard.Vmax_shade}\")\n", + "print(f\"Av_sun = {two_leaf_assim_standard.Av_sun}\")\n", + "print(f\"Av_shade = {two_leaf_assim_standard.Av_shade}\")\n", + "print(f\"Jmax25_sun = {two_leaf_assim_standard.Jmax25_sun}\")\n", + "print(f\"Jmax25_shade = {two_leaf_assim_standard.Jmax25_shade}\")\n", + "print(f\"J_sun = {two_leaf_assim_standard.J_sun}\")\n", + "print(f\"J_shade = {two_leaf_assim_standard.J_shade}\")\n", + "print(f\"Jmax_sun = {two_leaf_assim_standard.Jmax_sun}\")\n", + "print(f\"Jmax_shade = {two_leaf_assim_standard.Jmax_shade}\")\n", + "print(f\"Aj_sun = {two_leaf_assim_standard.Aj_sun}\")\n", + "print(f\"Aj_shade = {two_leaf_assim_standard.Aj_shade}\")\n", + "print(f\"Acanopy_sun = {two_leaf_assim_standard.Acanopy_sun}\")\n", + "print(f\"Acanopy_shade = {two_leaf_assim_standard.Acanopy_shade}\")\n", + "print(f\"gpp = {two_leaf_assim_standard.gpp}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6bf70634-b1ac-4810-a069-dc9c2c7cabc9", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "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.10.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tests/regression/two_leaf_canopy/test_two_leaf_model.py b/tests/regression/two_leaf_canopy/test_two_leaf_model.py new file mode 100644 index 00000000..8329766b --- /dev/null +++ b/tests/regression/two_leaf_canopy/test_two_leaf_model.py @@ -0,0 +1,108 @@ +"""TBC.""" + +import numpy as np +import pandas as pd +import pytest + +from pyrealm.pmodel import PModel +from pyrealm.pmodel.pmodel_environment import PModelEnvironment +from pyrealm.pmodel.two_leaf_irradience import TwoLeafAssimilation, TwoLeafIrradience + + +@pytest.fixture +def get_data(): + """Docstring.""" + forcing_data = pd.read_csv("pyrealm_build_data/two_leaf/merged_BE-Vie_data.csv") + forcing_data["time"] = pd.to_datetime(forcing_data["time"]) + + return forcing_data + + +@pytest.fixture +def solar_elevation(): + """Docstring.""" + # from Keiths code 0.982121, 0.9820922045388 + return np.array([0.98212117]) + + +@pytest.fixture +def solar_irradience(solar_elevation, get_data): + """Test.""" + + data = get_data.loc[(get_data["time"] == "2014-08-01 12:30:00")] + + PPFD = (data["ppfd"].to_numpy(),) + LAI = (data["LAI"].to_numpy(),) + PATM = data["patm"].to_numpy() + result = TwoLeafIrradience(solar_elevation, PPFD, LAI, PATM) + + return result + + +@pytest.fixture +def assimilation_single_day(solar_irradience, get_data): + """Test.""" + data = get_data.loc[(get_data["time"] == "2014-08-01 12:30:00")] + pmod_env = PModelEnvironment( + tc=data["tc"].to_numpy(), + patm=data["patm"].to_numpy(), + co2=data["co2"].to_numpy(), + vpd=data["vpd"].to_numpy(), + ) + + # Standard P Model + standard_pmod = PModel(pmod_env, kphio=1 / 8) + standard_pmod.estimate_productivity( + fapar=data["fapar"].to_numpy(), ppfd=data["ppfd"].to_numpy() + ) + + solar_irradience.calc_absorbed_irradience() + + assim = TwoLeafAssimilation( + standard_pmod, solar_irradience, leaf_area_index=data["LAI"].to_numpy() + ) + + assim.gpp_estimator() + + return assim + + +def test_two_leaf_irradience(solar_irradience, get_data): + """Docstring.""" + solar_irradience.calc_absorbed_irradience() + + assert np.allclose(solar_irradience.kb, np.array([0.6011949063494533])) + assert np.allclose(solar_irradience.kb_prime, np.array([0.55309931])) + assert np.allclose(solar_irradience.rho_h, np.array([0.04060739027615024])) + assert np.allclose(solar_irradience.rho_cb, np.array([0.03003319])) + assert np.allclose(solar_irradience.I_d, np.array([136.29450038])) + assert np.allclose(solar_irradience.I_b, np.array([714.70549962])) + assert np.allclose(solar_irradience.I_c, np.array([636.08707023])) + assert np.allclose(solar_irradience.Isun_beam, np.array([485.32862357])) + assert np.allclose(solar_irradience.Isun_diffuse, np.array([69.44255865])) + assert np.allclose(solar_irradience.Isun_scattered, np.array([25.43894747])) + assert np.allclose(solar_irradience.I_cshade, np.array([55.876940533221614])) + assert np.allclose(solar_irradience.I_csun, np.array([580.2101296936461])) + + +def test_two_leaf_assimilation(assimilation_single_day): + """Test.""" + assimilation = assimilation_single_day + assert np.allclose(assimilation.Vmax25_canopy, np.array([220.32138308])) + assert np.allclose(assimilation.Vmax25_sun, np.array([112.19495779])) + assert np.allclose(assimilation.Vmax25_shade, np.array([108.12642529])) + assert np.allclose(assimilation.Vmax_sun, np.array([88.33176716])) + assert np.allclose(assimilation.Vmax_shade, np.array([85.12858699])) + assert np.allclose(assimilation.Av_sun, np.array([23.7468972])) + assert np.allclose(assimilation.Av_shade, np.array([22.88576204])) + assert np.allclose(assimilation.Jmax25_sun, np.array([213.09973077])) + assert np.allclose(assimilation.Jmax25_shade, np.array([206.42733748])) + assert np.allclose(assimilation.Jmax_sun, np.array([181.16701827])) + assert np.allclose(assimilation.Jmax_shade, np.array([175.49447428])) + assert np.allclose(assimilation.J_sun, np.array([91.28498756])) + assert np.allclose(assimilation.J_shade, np.array([18.85937627])) + assert np.allclose(assimilation.Aj_sun, np.array([15.25224592])) + assert np.allclose(assimilation.Aj_shade, np.array([3.15109694])) + assert np.allclose(assimilation.Acanopy_sun, np.array([15.25224592])) + assert np.allclose(assimilation.Acanopy_shade, np.array([3.15109694])) + assert np.allclose(assimilation.gpp_estimate, np.array([186.34124706])) diff --git a/tests/unit/two_leaf_canopy/test_two_leaf_classes.py b/tests/unit/two_leaf_canopy/test_two_leaf_classes.py new file mode 100644 index 00000000..0c5e7d9f --- /dev/null +++ b/tests/unit/two_leaf_canopy/test_two_leaf_classes.py @@ -0,0 +1,181 @@ +"""TBC. + +TBC. +""" + +import numpy as np +import pytest + +from pyrealm.constants.core_const import CoreConst +from pyrealm.constants.two_leaf_canopy import TwoLeafConst +from pyrealm.pmodel.two_leaf_irradience import ( + TwoLeafAssimilation, + TwoLeafIrradience, +) + + +@pytest.fixture +def two_leaf_constants(): + """Fixture to provide the default constants.""" + return TwoLeafConst() + + +@pytest.fixture +def two_leaf_irradience(two_leaf_constants): + """Fixture to create a TwoLeafIrradience instance.""" + return TwoLeafIrradience( + beta_angle=np.array([0.6]), + ppfd=np.array([1000]), + leaf_area_index=np.array([2.0]), + patm=np.array([101325]), + constants=two_leaf_constants, + ) + + +def test_check_input_consistency(two_leaf_irradience): + """Test the _check_input_consistency method.""" + # Consistent shapes + assert two_leaf_irradience._check_input_consistency() is True + + # Inconsistent shapes + two_leaf_irradience.patm = np.array([101325, 101300]) + assert two_leaf_irradience._check_input_consistency() is False + + +def test_check_for_NaN(two_leaf_irradience): + """Test the _check_for_NaN method to identify NaN values.""" + # No NaNs + assert two_leaf_irradience._check_for_NaN() is True + + # NaN in beta_angle + two_leaf_irradience.beta_angle = np.array([np.nan]) + assert two_leaf_irradience._check_for_NaN() is False + + +def test_check_for_negative_values(two_leaf_irradience): + """Test the _check_for_negative_values method to identify negative values.""" + # No negative values + assert two_leaf_irradience._check_for_negative_values() is True + + # Negative value in leaf_area_index + two_leaf_irradience.leaf_area_index = np.array([-2.0]) + assert two_leaf_irradience._check_for_negative_values() is False + + +def test_initialization(two_leaf_irradience, two_leaf_constants): + """Test initialization of the TwoLeafIrradience class.""" + assert two_leaf_irradience.beta_angle.shape == (1,) + assert two_leaf_irradience.ppfd.shape == (1,) + assert two_leaf_irradience.leaf_area_index.shape == (1,) + assert two_leaf_irradience.patm.shape == (1,) + assert two_leaf_irradience.pass_checks is True + + +def test_calc_absorbed_irradience(two_leaf_irradience): + """Test the calc_absorbed_irradience method.""" + two_leaf_irradience.calc_absorbed_irradience() + + # Check if all attributes are calculated + attributes = [ + "kb", + "kb_prime", + "fd", + "rho_h", + "rho_cb", + "I_d", + "I_b", + "I_c", + "Isun_beam", + "Isun_diffuse", + "Isun_scattered", + "I_csun", + "I_cshade", + ] + for attr in attributes: + if attr == "rho_h": + pass + else: + assert hasattr(two_leaf_irradience, attr) + assert ( + getattr(two_leaf_irradience, attr).shape + == two_leaf_irradience.beta_angle.shape + ) + + +@pytest.fixture(scope="session") +def mock_pmodel(): + """Fixture to mock a PModel instance.""" + + class MockPModel: + vcmax = np.array([50.0]) + vcmax25 = np.array([45.0]) + optchi = type( + "MockOptimalChiABC", (), {"mc": np.array([0.9]), "mj": np.array([0.8])} + )() + env = type("MockEnv", (), {"tc": np.array([25.0]), "core_const": CoreConst()})() + core_const = type( + "MockCoreConst", (), {"k_c_molmass": 12.0, "core_const": CoreConst()} + )() + + # only required to allow test to work + subdaily_vcmax = np.array([50.0]) + subdaily_vcmax25 = np.array([45.0]) + optimal_chi = type( + "MockOptimalChiABC", (), {"mc": np.array([0.9]), "mj": np.array([0.8])} + )() + + return MockPModel() + + +@pytest.fixture +def two_leaf_assimilation(mock_pmodel, two_leaf_irradience): + """Fixture to create a TwoLeafAssimilation instance.""" + TLA = TwoLeafAssimilation( + pmodel=mock_pmodel, irrad=two_leaf_irradience, leaf_area_index=np.array([2.0]) + ) + return TLA + + +def test_initialization_assim(two_leaf_assimilation): + """Test initialization of the TwoLeafAssimilation class.""" + assert hasattr(two_leaf_assimilation, "vcmax_pmod") + assert hasattr(two_leaf_assimilation, "vcmax_pmod") + assert hasattr(two_leaf_assimilation, "vcmax25_pmod") + assert hasattr(two_leaf_assimilation, "optchi_obj") + assert hasattr(two_leaf_assimilation, "core_const") + assert isinstance(two_leaf_assimilation.irrad, TwoLeafIrradience) + + +def test_gpp_estimator(two_leaf_assimilation, two_leaf_irradience): + """Test the gpp_estimator method.""" + two_leaf_irradience.calc_absorbed_irradience() + two_leaf_assimilation.gpp_estimator() + + # Check if all GPP-related attributes are calculated + attributes = [ + "kv_Lloyd", + "Vmax25_canopy", + "Vmax25_sun", + "Vmax25_shade", + "Vmax_sun", + "Vmax_shade", + "Av_sun", + "Av_shade", + "Jmax25_sun", + "Jmax25_shade", + "Jmax_sun", + "Jmax_shade", + "J_sun", + "J_shade", + "Aj_sun", + "Aj_shade", + "Acanopy_sun", + "Acanopy_shade", + "gpp_estimate", + ] + for attr in attributes: + assert hasattr(two_leaf_assimilation, attr) + assert ( + getattr(two_leaf_assimilation, attr).shape + == two_leaf_assimilation.leaf_area_index.shape + ) diff --git a/tests/unit/two_leaf_canopy/test_two_leaf_irradience_functions.py b/tests/unit/two_leaf_canopy/test_two_leaf_irradience_functions.py new file mode 100644 index 00000000..2164cc70 --- /dev/null +++ b/tests/unit/two_leaf_canopy/test_two_leaf_irradience_functions.py @@ -0,0 +1,586 @@ +"""test_two_leaf_irradience_functions. + +This module contains pytest unit tests for validating the correctness of various +functions involved in the calculation of photosynthetic parameters and canopy +irradiance models. The functions tested include those related to beam extinction +coefficients, irradiance calculations, carboxylation scaling, electron transport +rates, and gross primary productivity (GPP). + +Each test function is parameterized with a variety of input scenarios to ensure +robustness and accuracy of the implemented algorithms. + +Tests are designed to cover edge cases, typical use cases, and mixed scenarios, +providing comprehensive coverage of the functionality. +""" + +import numpy as np +import pytest +from numpy.testing import assert_array_almost_equal + +from pyrealm.pmodel.two_leaf_irradience import ( + Jmax25, + Jmax25_temp_correction, + Vmax25_canopy, + Vmax25_shade, + Vmax25_sun, + assimilation_canopy, + assimilation_rate, + beam_extinction_coeff, + beam_irrad_unif_leaf_angle_dist, + beam_irradience, + beam_irradience_h_leaves, + canopy_extinction_coefficient, + canopy_irradience, + carboxylation_scaling_to_T, + diffuse_radiation, + electron_transport_rate, + fraction_of_diffuse_rad, + gross_primary_product, + photosynthetic_estimate, + scattered_beam_extinction_coeff, + scattered_beam_irradience, + shaded_absorbed_irrad, + sunlit_absorbed_irrad, + sunlit_beam_irrad, + sunlit_diffuse_irrad, + sunlit_scattered_irrad, +) + + +@pytest.mark.parametrize( + "beta_angle, k_sol_obs_angle, clip_angle, kb_numerator, expected_kb", + [ + # Test case 1: beta_angle > k_sol_obs_angle + ( + np.array([0.6, 0.7]), + 0.5, + 30, + 0.5, + np.array([0.5 / np.sin(0.6), 0.5 / np.sin(0.7)]), + ), + # Test case 2: beta_angle <= k_sol_obs_angle + (np.array([0.4, 0.3]), 0.5, 30, 0.5, np.array([30, 30])), + # Test case 3: Mixed case with some angles above and some below the threshold + ( + np.array([0.4, 0.6, 0.3, 0.7]), + 0.5, + 30, + 0.5, + np.array([30, 0.5 / np.sin(0.6), 30, 0.5 / np.sin(0.7)]), + ), + # Test case 4: Custom clip_angle and kb_numerator + (np.array([0.6, 0.4]), 0.5, 25, 0.6, np.array([0.6 / np.sin(0.6), 25])), + # Test case 5: All angles equal to the threshold + (np.array([0.5, 0.5]), 0.5, 30, 0.5, np.array([30, 30])), + ], +) +def test_beam_extinction_coeff( + beta_angle, k_sol_obs_angle, clip_angle, kb_numerator, expected_kb +): + """Test the beam_extinction_coeff function with various input scenarios.""" + result = beam_extinction_coeff( + beta_angle, k_sol_obs_angle, clip_angle, kb_numerator + ) + assert_array_almost_equal(result, expected_kb, decimal=5) + + +@pytest.mark.parametrize( + "beta_angle, k_sol_obs_angle, expected_kb_prime", + [ + # Test case 1: beta_angle > k_sol_obs_angle + (np.array([0.6, 0.7]), 0.5, np.array([0.46 / np.sin(0.6), 0.46 / np.sin(0.7)])), + ], +) +def test_scattered_beam_extinction_coeff( + beta_angle, k_sol_obs_angle, expected_kb_prime +): + """Test scattered_beam_extinction_coeff function with various input scenarios.""" + result = scattered_beam_extinction_coeff(beta_angle, k_sol_obs_angle) + assert_array_almost_equal(result, expected_kb_prime, decimal=5) + + +@pytest.mark.parametrize( + "patm, pa0, beta_angle, k_fa, expected_fd", + [ + # Test case 1: Basic case with standard pressure and moderate angles + ( + np.array([101325, 101300]), + 101325, + np.array([0.6, 0.7]), + 0.85, + np.array( + [ + (1 - 0.72 ** ((101325 / 101325) / np.sin(0.6))) + / ( + 1 + (0.72 ** ((101325 / 101325) / np.sin(0.6)) * (1 / 0.85 - 1)) + ), + (1 - 0.72 ** ((101300 / 101325) / np.sin(0.7))) + / ( + 1 + (0.72 ** ((101300 / 101325) / np.sin(0.7)) * (1 / 0.85 - 1)) + ), + ] + ), + ), + ], +) +def test_fraction_of_diffuse_rad(patm, pa0, beta_angle, k_fa, expected_fd): + """Test the fraction_of_diffuse_rad function with various input scenarios.""" + result = fraction_of_diffuse_rad(patm, pa0, beta_angle, k_fa) + assert_array_almost_equal(result, expected_fd, decimal=5) + + +@pytest.mark.parametrize( + "k_sigma, expected_rho_h", + [ + # Test case 1: Typical value for k_sigma + (0.5, 0.1715728752538099), + ], +) +def test_beam_irradience_h_leaves(k_sigma, expected_rho_h): + """Test the beam_irradience_h_leaves function with various input scenarios.""" + result = beam_irradience_h_leaves(k_sigma) + assert_array_almost_equal(result, expected_rho_h, decimal=5) + + +@pytest.mark.parametrize( + "rho_h, kb, expected_rho_cb", + [ + # Test case 1: Typical values for rho_h and kb + ( + np.array([0.1, 0.2]), + np.array([0.5, 1.0]), + np.array([0.064493015, 0.181269247]), + ), + ], +) +def test_beam_irrad_unif_leaf_angle_dist(rho_h, kb, expected_rho_cb): + """Test beam_irrad_unif_leaf_angle_dist function with various input scenarios.""" + result = beam_irrad_unif_leaf_angle_dist(rho_h, kb) + assert_array_almost_equal(result, expected_rho_cb, decimal=5) + + +@pytest.mark.parametrize( + "fd, ppfd, expected_I_d", + [ + # Test case 1: Typical values for fd and ppfd + (np.array([0.3, 0.5]), np.array([1000, 800]), np.array([300.0, 400.0])), + # Test case 2: mixed cases which should result in 0 diffuse radiation + ( + np.array([0.0, 0.5, -0.5, 0.5]), + np.array([1000, 0, 1000, -1000]), + np.array([0.0, 0.0, 0.0, 0.0]), + ), + ], +) +def test_diffuse_radiation(fd, ppfd, expected_I_d): + """Test the diffuse_radiation function with various input scenarios.""" + result = diffuse_radiation(fd, ppfd) + assert_array_almost_equal(result, expected_I_d, decimal=5) + + +@pytest.mark.parametrize( + "ppfd, fd, expected_I_b", + [ + # Test case 1: Typical values for ppfd and fd + (np.array([1000, 800]), np.array([0.3, 0.5]), np.array([700.0, 400.0])), + ], +) +def test_beam_irradience(ppfd, fd, expected_I_b): + """Test the beam_irradience function with various input scenarios.""" + result = beam_irradience(ppfd, fd) + assert_array_almost_equal(result, expected_I_b, decimal=5) + + +@pytest.mark.parametrize( + "I_b, kb, kb_prime, rho_cb, leaf_area_index, k_sigma, expected_I_bs", + [ + # Test case 1: Typical values for all parameters + ( + np.array([500.0, 300.0]), + np.array([0.7, 0.5]), + np.array([0.6, 0.4]), + np.array([0.2, 0.3]), + np.array([2.0, 1.5]), + np.array([0.1, 0.2]), + np.array([72.13125477, 45.91123081]), + ), + ], +) +def test_scattered_beam_irradience( + I_b, kb, kb_prime, rho_cb, leaf_area_index, k_sigma, expected_I_bs +): + """Test the scattered_beam_irradience function with various input scenarios.""" + result = scattered_beam_irradience( + I_b, kb, kb_prime, rho_cb, leaf_area_index, k_sigma + ) + assert_array_almost_equal(result, expected_I_bs, decimal=5) + + +@pytest.mark.parametrize( + "rho_cb, I_b, kb_prime, I_d, leaf_area_index, k_rho_cd, expected_I_c", + [ + # Test case 1: Typical values for all parameters + ( + np.array([0.2, 0.3]), + np.array([500.0, 300.0]), + np.array([0.6, 0.4]), + np.array([200.0, 150.0]), + np.array([2.0, 1.5]), + 0.5, + np.array([349.402894, 128.5886837]), + ), + ], +) +def test_canopy_irradience( + rho_cb, I_b, kb_prime, I_d, leaf_area_index, k_rho_cd, expected_I_c +): + """Test the canopy_irradience function with various input scenarios.""" + result = canopy_irradience(rho_cb, I_b, kb_prime, I_d, leaf_area_index, k_rho_cd) + assert_array_almost_equal(result, expected_I_c, decimal=5) + + +@pytest.mark.parametrize( + "I_b, k_sigma, kb, leaf_area_index, expected_Isun_beam", + [ + # Test case 1: Typical values for all parameters + ( + np.array([500.0, 300.0]), + 0.2, + np.array([0.7, 0.5]), + np.array([2.0, 1.5]), + np.array([301.3612144, 126.6320273]), + ), + ], +) +def test_sunlit_beam_irrad(I_b, k_sigma, kb, leaf_area_index, expected_Isun_beam): + """Test the sunlit_beam_irrad function with various input scenarios.""" + result = sunlit_beam_irrad(I_b, k_sigma, kb, leaf_area_index) + assert_array_almost_equal(result, expected_Isun_beam, decimal=5) + + +@pytest.mark.parametrize( + "I_d, k_rho_cd, k_kd_prime, kb, leaf_area_index, expected_Isun_diffuse", + [ + # Test case 1: Typical values for all parameters + ( + np.array([200.0, 150.0]), + np.array([0.3, 0.4]), + 0.5, + np.array([0.7, 0.5]), + np.array([2.0, 1.5]), + np.array([53.04145272, 34.95914279]), + ), + ], +) +def test_sunlit_diffuse_irrad( + I_d, k_rho_cd, k_kd_prime, kb, leaf_area_index, expected_Isun_diffuse +): + """Test the sunlit_diffuse_irrad function with various input scenarios.""" + result = sunlit_diffuse_irrad(I_d, k_rho_cd, k_kd_prime, kb, leaf_area_index) + assert_array_almost_equal(result, expected_Isun_diffuse, decimal=5) + + +@pytest.mark.parametrize( + "I_b, rho_cb, kb_prime, kb, leaf_area_index, k_sigma, expected_Isun_scattered", + [ + # Test case 1: Typical values for all parameters + ( + np.array([500.0, 300.0]), + np.array([0.2, 0.3]), + np.array([0.6, 0.5]), + np.array([0.7, 0.6]), + np.array([2.0, 1.5]), + 0.3, + np.array([6.545100366, -10.52110801]), + ), + ], +) +def test_sunlit_scattered_irrad( + I_b, rho_cb, kb_prime, kb, leaf_area_index, k_sigma, expected_Isun_scattered +): + """Test the sunlit_scattered_irrad function with various input scenarios.""" + result = sunlit_scattered_irrad(I_b, rho_cb, kb_prime, kb, leaf_area_index, k_sigma) + assert_array_almost_equal(result, expected_Isun_scattered, decimal=6) + + +@pytest.mark.parametrize( + "Isun_beam, Isun_diffuse, Isun_scattered, expected_I_csun", + [ + # Test case 1: Typical values for all components + ( + np.array([200.0, 150.0]), + np.array([50.0, 30.0]), + np.array([100.0, 80.0]), + np.array([350.0, 260.0]), + ), + ], +) +def test_sunlit_absorbed_irrad( + Isun_beam, Isun_diffuse, Isun_scattered, expected_I_csun +): + """Test the sunlit_absorbed_irrad function with various input scenarios.""" + result = sunlit_absorbed_irrad(Isun_beam, Isun_diffuse, Isun_scattered) + assert_array_almost_equal(result, expected_I_csun, decimal=5) + + +@pytest.mark.parametrize( + "beta_angle, k_sol_obs_angle, I_c, I_csun, expected_I_cshade", + [ + # Test case 1: Typical values where beta_angle > k_sol_obs_angle + ( + np.array([0.7, 0.8]), + 0.6, + np.array([400.0, 350.0]), + np.array([200.0, 150.0]), + np.array([200.0, 200.0]), + ), + # Test case 2: beta_angle <= k_sol_obs_angle, resulting in zero shaded + # irradiance + ( + np.array([0.5, 0.4]), + 0.6, + np.array([400.0, 350.0]), + np.array([200.0, 150.0]), + np.array([0.0, 0.0]), + ), + ], +) +def test_shaded_absorbed_irrad( + beta_angle, k_sol_obs_angle, I_c, I_csun, expected_I_cshade +): + """Test the shaded_absorbed_irrad function with various input scenarios.""" + result = shaded_absorbed_irrad(beta_angle, k_sol_obs_angle, I_c, I_csun) + assert_array_almost_equal(result, expected_I_cshade, decimal=5) + + +@pytest.mark.parametrize( + "vcmax_pmod, expected_kv_Lloyd", + [ + # Test case 1: Typical value for vcmax_pmod + (np.array([50.0, 60.0]), np.array([0.142487643, 0.156891625])), + ], +) +def test_canopy_extinction_coefficient(vcmax_pmod, expected_kv_Lloyd): + """Test the canopy_extinction_coefficient function with various input scenarios.""" + result = canopy_extinction_coefficient(vcmax_pmod) + assert_array_almost_equal(result, expected_kv_Lloyd, decimal=8) + + +@pytest.mark.parametrize( + "leaf_area_index, vcmax25_pmod, kv, expected_Vmax25_canopy", + [ + # Test case 1: Typical values for all parameters + ( + np.array([2.0, 3.0]), + np.array([50.0, 60.0]), + np.array([0.5, 0.6]), + np.array([78.69386806, 135.3565092]), + ), + ], +) +def test_Vmax25_canopy(leaf_area_index, vcmax25_pmod, kv, expected_Vmax25_canopy): + """Test the Vmax25_canopy function with various input scenarios.""" + result = Vmax25_canopy(leaf_area_index, vcmax25_pmod, kv) + assert_array_almost_equal(result, expected_Vmax25_canopy, decimal=6) + + +@pytest.mark.parametrize( + "leaf_area_index, vcmax25_pmod, kv, kb, expected_Vmax25_sun", + [ + # Test case 1: Typical values for all parameters + ( + np.array([2.0, 3.0]), + np.array([50.0, 60.0]), + np.array([0.5, 0.6]), + np.array([0.7, 0.8]), + np.array([44.75954636, 57.0127759]), + ), + ], +) +def test_Vmax25_sun(leaf_area_index, vcmax25_pmod, kv, kb, expected_Vmax25_sun): + """Test the Vmax25_sun function with various input scenarios.""" + result = Vmax25_sun(leaf_area_index, vcmax25_pmod, kv, kb) + assert_array_almost_equal(result, expected_Vmax25_sun, decimal=8) + + +@pytest.mark.parametrize( + "Vmax25_canopy, Vmax25_sun, expected_Vmax25_shade", + [ + # Test case 1: Typical values for Vmax25_canopy and Vmax25_sun + (np.array([100.0, 150.0]), np.array([60.0, 90.0]), np.array([40.0, 60.0])), + ], +) +def test_Vmax25_shade(Vmax25_canopy, Vmax25_sun, expected_Vmax25_shade): + """Test the Vmax25_shade function with various input scenarios.""" + result = Vmax25_shade(Vmax25_canopy, Vmax25_sun) + assert_array_almost_equal(result, expected_Vmax25_shade, decimal=8) + + +@pytest.mark.parametrize( + "Vmax25, tc, expected_Vmax", + [ + # Test case 1: Typical values for Vmax25 and temperature + ( + np.array([100.0, 150.0]), + np.array([25.0, 30.0]), + np.array([100.0, 230.9566403]), + ), + # Test case 2: Temperature below 25C + ( + np.array([100.0, 150.0]), + np.array([20.0, 15.0]), + np.array([63.99758172, 60.49060866]), + ), + ], +) +def test_carboxylation_scaling_to_T(Vmax25, tc, expected_Vmax): + """Test the carboxylation_scaling_to_T function with various input scenarios.""" + result = carboxylation_scaling_to_T(Vmax25, tc) + assert_array_almost_equal(result, expected_Vmax, decimal=6) + + +@pytest.mark.parametrize( + "Vmax, mc, expected_Av", + [ + # Test case 1: Typical values for Vmax and mc + (np.array([100.0, 150.0]), np.array([0.8, 0.9]), np.array([80.0, 135.0])), + ], +) +def test_photosynthetic_estimate(Vmax, mc, expected_Av): + """Test the photosynthetic_estimate function with various input scenarios.""" + result = photosynthetic_estimate(Vmax, mc) + assert_array_almost_equal(result, expected_Av, decimal=8) + + +@pytest.mark.parametrize( + "Vmax25, expected_Jmax25", + [ + # Test case 1: Typical values for Vmax25 + (np.array([50.0, 100.0]), np.array([111.1, 193.1])), + ], +) +def test_Jmax25(Vmax25, expected_Jmax25): + """Test the Jmax25 function with various input scenarios.""" + result = Jmax25(Vmax25) + assert_array_almost_equal(result, expected_Jmax25, decimal=8) + + +@pytest.mark.parametrize( + "Jmax25, tc, expected_Jmax", + [ + # Test case 1: Typical values for Jmax25 and temperature + ( + np.array([100.0, 150.0]), + np.array([25.0, 30.0]), + np.array([100.0, 201.0647139]), + ), + # Test case 2: Temperature below 25C + ( + np.array([100.0, 150.0]), + np.array([20.0, 15.0]), + np.array([73.86055721, 80.97433885]), + ), + ], +) +def test_Jmax25_temp_correction(Jmax25, tc, expected_Jmax): + """Test the Jmax25_temp_correction function with various input scenarios.""" + result = Jmax25_temp_correction(Jmax25, tc) + assert_array_almost_equal(result, expected_Jmax, decimal=8) + + +@pytest.mark.parametrize( + "Jmax, I_c, expected_J", + [ + # Test case 1: Typical values for Jmax and I_c + ( + np.array([100.0, 150.0]), + np.array([500.0, 400.0]), + np.array([59.02777778, 69.8630137]), + ), + ], +) +def test_electron_transport_rate(Jmax, I_c, expected_J): + """Test the electron_transport_rate function with various input scenarios.""" + result = electron_transport_rate(Jmax, I_c) + assert_array_almost_equal(result, expected_J, decimal=6) + + +@pytest.mark.parametrize( + "mj, J, expected_A", + [ + # Test case 1: Typical values for mj and J + ( + np.array([0.5, 0.7]), + np.array([100.0, 200.0]), + np.array( + [ + 12.5, # Calculated value: 0.5 * 100 / 4 + 35.0, # Calculated value: 0.7 * 200 / 4 + ] + ), + ), + ], +) +def test_assimilation_rate(mj, J, expected_A): + """Test the assimilation_rate function with various input scenarios.""" + result = assimilation_rate(mj, J) + assert_array_almost_equal(result, expected_A, decimal=8) + + +@pytest.mark.parametrize( + "Aj, Av, beta_angle, solar_obscurity_angle, expected_Acanopy", + [ + # Test case 1: Typical values with beta_angle > solar_obscurity_angle + ( + np.array([50.0, 100.0]), + np.array([60.0, 90.0]), + np.array([0.7, 0.8]), + 0.6, + np.array([50.0, 90.0]), + ), + # Test case 2: beta_angle < solar_obscurity_angle, resulting in zero + # assimilation + ( + np.array([50.0, 100.0]), + np.array([60.0, 90.0]), + np.array([0.5, 0.4]), + 0.6, + np.array([0.0, 0.0]), + ), + # Test case 3: Aj equal to Av + ( + np.array([70.0, 120.0]), + np.array([70.0, 120.0]), + np.array([0.7, 0.8]), + 0.6, + np.array([70.0, 120.0]), + ), + # Test case 4: Aj > Av + ( + np.array([80.0, 130.0]), + np.array([60.0, 110.0]), + np.array([0.7, 0.8]), + 0.6, + np.array([60.0, 110.0]), + ), + ], +) +def test_assimilation_canopy( + Aj, Av, beta_angle, solar_obscurity_angle, expected_Acanopy +): + """Test the assimilation_canopy function with various input scenarios.""" + result = assimilation_canopy(Aj, Av, beta_angle, solar_obscurity_angle) + assert_array_almost_equal(result, expected_Acanopy, decimal=8) + + +@pytest.mark.parametrize( + "k_c_molmass, Acanopy_sun, Acanopy_shade, expected_gpp", + [ + # Test case 1: Typical values for all inputs + (12.0, np.array([10.0, 20.0]), np.array([5.0, 10.0]), np.array([125.0, 250.0])) + ], +) +def test_gross_primary_product(k_c_molmass, Acanopy_sun, Acanopy_shade, expected_gpp): + """Test the gross_primary_product function with various input scenarios.""" + result = gross_primary_product(k_c_molmass, Acanopy_sun, Acanopy_shade) + assert_array_almost_equal(result, expected_gpp, decimal=8)