Getting Started with dplR
In dplR: Dendrochronology Program Library in R

knitr::opts_chunk$set(
  message = FALSE, warning = FALSE,
  fig.width = 7, fig.height = 4.5, fig.align = "center",
  dev = "png", dpi = 96
)

This vignette walks through the core of a tree-ring analysis in dplR: reading ring widths, describing them, detrending, building a chronology, and checking the crossdating. It uses only dplR and base R, and every example runs on data that ship with the package. For a longer treatment, with more on each step and on signal-free chronologies and time-series analysis, see Learning to Love dplR.

library(dplR)

Reading ring widths

Ring widths are read with read.rwl(), which handles Tucson (decadal), compact, Heidelberg, CSV-style spreadsheets and TRiDaS files. It guesses the format, but naming it is safer:

dat <- read.rwl("mysite.rwl", format = "tucson")

The result is an rwl object: a data frame with one column per series and one row per year, with the years as row names and NA where a series has no ring. Here we use co021, 35 Douglas-fir series from Mesa Verde, Colorado, which was read from the International Tree-Ring Data Bank this way.

data(co021)
class(co021)
dim(co021)
co021[1:5, 1:4]

Describing the data

rwl.report() gives an overview: number of series, span, mean length, mean interseries correlation, and any missing rings or suspicious values.

rwl.report(co021)

summary() returns per-series statistics, and plot() shows where each series sits in time.

head(summary(co021))
plot(co021, plot.type = "spag")

Detrending

Raw ring widths carry an age-related growth trend and differences in mean growth between trees. Detrending fits a curve to each series and divides the ring widths by it, giving dimensionless indices with a mean of about one. detrend() does this for every series; detrend.series() does one and plots the fit, which is a good way to choose a method.

x <- co021[, "641114"]
names(x) <- rownames(co021)
x.rwi <- detrend.series(x, method = c("Spline", "ModNegExp"),
                        make.plot = TRUE)

Here we use a cubic smoothing spline with a 50% frequency cutoff at two-thirds of each series' length (the default for "Spline").

co021.rwi <- detrend(co021, method = "Spline")

The result is an rwi object with the same shape as the rwl, which records how it was made. summary() describes the indices as a collection. rbar.eff is the mean interseries correlation and EPS the expressed population signal; an EPS above about 0.85 is the usual rule of thumb for a chronology that represents the population. These come from rwi.stats(), and here count each core as its own tree; pass ids to group cores by tree. The summary also correlates each series with the mean of the others and lists any that do not fit.

summary(co021.rwi)

plot() with plot.type = "image" shows every index at once, years across and series down, brown below 1 and green above. Vertical stripes are years the trees agree on, which is the signal a chronology is built from.

plot(co021.rwi, plot.type = "image")

The same plot is a quick check on the detrending. Dividing each series by its mean leaves the age trend in, and it shows as a green run at the start of nearly every series.

plot(detrend(co021, method = "Mean"), plot.type = "image")

Building a chronology

chron() averages the indices by year, using Tukey's biweight robust mean by default. The result is a crn object: the chronology and the number of series behind each year.

co021.crn <- chron(co021.rwi)
tail(co021.crn)
plot(co021.crn, add.spline = TRUE, nyrs = 32)

The early part of this chronology rests on few series (fewer than five before 1234), so check samp.depth before trusting a given year.

Checking the crossdating

Crossdating assigns each ring its exact calendar year. dplR does not replace visual crossdating, but it can check it statistically, the way COFECHA does. To see what an error looks like, we plant one: the 1500 ring of series 641143 is deleted, so every ring before 1500 is now dated one year too late.

dat <- co021
x <- dat[, "641143"]
names(x) <- rownames(dat)
dat[, "641143"] <- delete.ring(x, year = 1500)

corr.rwl.seg() correlates overlapping segments of each series against a master built from all the other series. Segments that do not correlate significantly are flagged, and with lag.max it also finds the shift at which each segment correlates best.

crs <- corr.rwl.seg(dat, seg.length = 50, pcrit = 0.01,
                    lag.max = 10, label.cex = 0.7)

Every tested segment of 641143 that ends before 1500 correlates best at a lag of -1, and every later one at a lag of 0. In dplR, as in COFECHA, a negative lag means the series is missing a ring, so this points to a missing ring near 1500.

ok <- !is.na(crs$best.lag["641143", ])
crs$best.lag["641143", ok]

ccf.series.rwl() looks at one series in more detail, plotting the full cross-correlation with the master for each segment.

ccf <- ccf.series.rwl(rwl = dat[, colnames(dat) != "641143"],
                      series = dat[, "641143"],
                      series.yrs = as.numeric(rownames(dat)),
                      seg.length = 50, bin.floor = 50)

xdate.report() puts all of this into a COFECHA-style report. It lists the flagged segments with their best lag and the gain in correlation, and records the settings and file checksum so the report can be reproduced. Printing it shows the report in the console.

rpt <- xdate.report(dat, title = "co021 with a planted fault")
rpt

## Split the Markdown version of the report into its sections so the
## vignette can show a couple of them as tables. Headings become bold
## labels, which keeps them out of the table of contents and out of
## pandoc's section structure.
md <- format(rpt, type = "markdown")
md <- sub("^#{1,2} (.*)$", "**\\1**", md)
sec <- cumsum(grepl("^\\*\\*", md))
sections <- split(md, sec)
names(sections) <- sapply(sections, function(s) gsub("\\*", "", s[1]))
show.section <- function(name) {
    cat(sections[[name]], sep = "\n")
    cat("\n")
}

The report opens with a summary of the collection:

show.section("Summary")

The flagged segments section is where to start. All nine flags are on 641143, all are B flags at a lag of -1, and all end before 1500:

show.section("Flagged segments")

The same table is in the report as a data frame, rpt$flagged, ready for further work:

rpt$flagged[, c("series", "from", "to", "flag", "best.lag", "gain")]

A B flag means the segment correlates better at another lag than at its dated position. An A flag, which does not appear here, means the dated position is the segment's best but falls short of significance. The report also has a correlation table for every segment of every series, descriptive statistics, the output of rwl.check(), and notes on how each number was computed:

The full report

wzxhzdk:21

write.xdate.report() saves the report as text, Markdown or HTML:

write.xdate.report(rpt, "co021-report.html")

The statistics say where to look; the wood says what happened. Having found the likely missing ring, you would go back to the sample, and once confirmed, fix the series with insert.ring().

Where to go next

Learning to Love dplR covers these steps and more in depth.
vignette("math-dplR") gives the mathematics behind the smoothing splines, detrending and other functions.
help(package = "dplR") lists every function. Good next stops are ?detrend, ?chron, ?ssf (signal-free chronologies), ?corr.rwl.seg and ?xdate.report.
dplR is part of openDendro.