R/markov_in_modifiers.R

In MartinHoldrege/precipr: Manipulating and examining precipitation intensity

# modifiers for markov .in files generated by rSOILWAT2

# adjust mkv_covar.in -----------------------------------------------------
# note this function could be updated to just work with the coefficients
# file (and calculate expected )

#' adjust weekly correction temp correction factors
#'
#' @param mkv_woy dataframe of weekly temp datagenerated by
#'     rSOILWAT2::dbW_weatherData_to_dataframe
#' @param data dataframe of daily temp and precip
#' @param mean_mult multiplier that mean precipitation event size is to
#' be multiplied by
#'
#' @return dataframe with same columns as mkv_woy, except the CF have been
#' adjusted (if mean_mult not 1)
#' @export
#'
#' @examples
#' data <-data.frame(rSOILWAT2::dbW_weatherData_to_dataframe(rSOILWAT2::weatherData))
#' mkv_woy <- rSOILWAT2::dbW_estimate_WGen_coefs(data)[[1]]
#' head(adjust_mkv_woy(mkv_woy, data))
#' head(adjust_mkv_woy(mkv_woy, data, mean_mult = 2))
adjust_mkv_woy <- function(mkv_woy, data, mean_mult = 1) {
  stopifnot(is.data.frame(mkv_woy),
            is.data.frame(data),
            # can't have NAs
            all(!is.na(data$Tmax_C)),
            all(!is.na(data$Tmin_C)))


  # creating hypothetical new mean temps to use to adjust the correction factors
  wk_smry <- wk_summary_wgen(data, mean_mult = mean_mult)

  # adjust correction factors so that mean temp stays the same when
  # the number of dry days changes
  out <- mkv_woy
  out$CF_Tmax_dry <- wk_smry$Tmax_dry - wk_smry$Tmax_new
  out$CF_Tmin_dry <- wk_smry$Tmin_dry - wk_smry$Tmin_new
  out$CF_Tmax_wet <- wk_smry$Tmax_wet - wk_smry$Tmax_new
  out$CF_Tmin_wet <- wk_smry$Tmin_wet - wk_smry$Tmin_new
  out
}

# probability of rain  ----------------------------------------------------

#' probability of rain each day based on mkv_doy
#'
#' @param mkv_doy dataframe of daily markov ppt transition probabilities
#' @param adjust_for_truncnorm logical, whether to adjust the number of
#' expected wet days based on fact that the weather generator uses a truncated
#' normal distribution (i.e. some 'wet' days then become dry)
#'
#' @return numeric vector giving probability of rain each day (non conditional)
#' @export
#'
#' @examples
#' data <-data.frame(rSOILWAT2::dbW_weatherData_to_dataframe(rSOILWAT2::weatherData))
#' mkv_doy <- rSOILWAT2::dbW_estimate_WGen_coefs(data)[[2]]
#' calc_p_W(mkv_doy = mkv_doy)
calc_p_W <- function(mkv_doy, adjust_for_truncnorm = FALSE) {

  n <- nrow(mkv_doy)
  stopifnot(n == 366, # code here meant to adjust for leap years
            is.logical(adjust_for_truncnorm))

  # probability of wet day on day 0,
  # taking the P_W_D, of DOY 366, which
  # will be a slight under estimate of P_W, so run loop again below
  p_W_day0 <- mkv_doy$p_W_D[n]
  p_W <- rep(NA, n) # probability doy is wet (unconditional)


  if(adjust_for_truncnorm) {
    # if wet day prob will actually be > 0 ppt
    prob_gt0 <- with(mkv_doy, pnorm(0, mean = PPT_avg, sd = PPT_sd,
                                    lower.tail = FALSE))
  }


  # probability that doy 1 is wet
  p_W[1] <- mkv_doy$p_W_W[1]*p_W_day0 + mkv_doy$p_W_D[1]*(1 - p_W_day0)

  if (adjust_for_truncnorm) {
    # not worrying that p_W_day0 wasn't adjusted for prob_gt0 b/ loop will
    # be run again anyways
    p_W[1] <- p_W[1]*prob_gt0[1]
    for (i in 2:n) {
      p_W[i] <- mkv_doy$p_W_W[i]*p_W[i-1] + mkv_doy$p_W_D[i]*(1 - p_W[i-1])
      # now adjust for truncated normal issue for next round in loop
      p_W[i] <- p_W[i]*prob_gt0[i]
    }
  } else {
    for (i in 2:n) {
      p_W[i] <- mkv_doy$p_W_W[i]*p_W[i-1] + mkv_doy$p_W_D[i]*(1 - p_W[i-1])
    }
  }


  # running again with p_W doy 366/365 now known better
  # corrects the first few days of the year (4 days with one test dataset),
  # after that values are identical to previous loop
  # taking weighted mean of p_W for doy 366 and 365 to account for leap years

  if (adjust_for_truncnorm) {
    # multiplying p_W by probability that drawn ppt will actually be >0
    p_W_day0 <- 0.25*p_W[n]*prob_gt0[n] +  0.75*p_W[n - 1]*prob_gt0[n - 1]

    # calc based on previous days wet status
    p_W[1] <- mkv_doy$p_W_W[1]*p_W_day0 + mkv_doy$p_W_D[1]*(1 - p_W_day0)
    p_W[1] <- p_W[1]*prob_gt0[1] # adjust for todays chance of drawing positive number
    for (i in 2:n) {
      p_W[i] <- mkv_doy$p_W_W[i]*p_W[i-1] + mkv_doy$p_W_D[i]*(1 - p_W[i-1])

      p_W[i] <- p_W[i]*prob_gt0[i]
    }
  } else {
    p_W_day0 <- 0.25*p_W[n] +  0.75*p_W[n - 1]
    p_W[1] <- mkv_doy$p_W_W[1]*p_W_day0 + mkv_doy$p_W_D[1]*(1 - p_W_day0)
    for (i in 2:n) {
      p_W[i] <- mkv_doy$p_W_W[i]*p_W[i-1] + mkv_doy$p_W_D[i]*(1 - p_W[i-1])
    }
  }
  p_W
}


