R/sedproxy-interface.R

In EarthSystemDiagnostics/spectraluncertainty:

#' Weights from Phi_c and tau_p
#'
#' @param n.wts Number of required weights
#' @inheritParams ProxyErrorSpectrum
#'
#' @return A named vector of 0/1 weights
#' @export
#'
#' @examples
#' n.wts <- 12
#' wts <- MonthlyWtsFromPhi_ctau_p(phi_c = 0, tau_p = 0.1/n.wts, n.wts = n.wts)
#' wts.col <- rep("Black", n.wts)
#' wts.col[wts == 1] <- "Red"
#' plot(1:n.wts, cos(seq(pi, 3*pi - 2*pi/n.wts, length.out = n.wts)),
#'      col = wts.col, pch = 16, type = "b")
MonthlyWtsFromPhi_ctau_p <- function(phi_c, tau_p, n.wts = 12){

  # phi_c = 0 -> abundance centred at maximum of seasonal cycle

  # phi_c = pi -> abundance centred at minimum of seasonal cycle

  # phi_c = +-pi/2 -> abundance centred at one of the zero crossings
  # of the seasonal cycle

  stopifnot(phi_c <= pi & phi_c >= -pi)
  stopifnot(tau_p > 0 & tau_p <= 1)

  mid.wt <- n.wts/2 + phi_c * n.wts / (2*pi) +1
  mid.wt[mid.wt == 0] <- 1

  n.1s <- round(n.wts * tau_p)
  n.1s[n.1s == 0] <- 1

  first.1 <- ceiling(mid.wt - n.1s / 2)
  last.1 <- first.1 + (n.1s-1)

  i <- first.1:last.1

  i <- ifelse(i%%n.wts == 0, n.wts, i%%n.wts)

  wts <- rep(0, n.wts)
  wts[i] <- 1
  names(wts) <- 1:n.wts
  return(wts)
}

#' PlotProxySeasonality
#'
#' @param res resolution
#' @inheritParams ProxyErrorSpectrum
#'
#' @return
#' @export
#'
#' @examples
#' PlotProxySeasonality(phi_c = 0, tau_p = 1/12, res = 365)
#' PlotProxySeasonality(phi_c = 0, tau_p = 6/12, res = 365)
#' PlotProxySeasonality(phi_c = -pi, tau_p = 3/12, res = 365)
#' PlotProxySeasonality(phi_c = -pi, tau_p = 1/12, res = 12)
#' PlotProxySeasonality(phi_c = pi, tau_p = 2/12, res = 12)
PlotProxySeasonality <- function(phi_c, tau_p, res = 365){
  
  n <- res
  sc <- cos(seq(pi, 3*pi - 2*pi/n, length.out = n))
  
  wts <- MonthlyWtsFromPhi_ctau_p(phi_c = phi_c, tau_p = tau_p, n.wts = res)
  
  i <- 1 + (((1:n)-(1)) %% (res))
  
  wts <- wts[i]
  
  dat <- data.frame(x = 1:n, y = sc, wts = as.character(wts))
  p <- ggplot(dat, aes(x = x, y = y, colour = wts, group = NA)) +
    geom_path(size = 1) +
    scale_colour_manual(values = c("1" = "Darkgreen", "0" = "Grey"),
                        guide = FALSE ) +
    labs(x = "1 year", y = "Seasonal cycle") +
    theme_classic() +
    theme(axis.text.x = element_blank())
  
  if (res <= 24) {
    p <- p + geom_point()
  }
  p
}

