Annual fAPAR_max and LAI_max #403
Conversation
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## develop #403 +/- ##
===========================================
- Coverage 96.26% 94.71% -1.55%
===========================================
Files 34 35 +1
Lines 2652 2708 +56
===========================================
+ Hits 2553 2565 +12
- Misses 99 143 +44 ☔ View full report in Codecov by Sentry. |
|
@davidorme , could you have a look and sense-check what's there? As you suggested in #407 (I hope I understood that right) I used |
There was a problem hiding this comment.
Some optimisations and a couple of issues. I think we're on the right track here. I did wonder about making the class __init__ accept daily/subdaily observations and moving the resampling in the from_Pmodel into there, but on balance I think it is right to have the straightforward usage as __init__ and then use factory methods to handle other core inputs.
| """ | ||
|
|
||
| # Extract years from datetimes | ||
| all_years = [np.datetime64(i, "Y") for i in datetimes] |
There was a problem hiding this comment.
| all_years = [np.datetime64(i, "Y") for i in datetimes] | |
| all_years = datetimes.astype("datetime64[Y]") |
| """Class to compute the fAPAR_max and annual peak Leaf Area Index (LAI).""" | ||
|
|
||
| import numpy as np | ||
| from numpy.ma.core import zeros |
There was a problem hiding this comment.
Why not just use np.zeros? I think ma is the masked array module and that has significant performance issues.
| z: float = 12.227, | ||
| k: float = 0.5, | ||
| ) -> None: |
There was a problem hiding this comment.
We can bundle these constants (and the hardcoded ones below) into a pyrealm.constants class - probably a new one for PhenologyConsts?
| self.faparmax = -9999 * np.ones(np.shape(fapar_waterlim)) | ||
| self.energylim = -9999 * np.ones(np.shape(fapar_waterlim)) | ||
| self.annual_precip_molar = annual_total_precip | ||
|
|
||
| for i in range(len(fapar_waterlim)): | ||
| if fapar_waterlim[i] < fapar_energylim[i]: | ||
| self.faparmax[i] = fapar_waterlim[i] | ||
| self.energylim[i] = False | ||
| else: | ||
| self.faparmax[i] = fapar_energylim[i] | ||
| self.energylim[i] = True |
There was a problem hiding this comment.
| self.faparmax = -9999 * np.ones(np.shape(fapar_waterlim)) | |
| self.energylim = -9999 * np.ones(np.shape(fapar_waterlim)) | |
| self.annual_precip_molar = annual_total_precip | |
| for i in range(len(fapar_waterlim)): | |
| if fapar_waterlim[i] < fapar_energylim[i]: | |
| self.faparmax[i] = fapar_waterlim[i] | |
| self.energylim[i] = False | |
| else: | |
| self.faparmax[i] = fapar_energylim[i] | |
| self.energylim[i] = True | |
| self.fapar_max = np.minimum(fapar_waterlim, fapar_energylim) | |
| """Maximum fapar given water or energy limitation.""" | |
| self.energy_limited = fapar_energylim < fapar_waterlim | |
| """Is fapar_max limited by water or energy.""" | |
| self.annual_precip_molar = annual_total_precip | |
| """The annual precipitation in moles.""" |
E.g.
In [8]: fapar_waterlim = np.arange(10)
In [10]: fapar_waterlim
Out[10]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [12]: fapar_energylim = np.arange(9, -1, -1)
In [13]: fapar_energylim
Out[13]: array([9, 8, 7, 6, 5, 4, 3, 2, 1, 0])
In [15]: np.minimum(fapar_waterlim, fapar_energylim)
Out[15]: array([0, 1, 2, 3, 4, 4, 3, 2, 1, 0])
In [16]: fapar_waterlim < fapar_energylim
Out[16]: array([ True, True, True, True, True, False, False, False, False,
False])|
|
||
| annual_total_potential_gpp = get_annual( | ||
| pmodel.gpp, datetimes, growing_season, "total" |
There was a problem hiding this comment.
I think this need to check we are dealing with whole years in the input?
| annual_x = zeros(len(years)) | ||
|
|
||
| if method == "total": | ||
| for i in range(len(years)): | ||
| annual_x[i] = sum( | ||
| daily_x[growing_season & (years_by_day == years[i])].astype(np.int64) | ||
| ) | ||
| elif method == "mean": | ||
| for i in range(len(years)): | ||
| annual_x[i] = np.mean( | ||
| daily_x[growing_season & (years_by_day == years[i])].astype(np.int64) | ||
| ) |
There was a problem hiding this comment.
I think we can avoid this loop by reshaping the inputs?
Oh no, we can't. Leap years. Hang on - does SubdailyScalar handle leap years? Must do? Need to check on this.
There was a problem hiding this comment.
I was mis-remembering this - SubdailyScalar reshapes evenly spaced (e.g.) hourly data into days by hours so is unaffected by leap years.
| / (1.6 * annual_mean_vpd * annual_total_potential_gpp) | ||
| ) | ||
|
|
||
| self.faparmax = -9999 * np.ones(np.shape(fapar_waterlim)) |
There was a problem hiding this comment.
Generally, I'd use np.nan in cases like this. Here though, I think we can make this more efficient - see suggestion below.
| leaves. [mol m^{-2} year^{-1}] | ||
| k: Light extinction coefficient. | ||
| """ | ||
|
|
There was a problem hiding this comment.
We should _check_shapes on the inputs. All should be equal (although aridity_index can be scalar, that use case is permitted by _check_shape).
| x: Array of values to be converted to annual values. Should be either daily ( | ||
| same datetimes as growing_season) or subdaily (same datetimes as datetimes | ||
| array) | ||
| datetimes: Datetimes of the measurements as np.datetime62 arrays. |
There was a problem hiding this comment.
| datetimes: Datetimes of the measurements as np.datetime62 arrays. | |
| datetimes: Datetimes of the measurements as np.datetime64 arrays. |
| scaler = SubdailyScaler(datetimes) | ||
| scaler.set_nearest(np.timedelta64(12, "h")) |
There was a problem hiding this comment.
We can't use the SubdailyScalar on data that isn't sampled at subdaily frequency. The code explicitly checks for the spacing of observations within a day and breaks with coarser data (daily, monthly). I think we can still use the same approach if we nest this within the if statements below and solve it for the daily+ mode in another way.
There was a problem hiding this comment.
We should add a failure mode to SubdailyScalar for that case!
Description
Implements the computation of fAPAR_max and LAI_max, and provides examples usages of the API (including test data).
Fixes #348 .
Type of change
Key checklist
pre-commitchecks:$ pre-commit run -a$ poetry run pytestFurther checks