# expected annual ppt -----------------------------------------------------


#' @title  expected annual ppt based on markov coefficients
#'
#' @description Calculates probability of rain each day and from that calculates
#' expected value of annual precip. Accounts for the fact that when wgen draws
#' negative value from normal distribution that day is converted to a dry day--
#' this is still in testing, and currently I have removed this feature (ie not
#' accounting for truncated normal???)
#' The problem with this function currently is that while observed/simulated
#'  and expected
#' ppt is normally very close, when the sd is increased, the mismatch grows
#'
#' @param coeffs --list that contains mkv_doy element
#' @param adjust_for_truncnorm logical, whether to adjust the number of
#' expected wet days based on fact that the weather generator uses a truncated
#' normal distribution (i.e. some 'wet' days then become dry)
#'
#' @return expected value of annual precipitation
#'
#' @export
#' @examples
#' data <-data.frame(rSOILWAT2::dbW_weatherData_to_dataframe(rSOILWAT2::weatherData))
#' coeffs <- rSOILWAT2::dbW_estimate_WGen_coefs(data)
#' ex <- expected_ppt(coeffs, adjust_for_truncnorm = TRUE)
#' ex
expected_ppt <- function(coeffs, adjust_for_truncnorm = FALSE) {
  # probability of ppt each day (based on wet and dry probs)
  p_W <- calc_p_W(coeffs$mkv_doy,
                  adjust_for_truncnorm = adjust_for_truncnorm)

  # Further investigation will be needed to determine
  # whether this approach appropriately deals with the changes in probabilities
  # and expected values when mean event sizes are increased
  if (adjust_for_truncnorm) {
    # expected value from truncated normal
    trunc_exp <- truncnorm::etruncnorm(a = 0, b = Inf,
                                      mean = coeffs$mkv_doy$PPT_avg,
                                      sd = coeffs$mkv_doy$PPT_sd)
    # Nan value returned when 0 mean and 0 sd
    trunc_exp <- ifelse(is.nan(trunc_exp), 0, trunc_exp)
    ex <- p_W*trunc_exp
  } else {
    ex <- p_W*coeffs$mkv_doy$PPT_a
  }

  # weighted sum of annual ppt during leap yrs and normal years
  out <- 0.25 * sum(ex) + 0.75*sum(ex[-366])
  out
}


# expected n wet ----------------------------------------------------------

#' expected number of wet days
#'
#' @param coeffs list that contains mkv_doy element
#' @param adjust_for_truncnorm logical, whether to adjust the number of
#' expected wet days based on fact that the weather generator uses a truncated
#' normal distribution (i.e. some 'wet' days then become dry)
#'
#' @return expected number of wet days per year
#' @export
expected_nwet <- function(coeffs, adjust_for_truncnorm = FALSE) {
  stopifnot("mkv_doy" %in% names(coeffs))
  p_W <- calc_p_W(coeffs$mkv_doy,
                  adjust_for_truncnorm = adjust_for_truncnorm)

  # if(adjust_for_truncnorm) {
  #   # want to get probability of getting value >0, ie true wet
  #   prob_gt0 <- with(coeffs$mkv_doy,
  #                    pnorm(0, mean = PPT_avg, sd = PPT_sd, lower.tail = FALSE))
  #   # note this adjustment doesn't actually appropriately account for
  #   # the non independence of subsequent probabilities in the markov chain.
  #   p_W <- p_W*prob_gt0
  # }

  # weighted mean accounting for leap yrs
  n_wet <- 0.25*sum(p_W) + 0.75*sum(p_W[-366])
  n_wet
}

# adjust mkv_prob.in ------------------------------------------------------

