-
Notifications
You must be signed in to change notification settings - Fork 0
/
trace_asco.R
339 lines (300 loc) · 13.2 KB
/
trace_asco.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
#' Simulates the spread of Ascochyta blight in a chickpea field
#'
#' `trace_asco` simulates the spatiotemporal development of Ascochyta blight in a chickpea field
#' over a growing season
#'
#' @param weather weather data, recorded by a local weather station, over a chickpea
#' growing season for the model operation
#' @param paddock_length length of a paddock in metres (y)
#' @param paddock_width width of a paddock in metres (x)
#' @param sowing_date a character string of a date value indicating sowing
#' date of chickpea seed and the start of the ascotraceR model. Preferably
#' in \acronym{ISO8601} format (YYYY-MM-DD), \emph{e.g.} \dQuote{2020-04-26}.
#' Assumes there is sufficient soil moisture to induce germination and start the
#' crop growing season.
#' @param harvest_date a character string of a date value indicating harvest date of
#' chickpea crop, which is also the last day to run the ascotraceR model. Preferably in
#' \acronym{ISO8601} format (YYYY-MM-DD), \emph{e.g.} \dQuote{2020-04-26}.
#' @param seeding_rate indicate the rate at which chickpea seed is sown per
#' square metre. Defaults to \code{40}
#' @param gp_rr refers to rate of increase in chickpea growing points
#' per degree Celsius per day. Defaults to \code{0.0065}
#' @param max_gp_lim Maximum number of chickpea growing points (meristems) allowed
#' per square meter. Defaults to \code{15000}.
#' @param max_new_gp Maximum number of new chickpea growing points (meristems)
#' which develop per day, per square meter. Defaults to \code{350}.
#' @param primary_infection_foci it refers to the inoculated quadrat
#' located at the centre of the paddock from where disease spreads
#' Defaults to \code{"centre"}
#' @param primary_infection_intensity The intensity of the starting epidemic as
#' described by the number of number of sporulating growing points.
#' @param latent_period_cdd latent period in cumulative degree days (sum of
#' daily temperature means) is the period between infection and production of
#' lesions on susceptible growing points. Defaults to \code{200}
#' @param initial_infection refers to initial or primary infection on seedlings,
#' resulting in the production of infected growing points
#' @param time_zone refers to time in Coordinated Universal Time (UTC)
#' @param spores_per_gp_per_wet_hour Number of spores produced per sporulating growing point each wet hour.
#' Also known as the 'spore_rate'. Value is dependent on the susceptibility of the host genotype.
#' @param n_foci only relevant when primary_infection_foci = "random" and notes the number
#' of primary_infection_foci at initial infection.
#'
#' @return a x y `data.frame` simulating the spread of Ascochyta blight in a
#' chickpea paddock
#' @export
#'
#' @examples
#' First weather data needs to be imported and formatted with `format_weather`
#' Billa_Billa <-
#' readr::read_csv(system.file("extdata",
#' "2020_Billa_Billa_weather.csv", package = "ascotraceR")) %>%
#' # convert wind direction recorded as text to degrees
#' dplyr::mutate(
#' wind_direction = dplyr::case_when(
#' max_wind_direction == "N" ~ "0",
#' max_wind_direction == "NbE" ~ "11.25",
#' max_wind_direction == "NNE" ~ "22.5",
#' max_wind_direction == "NEbN" ~ "33.75",
#' max_wind_direction == "NE" ~ "45",
#' max_wind_direction == "NEbE" ~ "56.25",
#' max_wind_direction == "ENE" ~ "67.5",
#' max_wind_direction == "EbN" ~ "73.5",
#' max_wind_direction == "E" ~ "90",
#' max_wind_direction == "EbS" ~ "101.2",
#' max_wind_direction == "ESE" ~ "112.5",
#' max_wind_direction == "SEbE" ~ "123.8",
#' max_wind_direction == "SE" ~ "135.1",
#' max_wind_direction == "SEbS" ~ "146.3",
#' max_wind_direction == "SSE" ~ "157.6",
#' max_wind_direction == "SbE" ~ "168.8",
#' max_wind_direction == "S" ~ "180",
#' max_wind_direction == "SbW" ~ "191.2",
#' max_wind_direction == "SSW" ~ "202.5",
#' max_wind_direction == "SWbS" ~ "213.8",
#' max_wind_direction == "SW" ~ "225",
#' max_wind_direction == "SWbW" ~ "236.2",
#' max_wind_direction == "WSW" ~ "247.5",
#' max_wind_direction == "WbS" ~ "258.8",
#' max_wind_direction == "W" ~ "270",
#' max_wind_direction == "WbN" ~ "281.2",
#' max_wind_direction == "WNW" ~ "292.5",
#' max_wind_direction == "NWbW" ~ "303.8",
#' max_wind_direction == "NW" ~ "315",
#' max_wind_direction == "NWbN" ~ "326.2",
#' max_wind_direction == "NNW" ~ "337.5",
#' max_wind_direction == "NbW" ~ "348.8",
#' TRUE ~ max_wind_direction
#' )
#' ) %>%
#' dplyr::mutate(wind_direction = as.numeric(wind_direction))
#'
#' station_data <-
#' system.file("extdata", "stat_dat.csv", package = "ascotraceR")
#'
#' weather_dat <- format_weather(x = Billa_Billa,
#' POSIXct_time = "local_time",
#' temp = "mean_daily_temperature",
#' ws = "avg_wind_speed (km/h)",
#' wd_sd = "wd_sd",
#' rain = "rain",
#' wd = "wind_direction",
#' station = "location",
#' time_zone = "Australia/Brisbane",
#' lonlat_file = station_data
#' )
#'
#'
#' traced <- trace_asco(
#' weather = weather_dat,
#' paddock_length = 20,
#' paddock_width = 20,
#' initial_infection = "2020-07-16",
#' sowing_date = as.POSIXct("2020-06-04"),
#' harvest_date = as.POSIXct("2020-11-27") + lubridate::ddays(30),
#' time_zone = "Australia/Perth",
#' primary_infection_foci = "center")
#'
trace_asco <- function(weather,
paddock_length,
paddock_width,
sowing_date,
harvest_date,
initial_infection,
seeding_rate = 40,
gp_rr = 0.0065,
max_gp_lim = 15000,
max_new_gp = 350,
latent_period_cdd = 200,
time_zone = "UTC",
primary_infection_foci = "random",
primary_infection_intensity = 1,
n_foci = 1,
spores_per_gp_per_wet_hour = 0.22){
x <- y <- sp_gp <- NULL
# check date inputs for validity -----------------------------------------
.vali_date <- function(x) {
tryCatch(
# try to parse the date format using lubridate
x <- lubridate::parse_date_time(x,
c(
"Ymd",
"dmY",
"mdY",
"BdY",
"Bdy",
"bdY",
"bdy"
)),
warning = function(c) {
stop(call. = FALSE,
"\n",
x,
" is not a valid entry for date. Enter as YYYY-MM-DD.\n")
}
)
return(x)
}
if (primary_infection_intensity > seeding_rate) {
stop(
"primary_infection_intensity exceeds the number of starting growing points - 'seeding_rate': ",
seeding_rate
)
}
# convert times to POSIXct -----------------------------------------------
initial_infection <-
lubridate::ymd(.vali_date(initial_infection), tz = time_zone) + lubridate::dhours(0)
sowing_date <-
lubridate::ymd(.vali_date(sowing_date), tz = time_zone) + lubridate::dhours(0)
harvest_date <-
lubridate::ymd(.vali_date(harvest_date), tz = time_zone) + lubridate::dhours(23)
# check epidemic start is after sowing date
if(initial_infection <= sowing_date){
stop("initial_infection occurs prior to sowing_date\n
please submit an initial_infection date which occurs after crop_sowing")
}
# makePaddock equivalent
paddock <- as.data.table(expand.grid(x = 1:paddock_width,
y = 1:paddock_length))
# sample a paddock location randomly if a starting foci is not given
if (is.data.table(primary_infection_foci) &
all(c("x", "y") %in% colnames(primary_infection_foci))) {
# Skip the rest of the tests
} else{
if (class(primary_infection_foci) == "character") {
if (primary_infection_foci == "random") {
primary_infection_foci <-
paddock[sample(seq_len(nrow(paddock)),
size = n_foci,
replace = TRUE),
c("x", "y")]
} else{
if (primary_infection_foci == "center") {
primary_infection_foci <-
paddock[x == as.integer(round(paddock_width / 2)) &
y == as.integer(round(paddock_length / 2)),
c("x", "y")]
}else{
stop("primary_infection_foci input not recognised")
}
}
} else{
if (is.vector(primary_infection_foci)) {
if (length(primary_infection_foci) != 2 |
is.numeric(primary_infection_foci) == FALSE) {
stop("primary_infection_foci should be supplied as a numeric vector of length two")
}
primary_infection_foci <-
as.data.table(as.list(primary_infection_foci))
setnames(x = primary_infection_foci,
old = c("V1", "V2"),
new = c("x", "y"))
}
if (is.data.table(primary_infection_foci) == FALSE &
is.data.frame(primary_infection_foci)) {
setDT(primary_infection_foci)
if (all(c("x", "y") %in% colnames(primary_infection_foci)) == FALSE) {
stop("primary_infection_foci data.table needs colnames 'x' and 'y'")
}
}
}
}
infected_rows <- which_paddock_row(paddock = paddock,
query = primary_infection_foci)
if(ncol(primary_infection_foci) == 2){
primary_infection_foci[,sp_gp := primary_infection_intensity]
}
# define paddock variables at time 1
#need to update so can assign a data.table of things primary infection foci!!!!!!!!!!!!!!!
paddock[, c(
"new_gp", # Change in the number of growing points since last iteration
"noninfected_gp",
"infected_gp",
"sporulating_gp" # replacing InfectiveElementList
) :=
list(
seeding_rate,
seeding_rate,
0,
0
)]
paddock[infected_rows, "noninfected_gp" :=
seeding_rate - primary_infection_foci[,3]]
paddock[infected_rows, "sporulating_gp" :=
primary_infection_foci[,3]]
# calculate additional parameters
spore_interception_parameter <-
0.00006 * (max_gp_lim/max_new_gp)
# define max_gp
max_gp <- max_gp_lim * (1 - exp(-0.138629 * seeding_rate))
# Notes: as area is 1m x 1m many computation in the mathematica
# code are redundant because they are being multiplied by 1.
# I will reduce the number of objects containing the same value,
# Below is a list of Mathematica values consolidated into 1
#
# refUninfectiveGPs <- minGrowingPoints <- seeding_rate
daily_vals_list <- list(
list(
paddock = paddock, # data.table each row is a 1 x 1m coordinate
i_date = sowing_date, # day of the simulation (iterator)
i_day = 1,
day = lubridate::yday(sowing_date), # day of the year
cdd = 0, # cumulative degree days
cwh = 0, # cumulative wet hours
cr = 0, # cumulative rainfall
gp_standard = seeding_rate, # standard number of growing points for 1m^2 if not inhibited by infection (refUninfectiveGrowingPoints)
new_gp = seeding_rate, # new number of growing points for current iteration (refNewGrowingPoints)
infected_coords = primary_infection_foci, # data.frame
newly_infected = data.table(x = numeric(),
y = numeric(),
spores_per_packet = numeric(),
cdd_at_infection = numeric()) # data.table of infected growing points still in latent period and not sporilating (exposed_gp)
)
)
time_increments <- seq(sowing_date,
harvest_date,
by = "days")
daily_vals_list <- rep(daily_vals_list,
length(time_increments)+1)
for(i in seq_len(length(time_increments))){
# update time values for iteration of loop
daily_vals_list[[i]][["i_date"]] <- time_increments[i]
daily_vals_list[[i]][["i_day"]] <- i
daily_vals_list[[i]][["day"]] <- lubridate::yday(time_increments[i])
# currently working on one_day
daily_vals_list[[i + 1]] <- one_day(i_date = time_increments[i],
daily_vals = daily_vals_list[[i]],
weather_dat = weather,
gp_rr = gp_rr,
max_gp = max_gp,
spore_interception_parameter = spore_interception_parameter,
spores_per_gp_per_wet_hour = spores_per_gp_per_wet_hour)
# temporary line of code to test building of daily_vals in loop
#daily_vals_list <- day_out
}
daily_vals_list[[length(daily_vals_list)]][["i_date"]] <-
daily_vals_list[[length(daily_vals_list)]][["i_date"]] + lubridate::ddays(1)
daily_vals_list[[length(daily_vals_list)]][["i_day"]] <- length(daily_vals_list)
daily_vals_list[[length(daily_vals_list)]][["day"]] <-
lubridate::yday(daily_vals_list[[length(daily_vals_list)]][["i_date"]])
return(daily_vals_list)
}