R/gen_month.R

In walkerjeffd/weathergen: Daily weather generator for performing a climate stress test.

#' Generate monthly climate timeseries using ARMA model
#'
#' @param x.reg.yr data.frame with regional average annual climate timeseries (YEAR, PRCP, TMAX, TMIN, TAVG)
#' @param x.loc.mon data.frame with annual climate timeseries at selected location (YEAR, MONTH, PRCP, TMAX, TMIN, TAVG)
#' @param n.iter number of iterations for arima simulation
#' @param n.year number of years for arima simulation
#' @export
#' @examples
#' gen_month_arma(x.reg.yr=clim.reg.yr, x.loc.mon=clim.loc.mon, n.iter=50, n.year=nrow(clim.reg.yr))
#'
gen_month_arma <- function(x.reg.yr, x.loc.mon, n.iter, n.year) {
  # fit arma model to regional average precipitation timeseries
  models <- list(PRCP=forecast::auto.arima(x.reg.yr$PRCP,
                                           max.p=2,max.q=2,max.P=0,max.Q=0,
                                           stationary=TRUE))
  # generate n.iter arima simulations of length n.year years
  sim.prcp.yr <- mc_arimas(models=models, n=n.year, n.iter=n.iter)

  # empty matrix of the sampled years with dims: <year, sim iteration>
  sampled.years <- array(NA, c(n.year, n.iter))

  # empty matrix of monthly simulations with dims: <month/year, variable, sim iteration>
  sim.mon <- array(NA, c(n.year*12, 8, n.iter))

  # loop through simulation iterations
  for (samp in 1:n.iter) {
    # loop through years
    for (j in 1:n.year) {
      # get total annual precipitation of current simulation
      target.prcp <- sim.prcp.yr$x[j, samp]

      # compute distances between simulated timeseries and observed regional precip for top 8 nearest neighbors
      distances <- dplyr::mutate(x.reg.yr,
                                 DISTANCE=sqrt((target.prcp - x.reg.yr$PRCP)^2),
                                 IN_TOP_8=rank(DISTANCE) %in% seq(1,8))
      distances <- dplyr::filter(distances, IN_TOP_8)
      distances <- dplyr::mutate(distances, WEIGHT=(1/DISTANCE)^2,
                                            WEIGHT=WEIGHT/sum(WEIGHT))

      # sample year from top 8 nearest neighbors using sampling weights
      sampled.years[j,samp] <- sample(distances$YEAR,
                                      size=1,
                                      prob=distances$WEIGHT,
                                      replace=FALSE)

      # copy monthly climate timeseries for location and sampled year to final generated timeseries
      row.idx <- ((1+12*(j-1)):(12+12*(j-1)))
      mtrx <- dplyr::filter(x.loc.mon, YEAR==sampled.years[j, samp])
      mtrx <- dplyr::mutate(mtrx, YEAR_SIM=j, YEAR_OBS=YEAR, PRCP_OBS_YEAR=target.prcp)
      mtrx <- dplyr::select(mtrx, YEAR_SIM, MONTH, PRCP, TMAX, TMIN, TAVG, YEAR_OBS, PRCP_OBS_YEAR)
      mtrx <- data.matrix(mtrx)
      sim.mon[row.idx, , samp] <- mtrx
    }
  }
  sim.mon.list <- apply(sim.mon, 3, function(x) {
    df <- as.data.frame(x)
    names(df) <- c('YEAR_SIM', 'MONTH', 'PRCP', 'TMAX', 'TMIN', 'TAVG', 'YEAR_OBS', 'PRCP_OBS_YEAR')
    return(df)
  })
  return(list(models=models, sim.prcp.yr=sim.prcp.yr, sim.mon=sim.mon.list))
}