#' adjust prob.in markov file to increase mean event size
#'
#' @param mkv_doy prob.in file from rSOILWAT2
#' @param mean_mult how much mean event size should be multiplied by
#' @param adjust_sd logical, whether to adjust the PPT_sd for a give DOY.
#' This option is used so that a draw from a normal distribution with a new
#' (larger) mean is no more/less likely draw a negative value
#'
#' @return mkv_doy file with frequency and mean event size adjusted
#' to maintain same total precipitation
#'
#' @export
#'
#' @examples
#' data <-data.frame(rSOILWAT2::dbW_weatherData_to_dataframe(rSOILWAT2::weatherData))
#' mkv_doy <- rSOILWAT2::dbW_estimate_WGen_coefs(data)[[2]]
#' adjust_mkv_doy(mkv_doy, mean_mult = 2)
adjust_mkv_doy <- function(mkv_doy,
                           mean_mult = 1,
                           adjust_sd = FALSE) {
  stopifnot(is.data.frame(mkv_doy),
            nrow(mkv_doy) == 366)

  out <- mkv_doy

  # calculate unconditional probabilities of ppt
  p_W <- calc_p_W(mkv_doy, adjust_for_truncnorm = TRUE)


  # reducing probabilities of ppt
  # out$p_W_D is not changed, because the probability of ppt the previous
  # day already was reduced by 1/mult

  n <- nrow(mkv_doy)
  # p_W for day 0 (ie last day of year, account for leap years)
  p_W_doy0 <- 0.25*p_W[n] +  0.75*p_W[n - 1]

  # the i-1 (ie previous doy's) p_w values
  p_W_offset <- c(p_W_doy0, p_W[-n])

  # this formula is being used so the the new p_W is equal to the old p_W/mean mult
  # (i.e by appropriately dividing p_W_D)
  W_D_mult <- mean_mult*(1 - p_W_offset/mean_mult)/(1 - p_W_offset)
  out$p_W_D <- mkv_doy$p_W_D/W_D_mult

  # increasing mean event size
  out$PPT_avg <- out$PPT_avg*mean_mult

  # now calc_pW(mkv_doy)/calc_p_W(out) should equal mean_mult

  if (adjust_sd) {
    # adjusting sd to account for the fact that annual ppt totals increase
    # when you increase the standard deviation (for a fix mean), because negative
    # draws are replaced with zero. This is my attempt to keep the expected value
    # of the turncated normal distribution

    # here using un-adjusted mean
    z_vec <- (0 - mkv_doy$PPT_avg)/mkv_doy$PPT_sd # z score of 0 value
    z_vec[!is.finite(z_vec)] <- NA

    # here using new adjusted mean
    out$PPT_sd <- (0-out$PPT_avg)/z_vec

    # in case above calculate led to NAs, replace with original sd
    bad_sd <- is.na(out$PPT_sd) | !is.finite(out$PPT_sd)
    out$PPT_sd[bad_sd] <- mkv_doy$PPT_sd[bad_sd]
  }
  out
}

#' adjust the two elements of the list for the markov wx generator
#'
#' @description The purpose of this function is to adjust the wgen input files
#' so that simulated weather has an increased mean event size,
#' decreased ppt frequency, and no change in mean temperature. Also no change
#' in total precipitation
#'
#' @param coeffs list created by rSOILWAT2::dbW_estimate_WGen_coefs
#' @param data  daily wather data (dataframe)
#' @param mean_mult what to multiply mean event sizes by
#' @param adjust_sd logical, whether to adjust the PPT_sd for a give DOY.
#' This option is used so that a draw from a normal distribution with a new
#' (larger) mean is no more/less likely draw a negative value
#'
#' @return a list with same elements as
#' @export
#'
#' @examples
#' # data <- precipr::wdata119 # this a more 'difficult' dataset to try
#' data <- data.frame(rSOILWAT2::dbW_weatherData_to_dataframe(rSOILWAT2::weatherData))
#' coeffs <- rSOILWAT2::dbW_estimate_WGen_coefs(data, imputation_type = "mean")
#' adjust_coeffs(coeffs, data, mean_mult = 2)
adjust_coeffs <- function(coeffs, data, mean_mult = 1, adjust_sd = TRUE) {
  stopifnot(
    is.list(coeffs),
    names(coeffs) == c("mkv_woy", "mkv_doy"),
    is.data.frame(data)
  )

  out <- coeffs

  # adjust temperature, to keep mean the same when fewer dry days
  out$mkv_woy <- adjust_mkv_woy(mkv_woy = coeffs$mkv_woy,
                                 data = data,
                                 mean_mult = mean_mult)

   # adjust mean event size and frequency
  out$mkv_doy <- adjust_mkv_doy(mkv_doy = coeffs$mkv_doy,
                                mean_mult = mean_mult,
                                adjust_sd = adjust_sd)

  out
}