#' Get Corresponding Sedproxy Parameters From Proxy Error Spectrum Parameters
#'
#' @param spec.pars
#'
#' @return A list of parameters for \code{ClimToProxyClim} from package \code{sedproxy}
#' @export
#'
#' @examples
#' SedproxyParsFromSpecPars(spectraluncertainty::GetSpecPars("Mg_Ca"))
SedproxyParsFromSpecPars <- function(spec.pars){

  n.bd <- 3
  bio.depth <- 10

  end_buffer <-  max(n.bd*spec.pars$tau_b, spec.pars$delta_t/2) + spec.pars$tau_s
  start_buffer <- spec.pars$tau_b + spec.pars$tau_s

  spec.pars.2 <- spec.pars
  spec.pars.2$T <- spec.pars$T + 2 * end_buffer

  if (spec.pars$clim.spec.fun == "ClimPowerFunction"){
    clim.in <- SimClimMat(spec.pars.2)$cl
  } else {clim.in <- NA}

  sedproxy.pars <- with(spec.pars, list(
    clim.signal = clim.in,
    timepoints = seq(0, T, delta_t) + start_buffer,
    plot.sig.res = tau_r,
    bio.depth = bio.depth,
    sed.acc.rate = (1000 * bio.depth) / tau_b,
    calibration.type = "identity",
    noise.type = "additive",
    habitat.weights = MonthlyWtsFromPhi_ctau_p(phi_c = phi_c,
                                               tau_p = tau_p),
    layer.width = tau_s / (tau_b / bio.depth),
    n.bd = n.bd,
    n.samples = N,
    n.replicates = 1,
    sigma.meas = sigma.meas,
    sigma.ind = sigma.ind
  ))

  return(sedproxy.pars)
}


#' Get Error Components from a Sedproxy PFM
#'
#' @param PFM Output from sedproxy::ClimToProxyClim
#'
#' @return A dataframe of proxy error components.
#' @export
#'
#' @examples
#' sed.pars <- SedproxyParsFromSpecPars(spectraluncertainty::GetSpecPars("Mg_Ca"))
#' PFM <- do.call(sedproxy::ClimToProxyClim, sed.pars)
#' err.comps <- GetErrorComponents(PFM)
GetErrorComponents <- function(PFM){

  if (class(PFM) == "sedproxy.pfm"){
    PFM <- PFM$simulated.proxy
  }

  err <- mutate(PFM,
                Climate = clim.timepoints.ssr,
                Bioturbation = proxy.bt - Climate,
                Seasonal.unc. = mean(proxy.bt.sb - proxy.bt),
                Orbital.unc. = (proxy.bt.sb - proxy.bt),
                Aliasing.Y = proxy.bt.sb.sampY - proxy.bt.sb,
                Aliasing.M = proxy.bt.sb.sampYM - proxy.bt.sb.sampY,
                Aliasing.YM = Aliasing.M + Aliasing.Y,
                Measurement = simulated.proxy - proxy.bt.sb.sampYM,
                Total.error = simulated.proxy - Climate)

  err <- select(err,
                -clim.signal.ann)
  return(err)
}


#' Get Spectra of Sedproxy PFM Error Components
#'
#' @param err A data.frame of error components. Output from GetErrorComponents.
#' @param res The resolution of the timeseries'
#'
#' @return A dataframe of error component power spectra
#' @export
#'
#' @examples
#' spec.pars <- spectraluncertainty::GetSpecPars("Mg_Ca")
#' sed.pars <- SedproxyParsFromSpecPars(spec.pars)
#' PFM <- do.call(sedproxy::ClimToProxyClim, sed.pars)
#' err.comps <- GetErrorComponents(PFM)
#' GetSpecFromErrComponents(err.comps, spec.pars$delta_t)
GetSpecFromErrComponents <- function(err, res){

  GetSpecDF <- function(df, res){

    spec <- spectrum(ts(df$Value, 1, frequency = 1/res),
                     detrend = FALSE, taper = 0, plot = F, fast = FALSE)

    data.frame(freq = spec$freq, spec = spec$spec)
  }

  spec <- err %>%
    select(timepoints, Climate, Bioturbation, Seasonal.unc., Orbital.unc.,
           Aliasing.Y, Aliasing.M, Aliasing.YM, Measurement, Total.error) %>%
    gather(stage, Value, -timepoints) %>%
    group_by(stage) %>%
    do({
      GetSpecDF(., res = res)
    }) %>%
    group_by(stage, freq) %>%
    summarise(spec = mean(spec))#%>%

  spec <- spec %>%
    ungroup() %>%
    mutate(stage = factor(stage, ordered = TRUE,
                          levels = c("Climate", "Bioturbation", "Seasonal.unc.", "Orbital.unc.",
                                     "Aliasing.Y", "Aliasing.M", "Aliasing.YM", "Measurement",
                                     "Total.error")),
           Component = "Error spectrum",
           period = 1/freq) %>%
    # Do the integration
    group_by(Component, stage) %>%
    mutate(Integral = (2*(cumsum(spec) * diff(freq)[1]))) %>%
    ungroup() %>%
    mutate(Type = "Numerical")
  return(spec)
}