#' Add climate trends to a single generated monthly climate timeseries
#'
#' @param sim.gen data.frame of a single monthly climate timeseries generated by gen_month_arma()
#' @param temp.factor additive temperature trend factor
#' @param prcp.factor multiplicative precipitation trend factor
#' @export
#' @examples
#' sim.gen <- gen_month_arma(x.reg.yr=clim.reg.yr, x.loc.mon=clim.loc.mon, n.iter=50, n.year=nrow(clim.reg.yr))
#' add_trend_mon(sim.gen=sim.gen$sim.mon[[1]], temp.factor=2, prcp.factor=1.5)
#'
add_trend_mon <- function(sim.gen, temp.factor, prcp.factor) {
  trend.factors <- list(PRCP=prcp.factor, # multiplicative trend factor
                        TEMP=temp.factor) # additive trend factor

  # create linear trend lines
  trend.lines <- list(PRCP=seq(1, trend.factors[['PRCP']], length.out=nrow(sim.gen)),
                      TEMP=seq(0, trend.factors[['TEMP']], length.out=nrow(sim.gen)))

  sim.trend <- dplyr::mutate(sim.gen,
                             PRCP.TREND=PRCP * trend.lines[['PRCP']],
                             TMAX.TREND=TMAX + trend.lines[['TEMP']],
                             TMIN.TREND=TMIN + trend.lines[['TEMP']],
                             TAVG.TREND=TAVG + trend.lines[['TEMP']])
  return(list(trend.factors=trend.factors,
              df=sim.trend))
}

#' Add multiple climate trends to generated monthly climate timeseries
#'
#' @param sim.gen.list list of data.frames containing multiple monthly climate timeseries generated by gen_month_arma()
#' @param temp.factors vector of additive temperature trend factors
#' @param prcp.factors vector multiplicative precipitation trend factors
#' @export
#' @examples
#' sim.gen <- gen_month_arma(x.reg.yr=clim.reg.yr, x.loc.mon=clim.loc.mon, n.iter=50, n.year=nrow(clim.reg.yr))
#' trends_mon(sim.mon.list=sim.gen$sim.mon, temp.factors=seq(0, 2, by=0.5), prcp.factors=seq(0.8, 1.2, by=0.1))
#'
trends_mon <- function(sim.mon.list, temp.factors, prcp.factors) {
  trends.factors <- expand.grid(TEMP=temp.factors,
                                PRCP=prcp.factors)

  sim.mon.trends <- apply(trends.factors, 1, function(factors) {
    temp.factor <- factors[["TEMP"]]
    prcp.factor <- factors[["PRCP"]]

    # select random simulation trial
    trial <- sample(length(sim.mon.list), 1)
    sim.mon <- sim.mon.list[[trial]]

    sim.trend <- add_trend_mon(sim.gen=sim.mon,
                               temp.factor=temp.factor,
                               prcp.factor=prcp.factor)
    sim.trend$trial <- trial
    return(sim.trend)
  })

  return(sim.mon.trends)
}

#' Generate monthly climate timeseries using ARMA model with trend factors
#'
#' @param x.reg.yr data.frame with regional average annual climate timeseries (YEAR, PRCP, TMAX, TMIN, TAVG)
#' @param x.loc.mon data.frame with annual climate timeseries at selected location (YEAR, MONTH, PRCP, TMAX, TMIN, TAVG)
#' @param n.iter number of iterations for arima simulation
#' @param n.year number of years for arima simulation
#' @param temp.factors vector of additive temperature trend factors
#' @param prcp.factors vector of multiplicative precipitation trend factors
#' @export
#' @examples
#' gen_month_arma_trends(x.reg.yr=clim.reg.yr, x.loc.mon=clim.loc.mon, n.iter=50, n.year=nrow(clim.reg.yr), sim.mon.list=sim.gen$sim.mon, temp.factor=2, prcp.factor=1.5)
#'
gen_month_arma_trend <- function(x.reg.yr, x.loc.mon, n.iter, n.year, temp.factors, prcp.factors) {
  sim.gen <- gen_month_arma(x.reg.yr=x.reg.yr, x.loc.mon=x.loc.mon, n.iter=n.iter, n.year=n.year)
  sim.gen.trends <- trends_mon(sim.mon.list=sim.gen$sim.mon, temp.factors=temp.factors, prcp.factors=prcp.factors)
  return(sim.gen.trends)
}