In EarthSystemDiagnostics/proxysnr: Separate Signal and Noise in Climate Proxy Records

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Overview

This vignette documents all the steps needed to obtain the main results presented in Münch and Laepple (2018) along with the plotting of the respective main figures.

Obtaining the spectral data

library(proxysnr)
library(RColorBrewer)

Produce the main spectral results for the DML1, DML2 and WAIS oxygen isotope data sets (i.e. the raw and corrected signal and noise spectra):

DWS <- WrapSpectralResults(
    dml1 = dml$dml1, dml2 = dml$dml2, wais = wais,
    diffusion = diffusion.tf,
    time.uncertainty = time.uncertainty.tf,
    df.log = c(0.15, 0.15, 0.1))

This function is only a wrapper for the main package functions ArraySpectra and SeparateSpectra which calls the two functions for all the data sets that are specified as input to WrapSpectralResults. The function ArraySpectra is used, for a specific data set, to calculate all individual spectra, the corresponding mean spectrum and the spectrum of the stacked record (thus, of the average isotope record in the time domain); SeparateSpectra is used to obtain the raw and corrected signal and noise spectra for this data set.

The data sets analysed in the paper are provided along proxysnr in the variables dml and wais; see ?dml and ?wais for details on these data sets.

The applied transfer functions to correct for the loss in high-frequency spectral power by the effects of diffusion and time uncertainty are provided in the variables diffusion.tf and time.uncertainty.tf, respectively; see ?diffusion.tf, ?time.uncertainty.tf, and the vignette vignette(topic = "calculate-transfer-functions", package = "proxysnr") for details on obtaining these functions.

The output from WrapSpectralResults is a list of the spectral results for each of the data sets providing the estimated signal, noise and signal-to-noise ratio spectra (i) without any correction applied ("raw"), (ii) for only applying the diffusion correction ("corr.diff.only"), (iii) for only applying the time uncertainty correction ("corr.t.unc.only"), and (iv) for applying both corrections ("corr.full"):

ls.str(DWS)

For all or only for some of the data sets, you can omit both or one of the two transfer functions from the call to WrapSpectralResults in which case only the raw, i.e. uncorrected, or only the partially corrected signal and noise spectra are returned for these data sets (see ?WrapSpectralResults for details).

Plot the main figures

Plot DML1 isotope array spectra (Figure 1)

PlotArraySpectra(ArraySpectra(dml$dml1, df.log = 0.12),
                 f.cutoff = DWS$dml1$corr.full$f.cutoff[2])

Plot DML and WAIS signal and noise spectra (Figure 2)

proxysnr:::muench_laepple_fig02(DWS, f.cut = TRUE)

Plot frequency dependence of signal-to-noise ratios (Figure 3)

Obtain the final signal-to-noise ratio "spectra" by combining both DML data sets and applying additional logarithmic smoothing for visual purposes:

SNR <- proxysnr:::PublicationSNR(DWS$dml1$corr.full, DWS$dml2$corr.full,
                                 DWS$wais$corr.full)

Plot the figure:

PlotSNR(SNR, f.cut = TRUE,
        names = c("DML", "WAIS"), col = c("black", "dodgerblue4"))

Plot estimated correlation with common signal (Figure 4)

Calculate the estimated correlation of a stacked isotope record with the underlying common signal as a function of records averaged and the temporal averaging period (i.e. resolution) of the records:

# for the DMl data
crl1 <- StackCorrelation(SNR$dml, N = 1 : 20,
                         freq.cut.lower = 1 / 100,
                         freq.cut.upper = SNR$dml$f.cutoff[2])

# for the WAIS data
crl2 <- StackCorrelation(SNR$wais, N = 1 : 20,
                         freq.cut.lower = 1 / 100,
                         freq.cut.upper = SNR$wais$f.cutoff[2])

Specify a function to create a colour palette for the contour plots:

palette <- colorRampPalette(rev(RColorBrewer::brewer.pal(10, "RdYlBu")))

Plot the figure for DML:

PlotStackCorrelation(freq = crl1$freq, correlation = crl1$correlation,
                     col.pal = palette, label = expression(bold("a.")~"DML"),
                     ylim = c(NA, log(50)))

Plot the figure for WAIS:

PlotStackCorrelation(freq = crl2$freq, correlation = crl2$correlation,
                     col.pal = palette, label = expression(bold("b.")~"WAIS"),
                     ylim = c(NA, log(50)))

Plot comparison of DML and Trench noise spectra (Figure 5)

Obtain the noise spectra from the trench oxygen isotope data (the data are supplied in the variable t15; see ?t15 for details on the data set):

TNS <- proxysnr:::TrenchNoise()

Plot the comparison of the noise spectra:

proxysnr:::muench_laepple_fig05(SNR, TNS, f.cut = TRUE)

Literature cited

Münch, T. and Laepple, T.: What climate signal is contained in decadal- to centennial-scale isotope variations from Antarctic ice cores?, Clim. Past, 14, 2053-2070, doi: 10.5194/cp-14-2053-2018, 2018.

EarthSystemDiagnostics/proxysnr documentation built on Sept. 15, 2024, 7:47 a.m.