In EarthSystemDiagnostics/spectraluncertainty:

knitr::opts_chunk$set(echo=TRUE, message = FALSE, warning = FALSE, cache = TRUE,
               fig.width = 7, fig.pos = "H", dpi = 300, autodep = TRUE)

#detach("package:spectraluncertainty", unload = TRUE)
#detach("package:sedproxy", unload = TRUE)
library(spectraluncertainty)
library(sedproxy)
library(tidyverse)

Relative seasonal cycle vs stochastic climate white noise

muc.info <- data.frame(
  Core = "SO213-84-2/MUC",
  Latitude = -45.1245,
  Longitude = 174.5865
)

spec.pars.muc <- GetSpecPars("Mg_Ca")

# PSD Climate
clim.spec.muc <- ModelSpectrum(freq = NULL, latitude = muc.info$Latitude,
                           variable = "T_deg_Mg_Ca", beta = 1)

p.clim.spec <- clim.spec.muc %>% 
  as.data.frame() %>% 
  ggplot(aes(x = freq, y = spec)) +
  geom_line() + 
  geom_segment(aes(x = exp(log(0.03) - 1), xend = exp(log(0.03)+1), y = 1, yend = 1),
               colour = "Black",
               arrow = arrow(length = unit(0.015, "npc"), ends = "both", type = "closed")) +
  annotate("text", x = exp(log(0.03)-1), y = 0.8, label = "Theoretical") + 
  annotate("text", x = exp(log(0.03)+1), y = 0.8, label = "Empirical") + 
  scale_x_continuous("Frequency", trans = "log10",
                     breaks = 1/c(2, 5, 10, 50, 100, 500, 1000),
                     labels = function(n) format(n, drop0trailing = T)) + 
  scale_y_continuous("Power spectral density", trans = "log10") +
  geom_vline(xintercept = 3/100, linetype = 4) +
  annotation_logticks() +
  theme_bw() +
  theme(panel.grid = element_blank())

p.clim.spec

breit.lons <- unique(breitkreuz.amp$longitude)
nearest.lon <- breit.lons[which.min(abs(muc.info$Longitude - breit.lons))]

breit.lats <- unique(breitkreuz.amp$latitude)
nearest.lat <- breit.lats[which.min(abs(muc.info$Latitude - breit.lats))]

seasonal.amp <- AmpFromLocation(longitude = muc.info$Longitude,
                                latitude = muc.info$Latitude,
                                proxy.type = "degC",
                                depth.upr = 0, depth.lwr = -50)

timescale.pars.muc <- crossing(
  T = 1e04 + 100, 
  delta_t = 100,
  s = exp(seq(log(1), log(500), length.out = 10)),
  N =1
)

orbital.pars.muc <- RelativeAmplitudeModulation(latitude = muc.info$Latitude,
                                            maxTimeKYear = 23,
                                            minTimeKYear = 1,
                                            bPlot = FALSE)

spec.pars.lst <- timescale.pars.muc %>% 
  plyr::dlply(., "s", function(x) {
    GetSpecPars(proxy.type = "T_deg_Mg_Ca",
            delta_t = x$delta_t,
            tau_r = x$delta_t,
            T = x$T,

            # Bioturbation
            tau_b = 1000 * 10 /x$s,
            tau_s = 1000 * 0/x$s,

            N = x$N,

            sig.sq_c = seasonal.amp$sig.sq_c,
            sigma.ind = spec.pars.muc$sigma.ind,
            sigma.meas = spec.pars.muc$sigma.meas,

            tau_p = 1,
            phi_c = 0,

            sig.sq_a = orbital.pars.muc$sig.sq_a,

            clim.spec.fun = "ModelSpectrum",
            clim.spec.fun.args = list(latitude = muc.info$Latitude, beta = 1,
                           variable = "T_deg_Mg_Ca")
            )
  })

spec.obj.lst <- plyr::llply(spec.pars.lst, function(x) do.call(ProxyErrorSpectrum, x))

spec.obj.df <- plyr::ldply(spec.obj.lst, .id = "s") %>% 
  mutate(s = as.numeric(as.character(s)))

var.obj.df <- spec.obj.df %>% 
 # mutate(s = as.character(s)) %>% 
  group_by(s) %>% 
  do({
    select(., -s) %>% 
    IntegrateErrorSpectra(err.spec = ., delta_t = spec.pars.muc$delta_t) 
  }) %>% 
  filter(smoothed.resolution == spec.pars.muc$delta_t) %>% 
  mutate(var.ratio = Aliasing.seasonal / Aliasing.stochastic)


var.obj.df %>% 
  ggplot(aes(x = 1000*10/s,
             #x = s,
             y = var.ratio)) +
  geom_line() +
  expand_limits(y = 0) +
  labs(x = "tau_b [years]") +
  #scale_x_continuous(expression("Sedimentation rate [cm kyr"^"-1"*"]"), breaks = c(0, 1, seq(10, 50, 10))) +
  scale_y_continuous("Variance ratio") +
  theme_bw() +
  theme(panel.grid.minor = element_blank())

var.obj.df %>% 
  ggplot(aes(#x = 1000*10/s,
             x = s,
             y = var.ratio)) +
  geom_line() +
  expand_limits(y = 0) +
  #labs(x = "tau_b [years]") +
  scale_x_continuous(expression("Sedimentation rate [cm kyr"^"-1"*"]"), breaks = c(seq(0, 500, 100))) +
  scale_y_continuous("Variance ratio") +
  theme_bw() +
  theme(panel.grid.minor = element_blank())

EarthSystemDiagnostics/spectraluncertainty documentation built on Jan. 8, 2020, 9:43 a.m.