R/one_day.R

In IhsanKhaliq/ascotraceR: Simulate the Spread of Ascochyta Blight in Chickpea

#' Simulates Ascochyta rabiei spore dispersal for a single day increment
#'
#' @param i_date A character string or POSIXct formatted string indicating an
#'   iteration date of the model. Preferably in \acronym{ISO8601} format
#'   (YYYY-MM-DD), \emph{e.g.}, \dQuote{2020-04-26}.
#' @param daily_vals A `list` of model variables which have been calculated for
#'   days prior to the `i_date` containing ...
#' @param weather_dat A `data.table` of weather observations which includes the
#'   query date, `i_date`, containing...
#' @param max_gp A numeric double value, a function of
#'   `max_gp_lim * (1 - exp(-0.138629 * seeding_rate))`.
#' @param spore_interception_parameter A value indicating the probability of a
#'   spore landing on a susceptible growing point.
#' @param spores_per_gp_per_wet_hour The number of spores produced per
#'   sporulating growing point each wet hour. Also known as the `spore_rate`.
#'   The value is dependent on the susceptibility of the host genotype.
#' @param splash_cauchy_parameter A parameter used in the cauchy distribution and
#'   describes the median distance spores travel due to rain splashes. default:
#'   `0.5`
#' @param wind_cauchy_multiplier A scaling parameter to estimate a cauchy distribution
#'  which resembles the possible distances a conidium travels due to wind dispersal.
#'  default: `0.015`
#' @param daily_rain_threshold Minimum cumulative rainfall required in a day to
#'  allow hourly spore spread events. See also `hourly_rain_threshold`
#'  Default: `2`
#' @param hourly_rain_threshold Minimum rainfall in an hour to trigger a spore
#'  spread event in the same hour (Assuming daily_rain_threshold is already met).
#'  Default: `0.1`
#' @param susceptible_days number of days a new growing point remains susceptible
#'  to infection from ascochyta spores. Default: `5`
#' @param rainfall_multiplier logical values will turn on or off rainfall multiplier
#'  default method. The default method increases the number of spores spread per
#'  growing point if the rainfall in the spore spread event hour is greater than
#'  one. Numeric values will scale the number of spores spread per growing point
#'  against the volume of rainfall in the hour. Default: `FALSE`
#'
#' @return A `list` detailing daily data for the day `i_date` generated by the
#'   model including:
#'   * paddock - an 'x' 'y' `data.table` containing
#'         * `x`
#'         * `y`
#'         * `new_gp`
#'         * `susceptible_gp`
#'         * `exposed_gp`
#'         * `infectious_gp`
#'   * `i_day` - iteration day
#'   * cumulative daily weather data
#'         * `cdd` - cumulative degree days
#'         * `cwh` - cumulative wet hours
#'         * `cr` - cumulative rainfall in mm
#'         * `gp_standard` - standard growing points assuming growth is not
#'           impeded by infection
#'   * `new_gp` - new growing points produced in the last 24 hours
#'   * `infected_coords` - a `data.table` of only infectious growing point
#'     coordinates
#'   * `new_infections` - a `data.table` of newly infected growing points
#'   * `exposed_gp` - a `data.table` of exposed growing points in the latent
#'     period phase of infection
#'
#' @keywords internal
#' @noRd
one_day <- function(i_date,
                    daily_vals,
                    weather_dat,
                    gp_rr,
                    max_gp,
                    spore_interception_parameter,
                    spores_per_gp_per_wet_hour,
                    splash_cauchy_parameter = 0.5,
                    wind_cauchy_multiplier = 0.015,
                    daily_rain_threshold = 2,
                    hourly_rain_threshold = 0.1,
                    susceptible_days = 5,
                    rainfall_multiplier = FALSE) {
  times <- temp <- rain <- new_gp <- infectious_gp <-
    cdd_at_infection <- susceptible_gp <- exposed_gp <-
    spores_per_packet <- NULL

  # expand time to be hourly
  i_time <- rep(i_date, 24) + lubridate::dhours(0:23)

  # subset weather data by day
  weather_day <-
    weather_dat[times %in% i_time, ]

  # obtain summary weather for i_day
  i_mean_air_temp <- mean(weather_day$temp)
  i_wet_hours <- weather_day[, sum(rain > 0,na.rm = TRUE)]
  i_rainfall <- sum(weather_day[, rain], na.rm = TRUE)

  # Start building a list of values for 'i'
  # NOTE: I may add this to after `make_some_infective`
  daily_vals$cdd <- sum(daily_vals$cdd, i_mean_air_temp)
  daily_vals$cwh <- sum(daily_vals$cwh, i_wet_hours)
  daily_vals$cr <- sum(daily_vals$cr, i_rainfall)

  # max new growing points are multiplied by five as growing points remain
  # susceptible for five days.
  max_interception_probability <-
    interception_probability(target_density = susceptible_days *
                               max(daily_vals[["paddock"]][, new_gp]),
                             k = spore_interception_parameter)

  # need to make a copy of the data.table otherwise it will modify all
  # data.tables in the following functions
  daily_vals$paddock <- copy(daily_vals$paddock)
  if (any(is.na(daily_vals[["paddock"]][, infectious_gp]))) {
    stop(
      call. = FALSE,
      "`NA` values in daily_vals[['paddock']][,infectious_gp]")
  }

  # Don't spread spores if there are no infected coordinates
  if (nrow(daily_vals[["infected_coords"]]) > 0) {
    # Spread spores and infect plants
    # Update growing points for paddock coordinates
    if (i_rainfall > daily_rain_threshold) {
      exposed_dt <-
        rbindlist(
          lapply(
            seq_len(nrow(weather_day[rain >= hourly_rain_threshold,])),
            FUN = spores_each_wet_hour,
            weather_hourly = weather_day[rain >= hourly_rain_threshold,],
            paddock = daily_vals$paddock,
            max_interception_probability = max_interception_probability,
            spore_interception_parameter = spore_interception_parameter,
            spores_per_gp_per_wet_hour = spores_per_gp_per_wet_hour,
            rainfall_multiplier = rainfall_multiplier,
            splash_cauchy_parameter = splash_cauchy_parameter,
            wind_cauchy_multiplier = wind_cauchy_multiplier
          )
        )
      exposed_dt[, cdd_at_infection := daily_vals[["cdd"]]]


      daily_vals$exposed_gps <-
        rbind(daily_vals$exposed_gps,
              exposed_dt)

    }

    # exposed gps which have undergone latent period are moved to sporulating gps
    daily_vals <- make_some_infective(daily_vals = daily_vals,
                                      latent_period = 150)


    daily_vals$paddock <- merge(
      x = daily_vals$paddock[, exposed_gp := NULL],
      y = daily_vals$exposed_gps[, list(exposed_gp = sum(spores_per_packet)), by = c("x", "y")],
      by = c("x", "y"),
      all.x = TRUE
    )
    # merge creates NA values, update these to zeros
    daily_vals$paddock[is.na(exposed_gp), exposed_gp := 0]
    setcolorder(daily_vals$paddock,
                c("x",
                  "y",
                  "new_gp",
                  "susceptible_gp",
                  "exposed_gp",
                  "infectious_gp"
                  ))

    # update infected coordinates
    daily_vals$infected_coords <-
      daily_vals$paddock[infectious_gp > 0,
                              c("x", "y")]

  }

  # Grow Plants
  # this code represents mathematica function `growth`; `updateRefUninfectiveGrowingPoints`

  # `updateGrowingPointsAllInfectiveElements`
  # Update Growing points for non-infected coords for time i
  daily_vals$new_gp <-
    calc_new_gp(
      current_growing_points = daily_vals$gp_standard,
      gp_rr = gp_rr,
      max_gp = max_gp,
      mean_air_temp = i_mean_air_temp
    )

  # this might be quicker if there was no fifelse statement
  daily_vals$paddock[, new_gp := fcase(
    susceptible_gp == 0,
    as.integer(0),
    susceptible_gp == daily_vals$gp_standard,
    daily_vals$new_gp,
    susceptible_gp < daily_vals$gp_standard,
    calc_new_gp(
      current_growing_points = susceptible_gp,
      gp_rr = gp_rr,
      max_gp = max_gp,
      mean_air_temp = i_mean_air_temp
    )
  )]

  daily_vals$gp_standard <- sum(daily_vals$gp_standard, daily_vals$new_gp)

  daily_vals[["paddock"]][, susceptible_gp := susceptible_gp + new_gp]

  return(daily_vals[])
}