1"""modvege.py
2
3"""
4
5import numpy as np
6
7import climag.modvege_consumption as cm
8import climag.modvege_lib as lm
9
10np.seterr("raise")
11
12
[docs]
13def sum_of_temperature_thresholds(timeseries, params) -> dict[str, float]:
14 """Calculate sum of temperatures
15
16 Parameters
17 ----------
18 timeseries : pandas.DataFrame
19 Input meteorological time series data
20 params : dict
21 A dictionary of input parameters
22
23 Returns
24 -------
25 dict[str, float]
26 Dictionary of the sum of temperature thresholds
27
28 Notes
29 -----
30 Calculate sum of temperatures at:
31 - the beginning of the reproductive period (ST₁) [°C d]
32 - the end of the reproductive period (ST₂) [°C d]
33 - the beginning of the grazing season (STg₁) [°C d]
34 - the end of the grazing season (STg₂) [°C d]
35 - the beginning of harvest (STh₁) [°C d]
36
37 Nolan and Flanagan (2020) define the thermal growing season length as the
38 number of days between the first occurrence of at least six consecutive
39 days with a daily mean temperature of > 5°C and the first occurrence of at
40 least six consecutive days with a daily mean temperature of < 5°C.
41
42 If the temperatures are too low (i.e. no six consecutive days > 5°C) to
43 calculate the start date of the growing season, assume it is the 15th
44 March, which is the median date for the midlands and part of northern
45 Ireland (Collins and Cummins, 1996; Connaughton, 1973). This is also
46 the date when cows are out full time according to Teagasc recommendations
47 (Kavanagh, 2016).
48
49 Calculating the end date of the growing season is not straightforward, as
50 the temperatues may be too high for there to be six consecutive days
51 < 5°C, or these six consecutive days may happen very early in the year,
52 and may even be before the growing season start date.
53
54 Grazing season calculations based solely on temperature do not consider
55 the delay before sufficient plant cover is available to support grazing
56 animals or the ability of animals and machinery to pass over land (Nolan
57 and Flanagan, 2020; Collins and Cummins, 1996; based on Keane, 1982).
58
59 The beginning of the grazing season has a delay of 5-15 days after the
60 start of the growing season based on Broad and Hough (1993).
61 A delay of 10 days is used to allow sufficient reproductive growth.
62 ~The end date of the grazing season is determined using the Smith formula
63 for calculating the grazing season length (Collins and Cummins, 1996;
64 based on Smith, 1976).~
65
66 The end date of the grazing season cannot exceed 1st December. Livestock
67 are assumed to be fully housed by 22nd November based on Teagasc
68 recommendations (Kavanagh, 2016). The mean latest autumn closing date is
69 3rd December, with a two-week interval of 26th November to 9th December
70 (Looney et al., 2021).
71
72 Grazing will only take place if there is sufficient biomass available; if
73 the residual biomass has a height of less than 5 cm, no grazing or
74 harvesting will take place. Therefore, the amount of ingested and
75 harvested biomass are mainly influenced by the environmental factors that
76 affect growth, such as temperature, radiation, and precipitation.
77
78 The beginning of harvest is assumed to be one day before the grazing
79 season ends. Grazing costs less than indoor feeding, so starting the
80 harvest just a day before the end of the grazing season ensures grazing is
81 maximised.
82 The end of harvest is the same as the end of the grazing season.
83 """
84 st_thresholds = {}
85
86 timeseries.sort_values(by=["time"], inplace=True)
87 timeseries.set_index("time", inplace=True)
88
89 # generate index
90 try:
91 # for multi-year time series, this should correspond to the day
92 # of the year
93 timeseries["idx"] = timeseries["doy"] - 1
94 except KeyError:
95 # some datasets do not use standard calendars, so day of the year
96 # should not be used
97 timeseries["idx"] = range(len(timeseries))
98
99 # return only mean values above t_0, and subtract by t_0
100 timeseries.loc[(timeseries["T"] >= params["t_0"]), "Tg"] = (
101 timeseries["T"] - params["t_0"]
102 )
103 # fill NaN values
104 timeseries[["Tg"]] = timeseries[["Tg"]].fillna(value=0)
105
106 for year in timeseries.index.year.unique():
107 st_thresholds[year] = {}
108
109 # # grazing season length using the Smith formula
110 # grazing_season = round(
111 # 29.3 * np.mean(timeseries.loc[str(year)]["T"]) -
112 # 0.1 * np.sum(timeseries.loc[str(year)]["PP"]) +
113 # 19.5
114 # )
115 # # adjust the length if the dataset has 360 days/year
116 # if len(timeseries.loc[str(year)]) == 360:
117 # grazing_season -= 5
118
119 # sum of temperatures at the start
120 try:
121 start = list(
122 timeseries.loc[str(year)]["T"]
123 .rolling(6)
124 .apply(lambda x: all(x > 5.0))
125 ).index(1.0)
126 # the lowest possible start index would be 5
127 # so force the lowest value to zero
128 if start == 5:
129 start = 0
130 except ValueError:
131 # if the temperatures are too low for the start of the growing
132 # season to be calculated, assume it is on 15th March
133 start = int(timeseries.loc[f"{year}-03-15"]["idx"])
134
135 # beginning of the reproductive period
136 st_thresholds[year]["st_1"] = timeseries.loc[str(year)]["Tg"].cumsum()[
137 start
138 ]
139
140 # beginning of the grazing season
141 st_thresholds[year]["st_g1"] = timeseries.loc[str(year)][
142 "Tg"
143 ].cumsum()[start + 10]
144
145 # end of the reproductive period
146 st_thresholds[year]["st_2"] = timeseries.loc[str(year)]["Tg"].cumsum()[
147 -1
148 ]
149
150 # end of the grazing and harvesting season
151 # use the calculated grazing season length
152 # if growing season continues in December, end grazing on 1st December
153 # grazing_end = min(
154 # int(timeseries.loc[f"{year}-12-01"]["idx"]),
155 # start + 10 + grazing_season
156 # )
157 grazing_end = int(timeseries.loc[f"{year}-12-01"]["idx"])
158 st_thresholds[year]["st_g2"] = timeseries.loc[str(year)][
159 "Tg"
160 ].cumsum()[grazing_end]
161
162 # beginning of harvest
163 st_thresholds[year]["st_h1"] = timeseries.loc[str(year)][
164 "Tg"
165 ].cumsum()[grazing_end - 1]
166
167 timeseries.reset_index(inplace=True)
168 timeseries.drop(columns=["Tg", "idx"], inplace=True)
169
170 return st_thresholds
171
172
[docs]
173def modvege(params, tseries, endday=365, t_init=None) -> dict[str, float]:
174 """ModVege model as a function
175
176 Parameters
177 ----------
178 params : dict
179 Parameters (constants)
180 tseries : pandas.DataFrame
181 Time series meteorological data
182 endday : int
183 Number of days of the year (default is 365)
184
185 Returns
186 -------
187 dict[str, float]
188 Dictionary of results
189
190 Notes
191 -----
192 Results:
193 - Green vegetative biomass [kg DM ha⁻¹]
194 - Dead vegetative biomass [kg DM ha⁻¹]
195 - Green reproductive biomass [kg DM ha⁻¹]
196 - Dead reproductive biomass [kg DM ha⁻¹]
197 - Total standing biomass [kg DM ha⁻¹]
198 - Potential growth [kg DM ha⁻¹]
199 - Total growth [kg DM ha⁻¹]
200 - Ingested biomass [kg DM ha⁻¹]
201 - Harvested biomass [kg DM ha⁻¹]
202 - Leaf area index [dimensionless]
203 - Water reserves [mm]
204 - Actual evapotranspiration [mm]
205 - Environmental limitation of growth [dimensionless]
206 - Reproductive function [dimensionless]
207 """
208 st_thresholds = sum_of_temperature_thresholds(
209 timeseries=tseries, params=params
210 )
211
212 # dictionary of outputs
213 outputs_dict = {
214 "time": [],
215 "bm_gv": [],
216 "bm_gr": [],
217 "bm_dv": [],
218 "bm_dr": [],
219 "age_gv": [],
220 "age_gr": [],
221 "age_dv": [],
222 "age_dr": [],
223 "bm": [],
224 "pgro": [],
225 "gro": [],
226 "i_bm": [],
227 "h_bm": [],
228 "c_bm": [],
229 "env": [],
230 "lai": [],
231 "aet": [],
232 "wr": [],
233 "sen_gv": [],
234 "sen_gr": [],
235 "abs_dv": [],
236 "abs_dr": [],
237 "omd_gv": [],
238 "omd_gr": [],
239 }
240
241 # dictionary to store intermediate time series values
242 ts_vals = {}
243
244 # daily loop
245 for i in range(endday):
246 # initialise starting parameters
247 # standing biomass, biomass age, and water reserves
248 # assume that the initial value of WR = WHC
249 if i == 0:
250 (
251 ts_vals["bm_gv"],
252 ts_vals["bm_gr"],
253 ts_vals["bm_dv"],
254 ts_vals["bm_dr"],
255 ts_vals["age_gv"],
256 ts_vals["age_gr"],
257 ts_vals["age_dv"],
258 ts_vals["age_dr"],
259 ts_vals["wr"],
260 ) = (
261 params["bm_gv"],
262 params["bm_gr"],
263 params["bm_dv"],
264 params["bm_dr"],
265 params["age_gv"],
266 params["age_gr"],
267 params["age_dv"],
268 params["age_dr"],
269 params["wr"],
270 )
271
272 # initialise ingested/harvested biomass and temperature sum
273 if tseries["time"][i].dayofyear == 1:
274 # reset to zero on the first day of the year
275 ts_vals["i_bm"], ts_vals["h_bm"], ts_vals["st"] = 0.0, 0.0, 0.0
276 # sum of temperature thresholds
277 for key in st_thresholds[tseries["time"][i].year]:
278 params[key] = st_thresholds[tseries["time"][i].year][key]
279
280 # 10-d moving average temperature (T_m10)
281 ts_vals["t_m10"] = lm.ten_day_moving_avg_temperature(
282 day=(i + 1), t_ts=tseries["T"], t_init=t_init
283 )
284
285 # sum of temperatures (ST)
286 ts_vals["st"] = lm.sum_of_temperatures(
287 params=params, ts_vals=ts_vals, t_ts=tseries["T"], day=(i + 1)
288 )
289
290 # temperature function (f(T))
291 ts_vals["f_t"] = lm.temperature_function(
292 ts_vals=ts_vals, params=params
293 )
294
295 # seasonal effect (SEA)
296 ts_vals["sea"] = lm.seasonal_effect(params=params)
297
298 # leaf area index (LAI)
299 ts_vals["lai"] = lm.leaf_area_index(ts_vals=ts_vals, params=params)
300
301 # actual evapotranspiration (AET)
302 ts_vals["aet"] = lm.actual_evapotranspiration(
303 pet=tseries["PET"][i], ts_vals=ts_vals
304 )
305
306 # water reserves (WR)
307 ts_vals["wr"] = lm.water_reserves(
308 ts_vals=ts_vals, params=params, precipitation=tseries["PP"][i]
309 )
310
311 # water stress (W)
312 ts_vals["w"] = lm.water_stress(ts_vals=ts_vals, params=params)
313
314 # water stress function (f(W))
315 ts_vals["f_w"] = lm.water_stress_function(
316 ts_vals=ts_vals, pet=tseries["PET"][i]
317 )
318
319 # environmental limitation of growth (ENV)
320 ts_vals["env"] = lm.environmental_limitation(
321 ts_vals=ts_vals, params=params, par_i=tseries["PAR"][i]
322 )
323
324 # potential growth (PGRO)
325 ts_vals["pgro"] = lm.potential_growth(
326 par_i=tseries["PAR"][i], ts_vals=ts_vals, params=params
327 )
328
329 # total growth (GRO)
330 ts_vals["gro"] = lm.total_growth(ts_vals=ts_vals)
331
332 # reproductive function (REP)
333 # grazing always takes place during the grazing season if the
334 # stocking rate is > 0
335 ts_vals["rep"] = lm.reproductive_function(
336 params=params, ts_vals=ts_vals
337 )
338
339 # senescence (SEN)
340 lm.senescence(
341 ts_vals=ts_vals, params=params, temperature=tseries["T"][i]
342 )
343
344 # abscission (ABS)
345 lm.abscission(
346 ts_vals=ts_vals, params=params, temperature=tseries["T"][i]
347 )
348
349 # standing biomass (BM) and biomass age (AGE)
350 lm.standing_biomass(ts_vals=ts_vals, params=params)
351 lm.biomass_age(
352 temperature=tseries["T"][i], params=params, ts_vals=ts_vals
353 )
354
355 # organic matter digestibility (OMD)
356 cm.organic_matter_digestibility(ts_vals=ts_vals, params=params)
357
358 # ingested biomass
359 cm.biomass_ingestion(ts_vals=ts_vals, params=params)
360
361 # harvested biomass
362 cm.biomass_harvest(ts_vals=ts_vals, params=params)
363
364 # recover output streams
365 outputs_dict["time"].append(tseries["time"][i])
366
367 # net standing biomass
368 outputs_dict["bm"].append(
369 ts_vals["bm_gv"]
370 + ts_vals["bm_gr"]
371 + ts_vals["bm_dv"]
372 + ts_vals["bm_dr"]
373 )
374
375 for out in [
376 "bm_gv",
377 "bm_gr",
378 "bm_dv",
379 "bm_dr",
380 "age_gv",
381 "age_gr",
382 "age_dv",
383 "age_dr",
384 "pgro",
385 "gro",
386 "i_bm",
387 "h_bm",
388 "c_bm",
389 "env",
390 "lai",
391 "aet",
392 "wr",
393 "sen_gv",
394 "sen_gr",
395 "abs_dv",
396 "abs_dr",
397 "omd_gv",
398 "omd_gr",
399 ]:
400 outputs_dict[out].append(ts_vals[out])
401
402 return outputs_dict