Package {dplR}

Dendrochronology Program Library in R

Description

This package contains functions for performing some standard tree-ring analyses.

Details

Main Functions

read.rwl

reads rwl files

detrend

detrends raw ring widths

chron

builds chronologies

corr.rwl.seg

crossdating function

Author(s)

Andy Bunn andy.bunn@wwu.edu with major additions from Mikko Korpela and other significant contributions from Franco Biondi, Filipe Campelo, Pierre Mérian, Fares Qeadan and Christian Zang. Function redfit is an improved translation of program REDFIT which is original work of Manfred Mudelsee and Michael Schulz. Jacob Cecile contributed a bug fix to detrend.series. Allan Buras came up with the revised formula of glk in dplR >= 1.6.1.

References

Cook, E. R. and Kairiukstis, L. A., editors (1990) Methods of Dendrochronology: Applications in the Environmental Sciences. Springer. ISBN-13: 978-0-7923-0586-6.

Fritts, H. C. (2001) Tree Rings and Climate. Blackburn. ISBN-13: 978-1-930665-39-2.

Subset an rwl Object

Description

Subset a rwl object by series (columns) and years (rows), keeping the class and the provenance record of the data.

Usage

## S3 method for class 'rwl'
x[i, j, drop]

## S3 method for class 'rwl'
subset(x, subset, select, drop = FALSE, ...)

Arguments

Details

Subsetting an rwl object works as it does for a data.frame, with i indexing years and j indexing series. Four things are different.

Dropping series returns the years the remaining series cover. A collection is ragged, and the years at the ends of an rwl object usually belong to a few long series, so dropping some of the series ordinarily leaves years in which nothing that remains was measured. The result runs from the first year to the last year in which any of the selected series has a measurement. Empty years between those two are kept: an rwl object holds one row per year, and an interior year cannot be dropped without breaking the sequence.

This applies to calls that drop series without indexing years – x[, j], x[j], and subset given only select. Two kinds of call leave the years alone. A call that indexes years returns exactly the years it asks for, measured or not: x[i, ], x[i, j], head, tail, common.interval, and subset with its subset argument. A year window is how an rwl object is aligned with something else, such as a climate series or a second collection, and returning fewer years than were asked for would put that alignment out. And a call that keeps every series, such as the reordering x[, order(names(x))] or x[], has dropped nothing to trim for, whatever empty years the object already held.

Two results are not trimmed. A subset holding no measurement in any year is returned whole rather than emptied, and rwl.check reports it. A single series taken with drop = TRUE is a numeric vector, which carries no years at all, so it is returned at full length and stays aligned with the object it came from; use drop = FALSE for a one-series rwl object trimmed to its own years.

subset drops empty series; [ keeps them. Taking years can leave a series with no values at all: in ca533[time(ca533) %in% 1800:1899, ], the series CAM132, CAM152, CAM161 and CAM201 end before 1800 and are left as columns of nothing but NA. An empty column is not a series. chron, rwi.stats and the plots skip it, and summary.rwi, interseries.cor, corr.rwl.seg and detrend either drop it or fail on it.

subset and window.rwl therefore drop any series left with no values, and a message names the series dropped. If no series is left with values, that is an error. Use one of these to take a span of years for analysis.

[ keeps every series it is asked for, empty or not, on purpose. x[i, j] always returns the columns that j names, in that order, so the result stays lined up with anything indexed alongside it: a vector of series names, an ids table from read.ids, or a second object with the same series. If [ dropped columns, those would fall out of step without a word, and code that indexes with [ – dplR's own included – could no longer count on names(x[, j]) being the series it asked for. So row subsets made with [, and head and tail, can return empty series. Drop them with window or subset, or find them with colSums(!is.na(x)) == 0.

Years must stay consecutive. An rwl object holds one row per year in increasing order, and dplR reads the year of a measurement from its row name, so time.rwl, plot.rwl, detrend, chron, rwl.stats and the crossdating functions all take the next row to be the next year. A subset such as x[c(1, 5, 9), ] leaves that untrue. Where row subsetting leaves years that are not consecutive and increasing, because rows were taken from the middle, repeated, or reordered, the result is returned as a plain data.frame, with a warning, and without the provenance record. The measurements are unchanged; what is withdrawn is the claim that they are a set of dated series.

The provenance record is carried and cut to fit. The "dplR.provenance" attribute attached by read.tucson follows the subset. What describes the file – its name, the reader, the header lines – is kept. What describes series – the precision of each, the renames, the interior gaps, the problems found while parsing – is kept for the series that remain, and gaps are dropped where they fall outside the years that remain. mixed.precision is recomputed, since a file measured at two precisions can be subset down to series measured at one. A subset element is added, holding the series in the result, its first and last year, and all.series, which is TRUE when every series the file held is still present. A record that cannot be matched to a series in the result is dropped rather than kept, so an object whose series have been renamed since it was read carries a record that says nothing about them in preference to one that says something wrong.

common.interval cuts an object down to the years its series share, which is a different question from the years they cover between them.

Value

An object of class c("rwl", "data.frame"), or a numeric vector when one column is selected and drop is TRUE, or a data.frame when row subsetting has left the years no longer consecutive. When series are dropped and years are not indexed, the result runs only over the years the remaining series cover. [ returns every series asked for, including any with no values; subset drops those, with a message.

Author(s)

Andy Bunn

See Also

read.rwl, read.tucson, rwl.check, common.interval, time.rwl, subset.data.frame

Examples

library(utils)
data(ca533)

## Series are columns and years are rows. Selecting series drops the years
## that none of them covers: ca533 begins in 626, but these five series do
## not.
ca533.sub <- ca533[, 1:5]
class(ca533.sub)
range(time(ca533))
range(time(ca533.sub))

## Naming years returns those years, whether or not they hold data.
dim(ca533[time(ca533) %in% 1000:1099, 1:5])

## subset() works the same way.
range(time(subset(ca533, select = 1:5)))

## Years are subset by name or by time().
ca533.1800s <- ca533[time(ca533) %in% 1800:1899, ]
range(time(ca533.1800s))

## Four series end before 1800. `[` keeps them as empty columns, so the
## result still has the 34 series asked for; subset() and window() drop
## them and say which.
ncol(ca533.1800s)
names(ca533.1800s)[colSums(!is.na(ca533.1800s)) == 0]
ncol(subset(ca533, time(ca533) %in% 1800:1899))
ncol(window(ca533, 1800, 1899))

## A single series comes back as a vector unless drop is FALSE.
class(ca533[, 1])
class(ca533[, 1, drop = FALSE])

## The provenance record of the data follows the subset.
data(wa082)
gap.series <- attr(wa082, "dplR.provenance")$gaps$series
wa082.sub <- wa082[, gap.series, drop = FALSE]
attr(wa082.sub, "dplR.provenance")$gaps
attr(wa082.sub, "dplR.provenance")$subset$all.series

## Rows taken out of the middle break the one-row-per-year promise that
## dplR relies on. This warns, and gives back a data.frame.
not.rwl <- ca533[c(1, 5, 9), ]
class(not.rwl)

Age-Dependent Spline

Description

Applies an age-dependent smoothing spline to y.

Usage

ads(y, nyrs0 = 50, pos.slope = TRUE)

Arguments

Details

This implements the age-dependent smoothing spline similar to that described by Melvin (2004). In this implementation a cubic smoothing spline (caps) is fit to y with an initial stiffness of nyrs0. For each succesive measurement, the siffness is incremented by that ring index. This results in a spline is nyrs0 flexible at the start of the series and grows progressively stiffer. This can help capture the initial fast growth of a juvinielle tree with a flexible spline that then progresses to a stiffer spline that can better model the constant growth commonly found in mature trees. In its details, the cubic smoothing spline follows the Cook and Peters (1981) spline with a 50% frequency cutoff. See Cook and Kairiukstis (1990) for more information.

The default setting for nyrs0 is 50 years which is approprite for trees with a classic growth model but a value of 10 or 20 might be more appropriate for a Hugershoff-like initial increase in growth. Cook (pers comm) suggests a value of 20 for RCS.

If pos.slope is FALSE, the function will attempt to prevent a positive slope at the end of the series. In some cases when ads is used for detrending, a positive slope can be considered biologically unlikely. In those cases, the user can constrain the positive slope in the spline. This works by calculating the spline and taking the first difference. Then the function finds the last (outside) index where the spline changes slope and fixes the spline values from the that point to the end. Finally the spline is rerun along this constrained curve. See examples for details. The wisdom of constraining the slope in this manner depends very much on expert knowledge of the system.

Value

A filtered vector.

Author(s)

Fortran code provided by Ed Cook. Ported and adapted for dplR by Andy Bunn.

References

Cook, E. R. and Kairiukstis, L. A., editors (1990) Methods of Dendrochronology: Applications in the Environmental Sciences. Springer. ISBN-13: 978-0-7923-0586-6.

Melvin, T. M. (2004) Historical Growth Rates and Changing Climatic Sensitivity of Boreal Conifers. PhD Thesis, Climatic Research Unit, School of Environmental Sciences, University of East Anglia.

Cook, E. R. and Peters, K. (1981) The Smoothing Spline: A New Approach to Standardizing Forest Interior Tree-Ring Width Series for Dendroclimatic Studies. Tree-Ring Bulletin, 41, 45-53.

See Also

caps, detrend

Examples

# fit a curve
data(co021)
aSeries <- na.omit(co021$`641114`)
plot(aSeries,type="l",col="grey50")
lines(ads(y = aSeries),col="blue",lwd=2)

# show an artificial series with a Hugershoff-like curve.
a <- 0.5
b <- 1
g <- 0.1
d <- 0.25

n <- 300
x <- 1:n
y <- I(a*x^b*exp(-g*x)+d)
# add some noise
y <- y + runif(n=length(y),min = 0,max = 0.5)
# Plot with two different splines.
plot(y,type="l",col="grey50")
lines(ads(y,50),col="darkgreen",lwd=2) # bad
lines(ads(y,10),col="darkblue",lwd=2) # good

# now repeat with a positive slope to constrain
y <- I(a*x^b*exp(-g*x)+d)
y[251:300] <- y[251:300] + seq(0,0.25,length.out=50)
y <- y + runif(n=length(y),min = 0,max = 0.5)
plot(y,type="l",col="grey50")
lines(ads(y,10),col="darkgreen",lwd=2) #bad?
lines(ads(y,10,pos.slope=FALSE),col="darkgreen",lwd=2,lty="dashed")

Rothenburg Tree Ring Widths

Description

This data set gives the raw ring widths for Norway spruce Picea abies at Rothenburg ob der Tauber, Bavaria, Germany. There are 20 series from 10 trees. Data set was created using read.rwl and saved to an .rda file using save.

Usage

data(anos1)

Format

A data.frame containing 20 tree-ring series from 10 trees in columns and 98 years in rows. The correct stc mask for use with read.ids is c(5, 2, 1).

References

Zang, C. (2010) Growth reaction of temperate forest tree species to summer drought – a multispecies tree-ring network approach. Ph.D. thesis, Technische Universität München.

Basal Area Increment (bai) Objects

Description

Class "bai" holds basal area increment: the area of each ring, with one series per column and one year per row, as made by bai.in and bai.out. as.bai labels a data.frame or matrix of ring areas as one.

Usage

as.bai(x)

## S3 method for class 'bai'
x[i, j, drop]

## S3 method for class 'bai'
subset(x, subset, select, drop = FALSE, ...)

## S3 method for class 'bai'
time(x, ...)

## S3 method for class 'bai'
summary(object, ...)

## S3 method for class 'bai'
plot(x, plot.type = c("seg", "spag"), ...)

Arguments

Details

A "bai" object has the same shape as an "rwl" object and makes the same promise: the row names are the years, consecutive and increasing. It inherits from neither "rwl" nor "rwi", so that areas, widths and indices can be told apart.

dplR's functions use the class to check what they are given. Those that want ring-width indices, such as chron and rwi.stats, take a "bai" object without a word: a mean of basal area increment by year is a BAI chronology. So do the crossdating functions, the plots, common.interval, and detrend, since fitting a curve to basal area increment is established practice. Functions that want ring widths, such as bai.in, rcs, cms and rwl.report, warn and treat the areas as widths. If the values really are widths, relabel them with as.rwl.

Two attributes travel with it. attr(x, "dplR.bai") says how the areas were made: the function (fun) and whether the distance to pith (d2pith, for bai.in) or the diameters (diam, for bai.out) were given rather than estimated from the widths. attr(x, "dplR.provenance") is the provenance record of the ring widths the areas were made from, if they had one (see read.tucson). Neither is set by as.bai, which only labels the object.

Subsetting works as for "rwl" objects ([.rwl), and window takes a span of years (window.rwl). summary gives the per-series statistics of rwl.stats, in the units of the areas.

Value

as.bai: an object of class c("bai", "data.frame"). x is returned unchanged if it already has that class. time: a numeric vector of years. summary: a data.frame, as from rwl.stats.

Author(s)

Andy Bunn

See Also

bai.in, bai.out, as.rwl, as.rwi, chron

Examples

library(utils)
data(gp.rwl)
data(gp.d2pith)
gp.bai <- bai.in(rwl = gp.rwl, d2pith = gp.d2pith)
class(gp.bai)
str(attr(gp.bai, "dplR.bai"))
head(summary(gp.bai))
## A BAI chronology
gp.bai.crn <- chron(gp.bai)
plot(gp.bai.crn, add.spline = TRUE, nyrs = 32)
plot(gp.bai[, 1:5])

Ring-Width Index (rwi) Objects

Description

Class "rwi" holds ring-width indices: detrended series with one series per column and one year per row, as made by detrend, rcs and cms. as.rwi labels a data.frame or matrix of indices as one.

Usage

as.rwi(x)

## S3 method for class 'rwi'
x[i, j, drop]

## S3 method for class 'rwi'
subset(x, subset, select, drop = FALSE, ...)

## S3 method for class 'rwi'
time(x, ...)

## S3 method for class 'rwi'
summary(object, ids = NULL, pcrit = 0.05, ...)

## S3 method for class 'summary.rwi'
print(x, max.print = 10, ...)

## S3 method for class 'summary.rwi'
as.data.frame(x, ...)

## S3 method for class 'rwi'
plot(x, plot.type = c("spag", "seg", "image"), ...)

Arguments

Details

An "rwi" object has the same shape as an "rwl" object and makes the same promise: the row names are the years, consecutive and increasing. It does not inherit from "rwl", so that indices and ring widths can be told apart.

dplR's functions use the class to check what they are given. Those that want ring widths – detrend, rcs, cms, i.detrend, bai.in, bai.out, pointer, strip.rwl, ssf and rwl.report – warn when given an "rwi" object, and say what goes wrong: detrend on indices, for instance, divides out a growth curve that has already been removed. Those that want indices – chron, chron.ars, chron.stabilized, rwi.stats, rwi.stats.running and sss – warn when given an "rwl" object: rwi.stats(ca533) gives an rbar.eff of 0.350 from the widths against 0.423 from the Spline indices, and nothing about the number says it is wrong. They still take a plain data.frame or matrix without a word, and a "bai" object of basal area increment (see as.bai). The crossdating functions, the plots, common.interval, sgc and rwl.stats take either class quietly. Each warning is only a warning, and the function goes on: if the class is what is wrong – indices read from a file with read.rwl come back as class "rwl" – relabel them with as.rwi, or relabel widths with as.rwl.

Two attributes travel with it. attr(x, "dplR.detrend") is a list that says how the indices were made: the function (fun), the method and its settings, and always difference, which is TRUE when the indices are differences from the fitted curve (centred on 0) rather than ratios (centred on 1). attr(x, "dplR.provenance") is the provenance record of the ring widths the indices were made from, if they had one (see read.tucson). Neither is set by as.rwi, which only labels the object.

Subsetting works as for "rwl" objects ([.rwl): the class and both records are kept, dropping series trims the years that none of the remaining series cover, and a row subset that leaves the years out of order or with holes returns a plain data.frame, with a warning. subset and window.rwi drop series left with no values and name them in a message; [ keeps them, on purpose, so that its columns stay lined up with anything indexed alongside them.

summary describes the indices as a collection. It gives how they were made (from attr(x, "dplR.detrend")), their span and common interval, the collection statistics from rwi.stats (including rbar.eff, EPS and SNR), and for each series its first and last years, length, mean, standard deviation, first-order autocorrelation, and its correlation with the mean of the other series, with the p-value, from interseries.cor with its defaults. Printing it shows the collection and lists the series that do not correlate with the others at pcrit; as.data.frame gives the per-series table. The collection statistics count each series as its own tree unless ids is given, and the print says which. The print ends by pointing to summary with ids (when it was not given), rwi.stats.running (rbar and EPS through time, which usually fall off where sample depth does) and corr.rwl.seg (where in a series the fit breaks down). A series with no values at all is listed and left out of the collection statistics, the correlations and the common interval. With fewer than two series with values there is no collection, and the statistics and correlations are NULL and NA.

plot draws, by default, spag.plot, which draws indices around 1 (or 0 for differences) and not around each series' mean, so the grey line under each series is the value the indices should sit at. plot.type = "image" draws each index as a coloured cell, with years across and series down, earliest-starting at the bottom. Colours are brown below 1 (or 0) and green above it. Each side is scaled separately and clipped at the clip quantile of its departures, since ratio indices cannot fall below 0 but can run well above 2; the key gives the clips. A run of one colour at the start or end of a series is a growth trend the detrending did not remove, and a vertical stripe is a year the series agree on.

Arithmetic on an "rwi" object (e.g. x + 1) returns a plain data.frame, as it does for any data.frame subclass. Use as.rwi on the result if it is still a set of indices.

Value

as.rwi: an object of class c("rwi", "data.frame"). x is returned unchanged if it already has that class. time: a numeric vector of years. summary: an object of class "summary.rwi", a list with how (the dplR.detrend record), n.series, first, last, common (the first and last years in which every series has a value, or NA), stats (from rwi.stats), ids.given, empty (the names of series with no values), series (the per-series data.frame) and pcrit.

Author(s)

Andy Bunn

See Also

detrend, rcs, cms, as.rwl, chron, rwi.stats

Examples

library(utils)
data(ca533)
ca533.rwi <- detrend(ca533, method = "Spline")
class(ca533.rwi)
str(attr(ca533.rwi, "dplR.detrend"))
summary(ca533.rwi)
## Group cores by tree, from the series IDs
summary(ca533.rwi, ids = autoread.ids(ca533))
head(as.data.frame(summary(ca533.rwi)))
plot(ca533.rwi[, 1:5])
data(co021)
## The growth trend "Mean" leaves in: a run of green at the start of
## nearly every series.
plot(detrend(co021, method = "Mean"), plot.type = "image")

as.rwl

Description

Attempts to turn its argument into a rwl object.

Usage

as.rwl(x)

Arguments

Details

This tries to coerce x into class c("rwl","data.frame"). Failable.

Value

An object of class c("rwl", "data.frame") with the series in columns and the years as rows. The series IDs are the column names and the years are the row names.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

Examples

library(graphics)
library(stats)
library(utils)
## Toy
n <- 100
## Make a data.frame that is tree-ring like
base.series <- 0.75 + exp(-0.2 * 1:n)
foo <- data.frame(x1 = base.series + abs(rnorm(n, 0, 0.25)),
                  x2 = base.series + abs(rnorm(n, 0, 0.25)),
                  x3 = base.series + abs(rnorm(n, 0, 0.25)),
                  x4 = base.series + abs(rnorm(n, 0, 0.25)),
                  x5 = base.series + abs(rnorm(n, 0, 0.25)),
                  x6 = base.series + abs(rnorm(n, 0, 0.25)))
# coerce to rwl and use plot and summary methods
foo <- as.rwl(foo)
class(foo)
plot(foo, plot.type="spag")
summary(foo)

Basal Area Increment (Inside Out)

Description

Convert multiple ring-width series to basal area increment (i.e., ring area) going from the pith to the bark.

Usage

bai.in(rwl, d2pith = NULL)

Arguments

Details

This converts ring-width series (mm) to ring-area series (mm squared) (aka basal area increments) based on the distance between the innermost measured ring and the pith of the tree. It is related to bai.out, which calculates each ring area starting from the outside of the tree and working inward. Both methods assume a circular cross section (Biondi 1999). See the references below for further details.

Value

An object of class c("bai", "data.frame") containing the ring areas for each series with the column names, row names and dimensions of rwl. It carries the provenance record of rwl, if it has one, and records how it was made in attr(x, "dplR.bai"): the function (fun) and whether d2pith was given (d2pith). See as.bai for which functions take it.

Note

DendroLab website: https://dendrolaborg.wordpress.com/

Author(s)

Code by Andy Bunn based on work from DendroLab, University of Nevada Reno, USA. Patched and improved by Mikko Korpela.

References

Biondi, F. (1999) Comparing tree-ring chronologies and repeated timber inventories as forest monitoring tools. Ecological Applications, 9(1), 216–227.

Biondi, F. and Qeadan, F. (2008) A theory-driven approach to tree-ring standardization: Defining the biological trend from expected basal area increment. Tree-Ring Research, 64(2), 81–96.

See Also

bai.out, as.bai

Examples

library(graphics)
library(stats)
library(utils)
## Toy
n <- 100
## Make three fake tree-ring series to show that these funcs work on rwl objects
base.series <- 0.75 + exp(-0.2 * 1:n)
rwl <- data.frame(x1 = base.series + abs(rnorm(n, 0, 0.05)),
                  x2 = base.series + abs(rnorm(n, 0, 0.05)),
                  x3 = base.series + abs(rnorm(n, 0, 0.05)))

## The inside out method
foo <- bai.in(rwl = rwl)
## The outside in method
bar <- bai.out(rwl = rwl)

## Identical
head(bar)
head(foo)

## Use gp data
data(gp.rwl)
data(gp.d2pith)
foo <- bai.in(rwl = gp.rwl, d2pith = gp.d2pith)
foo.crn <- chron(foo)
yrs <- time(foo.crn)
plot(yrs, foo.crn[, 1], type = "n",
     xlab = "Year", ylab = expression(mm^2))
lines(yrs, foo.crn[, 1], col = "grey", lty = "dashed")
lines(yrs, caps(foo.crn[, 1], nyrs = 32), col = "red", lwd = 2)

Basal Area Increment (Outside In)

Description

Convert multiple ring-width series to basal area increment (i.e., ring area) going from the bark to the pith.

Usage

bai.out(rwl, diam = NULL)

Arguments

Details

This converts ring-width series (mm) to ring-area series (mm squared) (aka basal area increments) based on the diameter of the tree and the width of each ring moving towards the pith of the tree. It is related to bai.in, which calculates each ring area starting from the inside of the tree and working outward. Both methods assume a circular cross section (Biondi 1999). See the references below for further details.

Value

An object of class c("bai", "data.frame") containing the ring areas for each series with the column names, row names and dimensions of rwl. It carries the provenance record of rwl, if it has one, and records how it was made in attr(x, "dplR.bai"): the function (fun) and whether diam was given (diam). See as.bai for which functions take it.

Note

DendroLab website: https://dendrolaborg.wordpress.com/

Author(s)

Code by Andy Bunn based on work from DendroLab, University of Nevada Reno, USA. Patched and improved by Mikko Korpela.

References

Biondi, F. (1999) Comparing tree-ring chronologies and repeated timber inventories as forest monitoring tools. Ecological Applications, 9(1), 216–227.

Biondi, F. and Qeadan, F. (2008) A theory-driven approach to tree-ring standardization: Defining the biological trend from expected basal area increment. Tree-Ring Research, 64(2), 81–96.

See Also

bai.in, as.bai

Examples

library(graphics)
library(utils)
## Not run: 
library(stats)
## Toy
n <- 100
## Make three fake tree-ring series to show that these funcs work on rwl objects
base.series <- 0.75 + exp(-0.2 * 1:n)
rwl <- data.frame(x1 = base.series + abs(rnorm(n, 0, 0.05)),
                  x2 = base.series + abs(rnorm(n, 0, 0.05)),
                  x3 = base.series + abs(rnorm(n, 0, 0.05)))

## The inside out method
foo <- bai.in(rwl = rwl)
## The outside in method
bar <- bai.out(rwl = rwl)

## Identical
head(bar)
head(foo)

## End(Not run)
## Use gp data
data(gp.rwl)
data(gp.dbh)
## dbh (minus the bark) from cm to mm 
gp.dbh2 <- gp.dbh[, 1:2]
gp.dbh2[, 2] <- (gp.dbh[, 2] - gp.dbh[, 3]) * 10
bar <- bai.out(rwl = gp.rwl, diam = gp.dbh2)
bar.crn <- chron(bar)
yrs <- time(bar.crn)
plot(yrs, bar.crn[, 1], type = "n",
     xlab = "Year", ylab = expression(mm^2))
lines(yrs, bar.crn[, 1], col = "grey", lty = "dashed")
lines(yrs, caps(bar.crn[, 1], nyrs = 32), col = "red", lwd = 2)

Basal Area Increment (Bakker)

Description

Convert multiple ring-width series to basal area increment (i.e., ring area) following the proportional method of Bakker (2005).

Usage

  bakker(rwl, ancillary)

Arguments

Details

This converts ring-width series (mm) to ring-area series (mm squared) (aka basal area increments) based on the diameter of the tree, the missing distance to the pith and the missing number of rings to the pith, following the proportional method for reconstructing historical tree diameters by Bakker (2005). It prevents bai.out transformations from producing negative increments when the sum of all ring widths in a series is larger than DBH/2. It prevents bai.in transformations from producing too small values when the sum of all ring widths in a series is smaller than DBH/2.

Value

A list containing the following objects:

Author(s)

Code by Stefan Klesse. Adapted for dplR by Andy Bunn.

References

Bakker, J.D., 2005. A new, proportional method for reconstructing historical tree diameters. Canadian Journal of Forest Research 35, 2515–2520. https://doi.org/10.1139/x05-136

Examples

  
data(zof.rwl)
data(zof.anc)
zof.bakker <- bakker(rwl = zof.rwl,ancillary = zof.anc)
zof.bai <- zof.bakker$baiBakker_raw

# first series bai
yrs <- time(zof.rwl)
plot(yrs,zof.bai[,1],type="l",
     xlab="Year",
     ylab=expression(BAI~(mm^2)),
     main = colnames(zof.bai)[1])
  

Campito Mountain Tree Ring Widths

Description

This data set gives the raw ring widths for bristlecone pine Pinus longaeva at Campito Mountain in California, USA. There are 34 series. Data set was created using read.rwl and saved to an .rda file using save.

Usage

data(ca533)

Format

A data.frame containing 34 tree-ring series in columns and 1358 years in rows.

Source

International tree-ring data bank, Accessed on 20-April-2021 at https://www.ncei.noaa.gov/pub/data/paleo/treering/measurements/northamerica/usa/ca533.rwl

References

Graybill, D. A. and LaMarche, Jr., V. C. (1983) Campito Mountain Data Set. IGBP PAGES/World Data Center for Paleoclimatology Data Contribution Series 1983-CA533.RWL. NOAA/NCDC Paleoclimatology Program, Boulder, Colorado, USA.

Examples

library(utils)
data(ca533)

## Where the data came from. read.tucson() records what it saw and attaches
## it to what it returns; this data set carries that record.
prov <- attr(ca533, "dplR.provenance")
prov$file          # the archived file it was read from
prov$precision     # the precision each series was measured at
prov$gaps          # interior gaps, if any -- this file has none
prov$events        # problems found while reading -- none here either

## rwl.report() puts the file and precision at the head of its report
rwl.report(ca533)

Twisted Tree Heartrot Hill Standard Chronology

Description

This data set gives the standard chronology for white spruce Picea glauca at Twisted Tree Heartrot Hill in Yukon, Canada. Data set was created using read.crn and saved to an .rda file using save.

Usage

data(cana157)

Format

A data.frame containing the standard chronology in column one and the sample depth in column two. There are 463 years (1530–1992) in the rows.

Source

International tree-ring data bank, Accessed on 20-April-2021 at https://www.ncei.noaa.gov/pub/data/paleo/treering/chronologies/northamerica/canada/cana157.crn

References

Jacoby, G., D’Arrigo, R. and Buckley, B. (1992) Twisted Tree Heartrot Hill Data Set. IGBP PAGES/World Data Center for Paleoclimatology Data Contribution Series 1992-CANA157.CRN. NOAA/NCDC Paleoclimatology Program, Boulder, Colorado, USA.

Cook and Peters Smoothing Spline with User-Specified Rigidity and Frequency Cutoff

Description

Applies a smoothing spline to y with rigidity determined by two parameters: frequency response f at a wavelength of nyrs years.

Usage

caps(y, nyrs = 32, f = 0.5)

Arguments

Details

This applies the classic smoothing spline from Cook and Peters (1981). The rigidity of the spline has a frequency response of 50% at a wavelength of nyrs. The references, of course, have more information.

This function was introduced to dplR in version 1.7.3 and replaces the now defunct ffcsaps. Where ffcsaps was written entirely in R, caps is a wrapper for a Fortran subroutine from Ed Cook's ARSTAN program that is hundreds of times faster.

The default value of nyrs was changed to 32 from 2/3 length of y in version 1.7.8 based on a suggestion from Klesse (2021).

Note: nyrs is passed to the Fortran subroutine as an integer, so a fractional nyrs is truncated. This matters when nyrs is given as a proportion of series length, which is rarely a whole number. The effect on the fitted curve is small but not zero; see the “Mathematical Details” vignette for a worked comparison.

Note: caps stops if there are any NA values in y. A smoothing spline through a series with missing values is not defined, and the caller has to decide what to do about the gap rather than have one chosen silently. Before version 1.8.0 such a call returned a vector of NA, which looked like an answer and was not. See examples, and fill.internal.NA.

Value

A filtered vector.

Note

The mathematical relationship between nyrs, f and the spline's smoothing parameter, together with an empirical check of the frequency response and a demonstration that caps reproduces the deprecated ffcsaps, is given in the vignette: vignette("math-dplR", package = "dplR").

Author(s)

Fortran code provided by Ed Cook and adapted for dplR by Andy Bunn.

References

Cook, E. R. and Kairiukstis, L. A., editors (1990) Methods of Dendrochronology: Applications in the Environmental Sciences. Springer. ISBN-13: 978-0-7923-0586-6.

Cook, E. R. and Peters, K. (1981) The Smoothing Spline: A New Approach to Standardizing Forest Interior Tree-Ring Width Series for Dendroclimatic Studies. Tree-Ring Bulletin, 41, 45-53.

Klesse, S. (2021) Critical Note on the Application of the "Two-Third"" Spline. Dendrochronologia Volume 65: 125786

See Also

ads

Examples

library(graphics)
library(utils)

## Use first series from the Mesa Verde data set
data(co021)
series <- co021[, 1]
series <- series[!is.na(series)]
plot(series, type = "l", ylab = "Ring Width (mm)", col = "grey")
lines(caps(series, nyrs = 10), col = "red", lwd = 2)
lines(caps(series, nyrs = 100), col = "green", lwd = 2)
# The default, a 32-year spline since version 1.7.8
lines(caps(series), col = "blue", lwd = 2)
legend("topright",
       c("Series", "nyrs=10", "nyrs=100", "Default nyrs (32)"),
       fill=c("grey", "red", "green", "blue"))


## caps() needs a series with no missing values. Note that the example
## above drops the NA padding before fitting. A gap inside the measured
## span has to be dealt with deliberately: drop the series, or close the
## gap with fill.internal.NA().
y <- c(NA, NA, rnorm(100))
try(caps(y))
## fit over the measured values instead
head(caps(y[!is.na(y)]))

Cross-Correlation between a Series and a Master Chronology

Description

Computes cross-correlations between a tree-ring series and a master chronology built from a rwl object at user-specified lags and segments.

Usage

ccf.series.rwl(rwl, series, series.yrs = as.numeric(names(series)),
               seg.length = 50, bin.floor = 100, n = NULL,
               nyrs = NULL, prewhiten = TRUE, ar.order.max = NULL,
               biweight = TRUE, pcrit = 0.05,
               lag.max = 5, make.plot = TRUE,
               floor.plus1 = FALSE, series.x = FALSE, ...)

Arguments

Details

This function calculates the cross-correlation function between a tree-ring series and a master chronology built from rwl looking at correlations lagged positively and negatively using ccf at overlapping segments set by seg.length. For instance, with lag.max set to 5, cross-correlations would be calculated at for each segment with the master lagged at k = -5:5 years.

The cross correlations are calculated calling ccf as
ccf(x=master, y=series, lag.max=lag.max, plot=FALSE) if series.x is FALSE and as ccf(x=series, y=master, lag.max=lag.max, plot=FALSE) if series.x is TRUE. This argument was introduced in dplR version 1.7.0. Different users have different expectations about how missing or extra rings are notated. If series.x = FALSE the behavior will be like COFECHA where a missing ring in a series produces a negative lag in the plot rather than a positive lag.

Correlations are calculated for the first segment, then the second segment and so on. Correlations are only calculated for segments with complete overlap with the master chronology.

Each series (including those in the rwl object) is optionally detrended as the residuals from a hanning filter with weight n. The filter is not applied if n is NULL. Detrending can also be done via prewhitening where the residuals of an ar model are added to each series mean. This is the default. The master chronology is computed as the mean of the rwl object using tbrm if biweight is TRUE and rowMeans if not. Note that detrending typically changes the length of the series. E.g., a hanning filter will shorten the series on either end by floor(n/2). The prewhitening default will change the series length based on the ar model fit. The effects of detrending can be seen with series.rwl.plot.

Value

A list containing matrices ccf and bins. Matrix ccf contains the correlations between the series and the master chronology at the lags window given by lag.max. Matrix bins contains the years encapsulated by each bin.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

References

Bunn, A. G. (2010) Statistical and visual crossdating in R using the dplR library. Dendrochronologia, 28(4), 251–258.

See Also

corr.rwl.seg, corr.series.seg, skel.plot, series.rwl.plot

Examples

library(utils)
data(co021)
dat <- co021
## Create a missing ring: delete the 1500 ring of series 641143 and
## drop the original from the master.  Dated from the bark, every ring
## before 1500 now sits one year late.
flagged <- dat$"641143"
names(flagged) <- rownames(dat)
flagged <- delete.ring(flagged, year = 1500)
dat$"641143" <- NULL
ccf.100 <- ccf.series.rwl(rwl = dat, series = flagged, seg.length = 100)
## Not run: 
flagged2 <- co021$"641143"
names(flagged2) <- rownames(dat)
ccf.100.1 <- ccf.series.rwl(rwl = dat, seg.length = 100,
                            series = flagged2)
## Select series by name or column position
ccf.100.2 <- ccf.series.rwl(rwl = co021, seg.length = 100,
                            series = "641143")
ccf.100.3 <- ccf.series.rwl(rwl = co021, seg.length = 100,
                            series = which(colnames(co021) == "641143"))
identical(ccf.100.1, ccf.100.2) # TRUE
identical(ccf.100.2, ccf.100.3) # TRUE

## End(Not run)

Check and Validate an rwl Object

Description

Checks that an object is a valid rwl and, if not, attempts coercion via as.rwl. Warns about internal NA values found within any series.

Usage

check.rwl(rwl, why = NULL, bai.ok = FALSE)

Arguments

Details

check.rwl is for functions that want ring widths. If rwl is an "rwi" object of ring-width indices (see as.rwi), a warning names the calling function, says that it wants widths, adds why, and suggests as.rwl if the values really are widths. The object is then relabelled as class "rwl" and checked as below. A "bai" object of basal area increment (see as.bai) is treated the same way, unless bai.ok is TRUE, when it is relabelled without a warning.

If rwl is not already of class "rwl", coercion is attempted using as.rwl. If coercion succeeds, a warning is issued. If coercion fails (e.g., because row names are not consecutive integers, columns are not numeric, or the input is not a data.frame or matrix), the function stops with a single informative message that describes both the class problem and the reason coercion failed.

After class validation, each series is checked for internal NA values, defined as NA values that are sandwiched between non-NA values within a series (as opposed to leading or trailing NAs, which are standard in rwl objects). If any are found, a warning names the affected series and suggests fill.internal.NA.

This function is called at the top of the dplR functions that want ring widths, providing uniform validation across the package. Functions that take widths or indices alike, and those that want indices, use internal checks of their own (see as.rwi).

Value

An object of class c("rwl", "data.frame") with series in columns and years as row names. The series IDs are the column names. The returned object is identical to the input when the input is already a valid rwl.

Author(s)

Andy Bunn

See Also

as.rwl, as.rwi, read.rwl, fill.internal.NA

Examples

library(utils)
data(ca533)

## A proper rwl passes silently and is returned unchanged
ca533.checked <- check.rwl(ca533)
identical(ca533, ca533.checked)

## A plain data.frame is coerced with a warning.
## suppressWarnings() is used here to keep R CMD check clean since the
## warning is expected. Remove it to see the coercion warning.
ca533.df <- ca533
class(ca533.df) <- "data.frame"
ca533.rwl <- suppressWarnings(check.rwl(ca533.df))
#ca533.rwl <- check.rwl(ca533.df)
class(ca533.rwl)

Build Mean Value Chronology

Description

This function builds a mean value chronology, typically from a data.frame of detrended ring widths as produced by detrend.

Usage

chron(rwi, biweight = TRUE, prewhiten = FALSE, ...)

Arguments

Details

This either averages the rows of the data.frame using a mean or a robust mean (the so-called standard chronology) or can do so from the residuals of an AR process (the residual chronology).

Note that the residual chronology in this function will return different values than the residual chronology from chron.ars which uses a slightly different method for determining AR order.

Value

An object of of class crn and data.frame with the standard chronology, residual chronology (if prewhitening was performed), and the sample depth. The years are stored as row numbers.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

References

Cook, E. R. and Kairiukstis, L. A., editors (1990) Methods of Dendrochronology: Applications in the Environmental Sciences. Springer. ISBN-13: 978-0-7923-0586-6.

Fritts, H. C. (2001) Tree Rings and Climate. Blackburn. ISBN-13: 978-1-930665-39-2.

See Also

read.rwl, detrend, ar, crn.plot

Examples

library(graphics)
library(utils)
data(ca533)
ca533.rwi <- detrend(rwl = ca533, method = "ModNegExp")
ca533.crn <- chron(ca533.rwi)
plot(ca533.crn,xlab="Year",ylab="RWI")
## With residual chron
ca533.crn2 <- chron(ca533.rwi, prewhiten = TRUE)
plot(ca533.crn2,xlab="Year",ylab="RWI")

Build ARSTAN Chronology

Description

This function builds three varieties of the mean-value chronology, including the ARSTAN chronology, typically from a data.frame of detrended ring widths as produced by detrend.

Usage

chron.ars(rwi, biweight=TRUE, maxLag=10, firstAICmin=TRUE,
  verbose=TRUE, prewhitenMethod=c("ar.yw","arima.CSS-ML"))

Arguments

Details

This produces three mean-value chronologies: standard, residual, and ARSTAN. Users unfamiliar with the concept behind the ARSTAN method should look to Cook (1985) for background and inspiration.

The standard chronology is the (biweight) mean value across rows and identical to chron.

The residual chronology is the prewhitened chronology as described by Cook (1985) and uses multivariate autoregressive modeling to determine the order of the AR process. It is important to note that the residual chronology produced here is different from the simple residual chronology produced by chron, which returns the residuals of an AR process using a naive call to ar. In practice the results will be similar. For more on the residual chronology in this function, see pp. 153-154 in Cook's 1985 dissertation.

The ARSTAN chronology builds on the residual chronology but returns a re-whitened chronology where the pooled AR coefficients from the multivariate autoregressive modeling are reintroduced. See references for details.

The order of the AR model is selected from the pooled AR coefficients by AIC using either the first local AIC minimum (firstAICmin = TRUE, the default) or the overall minimum across all lags up to maxLag (firstAICmin = FALSE). When firstAICmin = TRUE and the AIC decreases monotonically across all tested lags without reaching a local minimum, the function stops with an error. This can occur when maxLag is too small for the data or when the series have persistence structures that are not well-suited to this approach (e.g., non-stationary or very long-memory series). In such cases, consider increasing maxLag or reviewing the data preparation steps described in Cook (1985). When firstAICmin = FALSE the function always selects an order and does not stop in this way.

If the selected AR order is zero (no detectable common autocorrelation among the series), no prewhitening or post-reddening is applied and all three returned chronologies are identical to the standard chronology.

Once the AR order is determined an AR(p) model is fit to each series using either ar via the Yule-Walker method or by arima via conditional-sum-of-squares to find starting values, then maximum likelihood. It is possible that the model will not converge in which case a warning is produced. The AR fitting is determined via prewhitenMethod and defaults to using ar.

Value

A data.frame of class "crn" with the following columns:

Row names are the years taken from rownames(rwi).

Author(s)

Andy Bunn with contributions from Kevin Achukaitis and Ed Cook. Much of the function is a port of Cook's FORTRAN code.

References

Cook, E. R. and Kairiukstis, L. A., editors (1990) Methods of Dendrochronology: Applications in the Environmental Sciences. Springer. ISBN-13: 978-0-7923-0586-6.

Cook, E. R. (1985). A Time Series Analysis Approach to Tree Ring Standardization. PhD thesis, The University of Arizona.

See Also

chron, crn.plot, ar, arima

Examples

library(graphics)
library(utils)
data(co021)
co021.rwi <- detrend(rwl = co021, method = "AgeDepSpline")
co021.crn <- chron.ars(co021.rwi)
plot(co021.crn,xlab="Year",ylab="RWI",add.spline=TRUE,nyrs=20)
cor(co021.crn)

Build Mean Value Chronology with Confidence Intervals

Description

This function builds a mean value chronology with bootstrapped confidence intervals, typically from a data.frame of detrended ring widths as produced by detrend.

Usage

chron.ci(x, biweight=TRUE, conf=0.95, R=100)

Arguments

Details

This either averages the rows of the data.frame using a mean or a robust mean (the so-called standard chronology) and calculates boostrapped confidence intervals using the normal approximation. The function will fail if there are any rows in x that contain only one sample and in practice there should be several samples in a row. One of the guiding principles of bootstrapping is that the population is to the sample as the sample is to the bootstrap samples.

Value

An object of of class data.frame with the standard chronology, the upper and lower confidence interval, and the sample depth. The years are stored as row numbers.

Author(s)

Andy Bunn.

See Also

read.rwl, detrend, boot, boot.ci

Examples

library(graphics)
library(utils)
data(wa082)
# truncate to a sample depth of five
wa082Trunc <- wa082[rowSums(!is.na(wa082))>4,]
# detrend
wa082RWI <- detrend(wa082Trunc, method="AgeDepSpline")
# bootstrap the chronology and
wa082Crn <- chron.ci(wa082RWI, biweight = TRUE, R = 100, conf = 0.99)
head(wa082Crn)
# plot (this is so much easier in ggplot!)
wa082Crn$yrs <- time(wa082Crn)
xx <- c(wa082Crn$yrs,rev(wa082Crn$yrs))
yy <- c(wa082Crn$lowerCI,rev(wa082Crn$upperCI))
plot(wa082Crn$yrs,wa082Crn$std,type="n",ylim=range(yy),
     ylab="RWI",xlab="Year",main="Chronology with CI")
polygon(x=xx,y=yy,col = "grey",border = NA)
lines(wa082Crn$yrs,wa082Crn$std)

Build Mean Value Chronology with Stabilized Variance

Description

This function builds a variance stabilized mean-value chronology, typically from a data.frame of detrended ring widths as produced by detrend.

Usage

chron.stabilized(x, winLength, biweight = TRUE, running.rbar = FALSE)

Arguments

Details

The variance of a mean chronology depends on the variance of the individual samples, the number of series averaged together, and their interseries correlation (Wigley et al. 1984). As the number of series commonly decreases towards the beginning of a chronology averaging introduces changes in variance that are a solely an effect of changes in sample depth.

Additionally, time-dependent changes in interseries correlation can cause artificial variance changes of the final mean chronology. The function chron.stabilized accounts for both temporal changes in the interseries correlation and sample depth to produce a mean value chronology with stabilized variance.

The basic correction centers around the use of the effective independent sample size, Neff, which considers sample replication and mean interseries correlation between the samples at every time. This is defined as: Neff = n(t) / 1+(n(t)-1)rbar(t)

where n(t) is the number of series at time t, and rbar is the interseries correlation (see interseries.cor). Multiplication of the mean time series with the square root of Neff at every time t theoretically results in variance that is independent of sample size. In the limiting cases, when the rbar is zero or unity, Neff obtains values of the true sample size and unity, respectively.

Value

An object of of class crn and data.frame with the variance stabilized chronology, running interseries correlation ('if running.rbar=TRUE), and the sample depth.

Author(s)

Original code by David Frank and adapted for dplR by Stefan Klesse. Patched and improved by Andy Bunn.

References

Frank, D, Esper, J, Cook, E, (2006) On variance adjustments in tree-ring chronology development. Tree rings in archaeology, climatology and ecology, TRACE 4, 56–66

Frank, D, Esper, J, Cook, E, (2007) Adjustment for proxy number and coherence in a large-scale temperature reconstruction. Geophysical Research Letters 34

Wigley, T, Briffa K, Jones P (1984) On the Average Value of Correlated Time Series, with Applications in Dendroclimatology and Hydrometeorology. J. Climate Appl. Meteor., 23, 201–213

See Also

chron

Examples

library(graphics)
library(utils)
data(co021)
co021.rwi <- detrend(co021,method = "Spline")
co021.crn <- chron(co021.rwi)
co021.crn2 <- chron.stabilized(co021.rwi,
                                winLength=101,
                                biweight = TRUE,
                                running.rbar = FALSE)
yrs <- time(co021)
plot(yrs,co021.crn$std,type="l",col="grey")
lines(yrs,co021.crn2$adj.crn,col="red")

C-Method Standardization

Description

Detrend multiple ring-width series simultaneously using the C-method.

Usage

cms(rwl, po, c.hat.t = FALSE, c.hat.i = FALSE)

Arguments

Details

This method detrends and standardizes tree-ring series by calculating a growth curve based on constant annual basal area increment. The method is based on the “assumption that constant growth is expressed by a constant basal area increment distributed over a growing surface” (Biondi and Qeadan 2008). The detrending is the estimation and removal of the tree’s natural biological growth trend. The standardization is done by dividing each series by the growth trend to produce units in the dimensionless ring-width index (RWI).

This attempts to remove the low frequency variability that is due to biological or stand effects.

A series with no values at all is dropped, with a message naming it, and the rest are standardized as if it had never been there. po must still have a row for every series in rwl, empty or not.

See the reference below for further details.

Value

An object of class c("rwi", "data.frame") (see as.rwi) containing the dimensionless and detrended ring-width indices with column names, row names and dimensions of rwl if c.hat.t is FALSE and c.hat.i is FALSE.

Otherwise a list of length 2 or 3 containing the RWI data.frame, a data.frame containing the C-curves for each tree (c.hat.t), and/or a vector containing the C-values for each tree (c.hat.i) depending on the output flags. See Eq. 12 in Biondi and Qeadan (2008) for more detail on c.hat.t, and c.hat.i.

Note

DendroLab website: https://dendrolaborg.wordpress.com/

Author(s)

Code provided by DendroLab based on programming by F. Qeadan and F. Biondi, University of Nevada Reno, USA and adapted for dplR by Andy Bunn. Patched and improved by Mikko Korpela.

References

Biondi, F. and Qeadan, F. (2008) A theory-driven approach to tree-ring standardization: Defining the biological trend from expected basal area increment. Tree-Ring Research, 64(2), 81–96.

See Also

detrend, chron, rcs

Examples

library(graphics)
library(utils)
data(gp.rwl)
data(gp.po)
gp.rwi <- cms(rwl = gp.rwl, po = gp.po)
gp.crn <- chron(gp.rwi)
crn.plot(gp.crn, add.spline = TRUE)
## c.hat
gp.rwi <- cms(rwl = gp.rwl, po = gp.po, c.hat.t = TRUE, c.hat.i = TRUE)
dotchart(gp.rwi$c.hat.i, ylab = "Series", xlab = expression(hat(c)[i]))
tmp <- gp.rwi$c.hat.t
plot(tmp[, 1], type = "n", ylim = range(tmp, na.rm = TRUE),
     xlab = "Cambial Age", ylab = expression(hat(c)[t]))
apply(tmp, 2, lines)

Schulman Old Tree No. 1, Mesa Verde

Description

This data set gives the raw ring widths for Douglas fir Pseudotsuga menziesii at Mesa Verde in Colorado, USA. There are 35 series. Data set was created using read.rwl and saved to an .rda file using save.

Usage

data(co021)

Format

A data.frame containing 35 tree-ring series in columns and 788 years in rows.

Source

International tree-ring data bank, Accessed on 20-April-2021 at https://www.ncei.noaa.gov/pub/data/paleo/treering/measurements/northamerica/usa/co021.rwl

References

Schulman, E. (1963) Schulman Old Tree No. 1 Data Set. IGBP PAGES/World Data Center for Paleoclimatology Data Contribution Series 1983-CO021.RWL. NOAA/NCDC Paleoclimatology Program, Boulder, Colorado, USA.

Examples

library(utils)
data(co021)

## Where the data came from. read.tucson() records what it saw and attaches
## it to what it returns; this data set carries that record.
prov <- attr(co021, "dplR.provenance")
prov$file          # the archived file it was read from
prov$precision     # the precision each series was measured at
prov$gaps          # interior gaps, if any -- this file has none
prov$events        # problems found while reading -- none here either

## rwl.report() puts the file and precision at the head of its report
rwl.report(co021)

Combine Tree-Ring Data Sets

Description

This function combines any number of data.frames of tree-ring data into one data.frame.

Usage

combine.rwl(x, y = NULL)

Arguments

Details

The sequence of years in each data.frame must be increasing and continuous. The output produced by the function also fulfills this condition. If the input is differently formatted, the result will be wrong.

Value

An object of class c("rwl", "data.frame") with the series in columns and the years as rows. The keycodes are the column names and the years are the row names.

Author(s)

Christian Zang. Patched by Mikko Korpela.

Examples

library(utils)
data(ca533)
data(co021)
combi1 <- combine.rwl(list(ca533, co021))
## or alternatively for data.frames to combine
combi2 <- combine.rwl(ca533, co021)
identical(combi1, combi2) # TRUE

Common Interval

Description

This function finds the common interval on a set of tree-ring widths such as that produced by read.rwl.

Usage

common.interval(rwl, type=c("series", "years", "both"),
                make.plot=TRUE)

Arguments

Details

This trims an rwl object to a common interval that maximizes the number of series (type="series"), the number of years (type="years"), or a compromise between the two (type="both"). A modified seg.plot can be drawn as well.

Series with no values are left out. A single series is its own common interval, and its measured years are returned. If there are two or more series and no year in which any two of them overlap, or no series has any values, there is no common interval and that is an error.

Value

A data.frame with colnames(x) and rownames(x), with no missing values. It has the class of rwl: ring-width indices (class "rwi", see as.rwi) come back as indices. If the years kept are not consecutive, which type = "years" can give when a series has an interior gap, the result is a plain data.frame, with a warning (see [.rwl).

Author(s)

Filipe Campelo, Andy Bunn and Mikko Korpela

See Also

seg.plot

Examples

library(utils)
data(co021)
co021.s <- common.interval(co021, type="series", make.plot=TRUE)
co021.y <- common.interval(co021, type="years", make.plot=TRUE)
co021.b <- common.interval(co021, type="both", make.plot=TRUE)

dim(co021)
dim.s <- dim(co021.s)
dim.s       # the highest number of series
prod(dim.s) #   (33 series x 288 years = 9504)
dim.y <- dim(co021.y)
dim.y       # the highest number of years
prod(dim.y) #   (27 series x 458 years = 12366)
dim.b <- dim(co021.b)
dim.b       # compromise solution
prod(dim.b) #   (28 series x 435 years = 12180)

Compute Correlations between Series

Description

Computes the correlation between each tree-ring series in a rwl object.

Usage

corr.rwl.seg(rwl, seg.length = 50, bin.floor = 100, n = NULL,
             nyrs = NULL, prewhiten = TRUE, ar.order.max = NULL,
             pcrit = 0.05, biweight = TRUE,
             method = c("spearman", "pearson","kendall"),
             make.plot = TRUE, label.cex = 1, floor.plus1 = FALSE,
             master = NULL, lag.max = 0,
             master.yrs = as.numeric(if (is.null(dim(master))) {
                              names(master)
                          } else {
                              rownames(master)
                          }),
             ...)

Arguments

Details

This function calculates correlation serially between each tree-ring series and a master chronology built from all the other series in the rwl object (leave-one-out principle). Optionally, the user may give a master chronology (a vector) as an argument. In the latter case, the same master chronology is used for all the series in the rwl object. The user can also choose to give a master data.frame (series as columns, years as rows), from which a single master chronology is built.

Correlations are done for each segment of the series where segments are lagged by half the segment length (e.g., 100-year segments would be overlapped by 50-years). The first segment is placed according to bin.floor. The minimum bin year is calculated as ceiling(min.yr/bin.floor)*bin.floor where min.yr is the first year in either the rwl object or the user-specified master chronology, whichever is smaller. For example if the first year is 626 and bin.floor is 100 then the first bin would start in 700. If bin.floor is 10 then the first bin would start in 630.

Correlations are calculated for the first segment, then the second segment and so on. Correlations are only calculated for segments with complete overlap with the master chronology. For now, correlations are Spearman’s rho calculated via cor.test using method = "spearman".

Each series (including those in the rwl object) is optionally detrended as the residuals from a hanning filter with weight n. The filter is not applied if n is NULL. Detrending can also be done via prewhitening where the residuals of an ar model are added to each series mean. This is the default. The master chronology is computed as the mean of the rwl object using tbrm if biweight is TRUE and rowMeans if not. Note that detrending can change the length of the series. E.g., a hanning filter will shorten the series on either end by floor(n/2). The prewhitening default will change the series length based on the ar model fit. The effects of detrending can be seen with series.rwl.plot.

As an alternative to the hanning filter, nyrs divides each series by a smoothing spline, which removes no years from the ends of the series. COFECHA uses a 32-year spline. The spline alone does not usually save years at the start of a series, because the ar model is then fitted to a high-pass filtered series, and AIC tends to choose a high order for it. Setting ar.order.max to a small value (e.g., 3) limits the years lost to prewhitening, and after a spline a low-order model loses little of the fit.

With neither filter, each series is divided by its mean before the master is built, so every series counts equally. Data that can be negative, such as indices from detrend with difference = TRUE, log widths or isotope values, have the mean subtracted instead, since dividing by a negative mean would flip the series. The n and nyrs filters divide by a smooth curve, so they refuse such data. The same applies to the other crossdating functions.

The function is typically invoked to produce a plot where each segment for each series is colored by its correlation to the master chronology. Green segments are those that do not overlap completely with the width of the bin. Blue segments are those that correlate above the user-specified critical value. Red segments are those that correlate below the user-specified critical value and might indicate a dating problem.

A segment that correlates well where it is dated can still correlate better somewhere else. With lag.max greater than 0, each segment is also correlated against the master at every shift from -lag.max to lag.max years. Each year of the segment is paired with the master value k years away, so the correlation at lag k is that of series[t] with master[t + k]. The sign follows ccf.series.rwl with series.x = FALSE: a negative lag means missing rings in the series, a positive lag means false (extra) rings. The dated position wins ties.

This gives the two kinds of flag COFECHA reports. An ‘A’ segment correlates below the critical value, but its dated position is still the best one tested: it is weak, not misdated. A ‘B’ segment correlates better at some other position, whether or not it is significant as dated. The letters are not stored; they follow from the returned matrices as best.lag != 0 for B and best.lag == 0 & p.val >= pcrit for A (see Examples). When lag.max is greater than 0, the plot draws B segments in purple and leaves red for the rest of the segments below the critical value.

A B flag is a hypothesis about the dating, not a verdict. A one-year shift that gains 0.05 in correlation is noise as often as it is a missing ring, which is why the size of the gain is returned (best.rho - spearman.rho) alongside the lag, so that it can be thresholded. Check any B segment against the wood, for example with ccf.series.rwl and xskel.ccf.plot.

Read the lags along a series rather than one bin at a time. Dating runs from the bark inward, so a missing ring moves every ring before it one year late: every complete bin before the missing ring has best.lag of -1 and every bin after it has 0. A false ring does the same with +1. The error therefore lies where the lag changes, not in the bin with the largest gain. The bin that straddles it holds some shifted and some unshifted years, and may show either lag with a small gain. A run of -1 with correctly dated bins on both sides is a different fault: the ring count is right, but a ring is missing at the later end of the run and there is an extra ring at the earlier end. A locally absent ring entered as a zero in too early a year looks exactly like this. A test on the whole series, such as RWL_DATING_LAG in rwl.check, usually cannot see it, because most of the series is dated correctly. The Examples show both patterns.

The shifted correlations follow the same complete-overlap rule as the dated one: a lag is tested only if the master has a value for every year of the shifted window. Near the start and end of the record (and of the master), the lags that run off the record cannot be tested, and a segment that lies within lag.max years of an end is only searched in one direction. A segment that is not complete at its dated position is not tested at all. Floating series and series with dating errors near their ends sit exactly where the test is blind, so the absence of a flag there says little. Worse, a segment misdated in the direction that runs off the record can still be flagged B, at whichever of the remaining lags happens to correlate best, so the lag reported for a B near either end can be noise even when the segment really is misdated.

Value

A list containing matrices spearman.rho, p.val, overall, bins, rwi, vector avg.seg.rho, numeric seg.lag, seg.length, pcrit, label.cex, matrices best.lag, best.rho and numeric lag.max. An additional character flags is also returned if any segments fall below the critical value. Matrix spearman.rho contains the correlations for each series by bin. Matrix p.val contains the p-values on the correlation for each series by bin. Matrix overall contains the average correlation and p-value for each series. Matrix bins contains the years encapsulated by each bin. The vector avg.seg.rho contains the average correlation for each bin. Matrix rwi contains the detrended rwl data, the numerics seg.lag, seg.length, pcrit, label.cex are from the oroginal call and used to pass into plot.crs.

Matrices best.lag and best.rho have the same shape and names as spearman.rho. best.lag gives, for each series and bin, the lag (in years) at which the segment correlates best with the master, and best.rho gives that correlation. Both are NA where spearman.rho is. With lag.max = 0, best.lag is 0 and best.rho equals spearman.rho. lag.max is from the original call. Flags in flags depend on the p-value only and are the same whatever lag.max is.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

corr.series.seg, skel.plot, series.rwl.plot, ccf.series.rwl, plot.crs

Examples

library(utils)
data(co021)
crs <- corr.rwl.seg(co021, seg.length = 100, label.cex = 1.25)
names(crs)
## Average correlation and p-value for the first few series
head(crs$overall)
## Average correlation for each bin
crs$avg.seg.rho

## Plant a missing ring and follow it through the bins: delete the
## 1500 ring of series 641143, the same fault as in the examples for
## ccf.series.rwl(), corr.series.seg() and xskel.ccf.plot().  Dated
## from the bark, every ring before 1500 now sits one year late.
dat <- co021
x <- dat$"641143"
names(x) <- rownames(dat)
dat$"641143" <- delete.ring(x, year = 1500)
crs <- corr.rwl.seg(dat, lag.max = 5)
## -1 in every complete bin before 1500 (drawn in purple), 0 from 1500
## on, and a large gain in each shifted bin
ok <- !is.na(crs$best.lag["641143", ])
rbind(lag = crs$best.lag["641143", ok],
      rho = round(crs$spearman.rho["641143", ok], 2),
      best.rho = round(crs$best.rho["641143", ok], 2))

## The same missing ring with a false ring inserted at 1400.  The ring
## count is right again, and only the rings between the two errors are
## one year late.
dat$"641143" <- insert.ring(delete.ring(x, year = 1500), year = 1400)
crs <- corr.rwl.seg(dat, lag.max = 5, make.plot = FALSE)
## -1 between the two errors only, 0 on both sides
crs$best.lag["641143", ok]

## COFECHA-style A and B flags for every series
flag <- ifelse(crs$best.lag != 0, "B",
               ifelse(crs$p.val >= crs$pcrit, "A", ""))
flag[is.na(flag)] <- ""
## B segments with their lag and the gain in correlation
idx <- which(flag == "B", arr.ind = TRUE)
data.frame(series = rownames(flag)[idx[, 1]],
           bin = colnames(flag)[idx[, 2]],
           lag = crs$best.lag[idx],
           rho = round(crs$spearman.rho[idx], 3),
           best.rho = round(crs$best.rho[idx], 3))

Compute Correlation between a Series and a Master Chronology

Description

Compute correlation between a tree-ring series and a master chronology by segment.

Usage

corr.series.seg(rwl, series, series.yrs = as.numeric(names(series)),
                seg.length = 50, bin.floor = 100, n = NULL,
                nyrs = NULL, prewhiten = TRUE, ar.order.max = NULL,
                biweight = TRUE,
                method = c("spearman", "pearson","kendall"),
                pcrit = 0.05,
                make.plot = TRUE, floor.plus1 = FALSE, ...)

Arguments

Details

This function calculates the correlation between a tree-ring series and a master chronology built from a rwl object. Correlations are done by segment (see below) and with a moving correlation with length equal to the seg.length. The function is typically invoked to produce a plot.

Value

A list containing matrices bins, moving.rho, and vectors spearman.rho, p.val, and overall.

Matrix bins contains the years encapsulated by each bin (segments). Matrix moving.rho contains the moving correlation and p-value for a moving average equal to seg.length. Vector spearman.rho contains the correlations by bin and p.val contains the p-values. Vector overall contains the average correlation and p-value.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

corr.series.seg, skel.plot, series.rwl.plot, ccf.series.rwl

Examples

library(utils)
data(co021)
dat <- co021
## Create a missing ring: delete the 1500 ring of series 641143 and
## drop the original from the master.  Dated from the bark, every ring
## before 1500 now sits one year late.
flagged <- dat$"641143"
names(flagged) <- rownames(dat)
flagged <- delete.ring(flagged, year = 1500)
dat$"641143" <- NULL
seg.100 <- corr.series.seg(rwl = dat, series = flagged,
                           seg.length = 100, biweight = FALSE)
## Not run: 
flagged2 <- co021$"641143"
names(flagged2) <- rownames(dat)
seg.100.1 <- corr.series.seg(rwl=dat, seg.length=100, biweight=FALSE,
                             series = flagged2)
## Select series by name or column position
seg.100.2 <- corr.series.seg(rwl=co021, seg.length=100, biweight=FALSE,
                             series = "641143")
seg.100.3 <- corr.series.seg(rwl=co021, seg.length=100, biweight=FALSE,
                             series = which(colnames(co021) == "641143"))
identical(seg.100.1, seg.100.2) # TRUE
identical(seg.100.2, seg.100.3) # TRUE

## End(Not run)

Detrend Multiple Ring-Width Series Simultaneously

Description

This is a wrapper for detrend.series to detrend many ring-width series at once.

Usage

detrend(rwl, y.name = names(rwl), make.plot = FALSE,
        method = "Spline",
        nyrs = NULL, f = 0.5, pos.slope = FALSE,
        constrain.nls = c("never", "when.fail", "always"),
        verbose = FALSE, return.info = FALSE,
        wt, span = "cv", bass = 0, difference = FALSE)

Arguments

Details

See detrend.series for details on detrending methods. Setting make.plot = TRUE will cause plots of each series to be produced. These could be saved using Devices if desired.

A series with no values at all is dropped, with a message naming it, and the rest are detrended as if it had never been there; its name is dropped from y.name too. x[rows, ] and head can leave such a series (see [.rwl). If no series has any values, that is an error. A series with NA inside its measured span stops the call, with a message naming the series; see fill.internal.NA.

Value

If one detrending method is used, an object of class c("rwi", "data.frame") containing the dimensionless detrended ring widths with the column names, row names and dimensions of rwl, less any series with no values. It records how it was made in attr(x, "dplR.detrend") and carries the provenance record of rwl, if it has one; see as.rwi. If more methods are used, a list with ncol(rwl) elements each containing a data.frame with the detrended ring widths in each column.

If return.info is TRUE, the return value is a list with four parts:

Note

This function uses the foreach looping construct with the %dopar% operator. For parallel computing and a potential speedup, a parallel backend must be registered before running the function. If verbose is TRUE, parallel computation is disabled.

Author(s)

Andy Bunn. Improved by Mikko Korpela.

See Also

detrend.series, as.rwi

Examples

library(utils)
data(ca533)
## Detrend using modified exponential decay. Returns an rwi object
ca533.rwi <- detrend(rwl = ca533, method = "ModNegExp")
## Detrend using splines and compute
## residuals via subtraction
ca533.rwi <- detrend(rwl = ca533, method = "Spline",
                     difference = TRUE)

## Detrend using modified Hugershoff curve and return info on the model
## fits. Returns a list with: series, curves, modelinfo and data.info
data(co021)
co021.rwi <- detrend(rwl = co021, method = "ModHugershoff",
                     return.info=TRUE)

## Not run: 
library(grDevices)
## Detrend using all methods. Returns a list
ca533.rwi <- detrend(rwl = ca533,
                     method = c("Spline", "ModNegExp", "Mean", "Ar",
                                "Friedman", "ModHugershoff", "AgeDepSpline"))
## Save a pdf of all series
fname <- tempfile(fileext=".pdf")
print(fname) # tempfile used for output
pdf(fname)
ca533.rwi <- detrend(rwl = ca533, method = c("Spline", "ModNegExp"),
                     make.plot = TRUE)
dev.off()

unlink(fname) # remove the file

## End(Not run)

Detrend a Ring-Width Series

Description

Detrend a tree-ring series by one of two methods, a smoothing spline or a statistical model. The series and fits are plotted by default.

Usage

detrend.series(y, y.name = "", make.plot = TRUE,
               method = "Spline",
               nyrs = NULL, f = 0.5, pos.slope = FALSE,
               constrain.nls = c("never", "when.fail", "always"),
               verbose = FALSE, return.info = FALSE,
               wt, span = "cv", bass = 0, difference = FALSE)

Arguments

Details

This detrends and standardizes a tree-ring series. The detrending is the estimation and removal of the tree’s natural biological growth trend. The default standardization is done by dividing each series by the growth trend to produce units in the dimensionless ring-width index (RWI). If difference is TRUE, the index is calculated by subtraction. Values of zero (typically missing rings) in y are set to 0.001 when dividing.

A curve that is divided by must be positive, so when difference is FALSE a fit that is not all positive is replaced by a straight line or the mean (or, for "Ar", its negative values are set to zero) with a warning. Subtraction has no such need, so with difference = TRUE the fits are used as they are, zeros are left alone, and y may be negative, as log widths or isotope values are.

There are currently seven methods available for detrending although more are certainly possible. The methods implemented are an age-dependent spline via ads (method = "AgeDepSpline"), the residuals of an AR model (method = "Ar"), Friedman's Super Smoother (method = "Friedman"), a simple horizontal line (method = "Mean"), or a modified Hugershoff curve (method = "ModHugershoff"), a modified negative exponential curve (method = "ModNegExp"), and a smoothing spline via caps (method = "Spline").

The "AgeDepSpline" approach uses an age-dependent spline with an initial stiffness of 50 (nyrs=50). See ads. If some of the fitted values are not positive then method "Mean" is used. However, this is extremely unlikely.

The "Ar" approach is also known as prewhitening where the detrended series is the residuals of an ar model divided by the mean of those residuals to yield a series with white noise and a mean of one. This method removes all but the high frequency variation in the series and should only be used as such.

The "Friedman" approach uses Friedman’s ‘super smoother’ as implemented in supsmu. The parameters wt, span and bass can be adjusted, but periodic is always set to FALSE. If some of the fitted values are not positive then method "Mean" is used.

The "Mean" approach fits a horizontal line using the mean of the series. This method is the fallback solution in cases where the "Spline" or the linear fit (also a fallback solution itself) contains zeros or negative values, which would lead to invalid ring-width indices.

The "ModHugershoff" approach attempts to fit a Hugershoff model of biological growth of the form f(t) = a t^b e^{-g t} + d, where the argument of the function is time, using nls. See Fritts (2001) for details about the parameters. Option constrain.nls gives a possibility to constrain the parameters of the modified negative exponential function. If the constraints are enabled, the nonlinear optimization algorithm is instructed to keep the parameters in the following ranges: a \ge 0, b \ge 0 and d \ge 0. The default is to not constrain the parameters (constrain.nls = "never") for nls but warn the user when the parameters go out of range.

If a suitable nonlinear model cannot be fit (function is non-decreasing or some values are not positive) then a linear model is fit. That linear model can have a positive slope unless pos.slope is FALSE in which case method "Mean" is used.

The "ModNegExp" approach attempts to fit a classic nonlinear model of biological growth of the form f(t) = a e^{b t} + k, where the argument of the function is time, using nls. See Fritts (2001) for details about the parameters. Option constrain.nls gives a possibility to constrain the parameters of the modified negative exponential function. If the constraints are enabled, the nonlinear optimization algorithm is instructed to keep the parameters in the following ranges: a \ge 0, b \le 0 and k \ge 0. The default is to not constrain the parameters (constrain.nls = "never") for nls but warn the user when the parameters go out of range.

If a suitable nonlinear model cannot be fit (function is non-decreasing or some values are not positive) then a linear model is fit. That linear model can have a positive slope unless pos.slope is FALSE in which case method "Mean" is used.

The "Spline" approach uses a spline where the frequency response is 0.50 at a wavelength of 0.67 * “series length in years”, unless specified differently using nyrs and f in the function caps. If some of the fitted values are not positive then method "Mean" is used.

These methods are chosen because they are commonly used in dendrochronology. There is a rich literature on detrending and many researchers are particularly skeptical of the use of the classic nonlinear model of biological growth (f(t) = a e^{b t} + k) for detrending. It is, of course, up to the user to determine the best detrending method for their data.

Note that the user receives a warning if a series has negative values in the fitted curve. This happens fairly commonly with with the ‘Ar’ method on high-order data. It happens less often with method ‘Spline’ but isn't unheard of (see ‘Examples’). If this happens, users should look carefully at their data before continuing. Automating detrending and not evaluating each series individually is folly. Remember, frustration over detrending is the number one cause of dendros going to live as hermits in the tallgrass prairie, where there are no trees to worry about.

See the references below for further details on detrending. It's a dark art.

Value

If several methods are used, returns a data.frame containing the detrended series (y) according to the methods used. The columns are named and ordered to match method. If only one method is selected, returns a vector.

If return.info is TRUE, the return value is a list with four parts:

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela. A bug fix related to negative output values is based on work by Alice Cecile.

References

Cook, E. R. and Kairiukstis, L. A., editors (1990) Methods of Dendrochronology: Applications in the Environmental Sciences. Springer. ISBN-13: 978-0-7923-0586-6.

Fritts, H. C. (2001) Tree Rings and Climate. Blackburn. ISBN-13: 978-1-930665-39-2.

See Also

detrend

Examples

library(stats)
library(utils)
## Use series CAM011 from the Campito data set
data(ca533)
series <- ca533[, "CAM011"]
names(series) <- rownames(ca533)
# all seven methods
series.rwi <- detrend.series(y = series, y.name = "CAM011", verbose=TRUE,
                             method = c("Spline", "ModNegExp", "Mean", "Ar",
                                        "Friedman", "ModHugershoff",
                                        "AgeDepSpline"))
# see plot with three methods
series.rwi <- detrend.series(y = series, y.name = "CAM011",
                             method=c("Spline", "ModNegExp","Friedman"),
                             difference=TRUE)
# see plot with two methods
# interesting to note difference from ~200 to 250 years 
# in terms of what happens to low frequency growth
series.rwi <- detrend.series(y = series, y.name = "CAM011",
                             method=c("Spline", "ModNegExp"))
# see plot with just one method and change the spline
# stiffness to 50 years (which is not neccesarily a good choice!)
series.rwi <- detrend.series(y = series, y.name = "CAM011",
                             method="Spline",nyrs=50)
                             
# note that method "Ar" doesn't get plotted in first panel
# since this approach doesn't approximate a growth curve.
series.rwi <- detrend.series(y = series, y.name = "CAM011",
                             method="Ar")
                             
# note the difference between ModNegExp and ModHugershoff at the 
# start of the series. Also notice how curves, etc. are returned
# via return.info
data(co021)
series <- co021[, 4]
names(series) <- rownames(co021)
series.rwi <- detrend.series(y = series, y.name = names(co021)[4],
                             method=c("ModNegExp", "ModHugershoff"),
                             verbose = TRUE, return.info = TRUE, 
                             make.plot = TRUE)
                             
# A dirty dog.
# In the case of method=="Spline" the function carries-on
# and applies method=="Mean" as an alternative. 
data(nm046)
series <- nm046[,8]
names(series) <- rownames(nm046)
series.rwi <- detrend.series(y = series, y.name = names(nm046)[8],
                             method="Spline",
                             make.plot = FALSE)

Defunct functions in dplR

Description

These are functions no longer included in dplR

Details

Fill Internal NA

Description

This function fills internal NA values (i.e., those with numeric data above and below a small data gap) in each column of a data.frame such as a data set of ring widths as produced by read.rwl.

Usage

fill.internal.NA(x, fill = c("Mean", "Spline", "Linear"))

Arguments

Details

There are occasionally data gaps within a tree-ring series. Some of the functions in dplR will fail when an internal NA is encountered (e.g. caps). This function fills internal NA values with either a given numeric value (e.g., 0) or through crude imputation. The latter can be calculated as the mean of the series (fill="Mean") or calculated by fitting a cubic spline (fill="Spline") using the spline function or calculated by linear approximation (fill="Linear") using the function approx.

Editorial: Having internal NA in a tree-ring series is often bad practice and filling those values should be done with caution. For instance, some users code missing rings as NA instead of 0. And missing values (i.e., NA) are sometimes present in maximum latewood density data when the rings are small. A common, but not recommended, practice is to leave stretches of NA values in places where it has been impossible to accurately measure rings (perhaps because of a break in the core). It is often better to treat that core as two separate series (e.g., "01A" and "01B" rather than have internal NA values. As with all processing, the analyst should make a decision based on their experience with the wood and not rely on software to make a choice for them!

Value

A data.frame with colnames(x) and rownames(x). Internal NAs filled as above.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

spline, approx

Examples

library(graphics)
library(stats)
foo <- data.frame(x1=c(rnorm(5), NA, NA, rnorm(3)),
                  x2=c(rnorm(10)),
                  x3=c(NA, NA, rnorm(3), NA, rnorm(4)),
                  x4=c(NA, NA, rnorm(3), NA, rnorm(3), NA),
                  x5=c(NA, NA, rnorm(8)),
                  x6=c(NA, rnorm(9)),
                  x7=c(NA, rnorm(5), NA, rnorm(3)),
                  x8=c(rnorm(8), NA, NA),
                  x9=c(rnorm(5), NA, rnorm(3), NA))
row.names(foo) <- 1901:1910
class(foo) <- c("rwl","data.frame")
fill.internal.NA(foo, fill=0)

bar <- fill.internal.NA(foo, fill="Spline")
baz <- fill.internal.NA(foo, fill="Linear")

## note differences in method "Spline" vs. "Linear"
yrs <- time(foo)
plot(yrs, foo$x7, type="b", lwd=3)
lines(yrs, bar$x7, col="red", lwd=2)
lines(yrs, baz$x7, col="green", lwd=1)

Calculate the Gini Coefficient

Description

This function calculates the Gini coefficient on raw or detrended ring-width series.

Usage

gini.coef(x)

Arguments

Details

This calculates the Gini coefficient of inequality which is used as an all-lag measure of diversity in tree-ring records – typically detrended series. Lower values indicate lower diversity. The use of the Gini coefficient in dendrochronology is described by Biondi and Qeadan (2008). See Handcock and Morris (1999) for more information.

Value

the Gini coefficient.

Note

The equivalence between the Lorenz-curve formula used in the C code and the relative-mean-difference definition of the Gini coefficient is derived in the vignette: vignette("math-dplR", package = "dplR").

Author(s)

Mikko Korpela, based on original by Andy Bunn

References

Biondi, F. and Qeadan, F. (2008) Inequality in Paleorecords. Ecology, 89(4), 1056–1067.

Handcock, M. S. and Morris, M. (1999) Relative Distribution Methods in the Social Sciences. Springer. ISBN: 0-387-98778-9.

See Also

rwl.stats

Examples

library(utils)
data(ca533)
ca533.rwi <- detrend(rwl = ca533, method = "ModNegExp")
ca533.crn <- chron(ca533.rwi)
gini.coef(ca533.crn)

Calculate Gleichläufigkeit

Description

This function calculates the Gleichläufigkeit and related measures for a given set of tree-ring records.

Usage

  glk(x, overlap = 50, prob = TRUE)
  glk.legacy(x)

Arguments

Details

Gleichläufigkeit is a classical agreement test based on sign tests (Eckstein and Bauch, 1969). This function implements Gleichläufigkeit as the pairwise comparison of all records in data set. This vectorized implementation is faster than the previous version and follows the original definition (Huber 1942), instead of the incorrect interpretation that has been used in the past (Schweingruber 1988, see Buras/Wilmking 2015 for the correction).

The probability of exceedence (p) for the Gleichläufigkeit expresses the chance that the Gleichläufigkeit is incorrect. The observed value of the Gleichläufigkeit is converted to a z-score and based on the standard normal curve the probability of exceedence is calculated. The result is a matrix of all p-values (Jansma 1995, 60-61, see also Visser 2020).

Note that prior to dplR version 1.7.2, glk did not have the overlap or prob and returned a matrix with just the Gleichläufigkeit for all possible combinations of records. That function can still be accessed via glk.legacy.

Value

The funtions returns a named list of two or three matrices (p_mat is optional if prob = TRUE):

1. glk_mat: matrix with Gleichläufigkeit

2. overlap: matrix with number of overlapping growth changes.This is the number of overlapping years minus one.

3. p_mat: matrix of all probabilities of exceedence for all observed Gleichläufigkeit values.

The matrices can be extracted from the list by selecting the name or index number. If two curves have less than 3 years of overlap, Gleichläufigkeit cannot be computed, and NA is returned. To calculate the global glk of the dataset mean(x$glk_mat, na.rm = TRUE). See examples for calculating the the global glk without self-similarities (i.e., the ones on the diagonal).

Author(s)

Christian Zang. Patched and improved by Mikko Korpela. Improved by Allan Buras. Further improved and expanded by Ronald Visser and Andy Bunn

References

Buras, A. and Wilmking, M. (2015) Correcting the calculation of Gleichläufigkeit, Dendrochronologia 34, 29-30. DOI: https://doi.org/10.1016/j.dendro.2015.03.003

Eckstein, D. and Bauch, J. (1969) Beitrag zur Rationalisierung eines dendrochronologischen Verfahrens und zur Analyse seiner Aussagesicherheit. Forstwissenschaftliches Centralblatt, 88(1), 230-250.

Huber, B. (1943) Über die Sicherheit jahrringchronologischer Datierung. Holz als Roh- und Werkstoff 6, 263-268. DOI: https://doi.org/10.1007/BF02603303

Jansma, E., 1995. RemembeRINGs; The development and application of local and regional tree-ring chronologies of oak for the purposes of archaeological and historical research in the Netherlands, Nederlandse Archeologische Rapporten 19, Rijksdienst voor het Oudheidkundig Bodemonderzoek, Amersfoort

Schweingruber, F. H. (1988) Tree rings: basics and applications of dendrochronology, Kluwer Academic Publishers, Dordrecht, Netherlands, 276 p.

Visser, R.M. (2020) On the similarity of tree-ring patterns: Assessing the influence of semi-synchronous growth changes on the Gleichläufigkeit for big tree-ring data sets,Archaeometry, 63, 204-215 DOI: https://doi.org/10.1111/arcm.12600

See Also

sgc sgc is an alternative for glk)

Examples

library(utils)
  data(ca533)
  ca533.glklist <- glk(ca533)
  ca533.glk_mat <- ca533.glklist$glk_mat
  # calculating the mean GLK including self-similarities
  mean(ca533.glk_mat, na.rm = TRUE) 
  # calculating the mean GLK excluding self-similarities
  mean(ca533.glk_mat[upper.tri(ca533.glk_mat)], na.rm = TRUE) 

Ponderosa Pine Distance to Pith Corresponding to gp.rwl

Description

This data set gives the distance to pith for each series (in mm) that matches the ring widths for gp.rwl – a data set of ponderosa pine (Pinus ponderosa) from the Gus Pearson Natural Area (GPNA) in northern Arizona, USA. Data are further described by Biondi and Qeadan (2008) and references therein.

Usage

data(gp.d2pith)

Format

A data.frame containing series IDs in column 1 (series) and the distance (in mm) from the innermost ring to the pith of the tree (d2pith). This can be used together with the ring widths to calculate the area of each ring.

Source

DendroLab, University of Nevada Reno, USA. https://dendrolaborg.wordpress.com/

References

Biondi, F. and Qeadan, F. (2008) A theory-driven approach to tree-ring standardization: Defining the biological trend from expected basal area increment. Tree-Ring Research, 64(2), 81–96.

Ponderosa Pine Stem Diameters and Bark Thickness (gp.rwl)

Description

This data set gives the diameter at breast height for each series that matches the series in gp.rwl – a data set of ponderosa pine (Pinus ponderosa) from the Gus Pearson Natural Area (GPNA) in northern Arizona, USA. Data are further described by Biondi and Qeadan (2008) and references therein.

Usage

data(gp.dbh)

Format

A data.frame containing series IDs in column 1 (series), tree diameter (in cm) at breast height (dbh), and the bark thickness (in cm). This can be used together with the ring widths to calculate the area of each ring.

Source

DendroLab, University of Nevada Reno, USA. https://dendrolaborg.wordpress.com/

References

Biondi, F. and Qeadan, F. (2008) A theory-driven approach to tree-ring standardization: Defining the biological trend from expected basal area increment. Tree-Ring Research, 64(2), 81–96.

Ponderosa Pine Pith Offsets Corresponding to gp.rwl

Description

This data set gives the pith offsets that match the ring widths for gp.rwl – a data set of ponderosa pine (Pinus ponderosa) from the Gus Pearson Natural Area (GPNA) in northern Arizona, USA. Data are further described by Biondi and Qeadan (2008) and references therein.

Usage

data(gp.po)

Format

A data.frame containing series IDs in column 1 (series) and the number of years between the beginning of that series in gp.rwl and the pith of the tree (pith.offset). This can be used together with the ring widths to calculate the cambial age of each ring.

Source

DendroLab, University of Nevada Reno, USA. https://dendrolaborg.wordpress.com/

References

Biondi, F. and Qeadan, F. (2008) A theory-driven approach to tree-ring standardization: Defining the biological trend from expected basal area increment. Tree-Ring Research, 64(2), 81–96.

Ponderosa Pine Ring Widths from Gus Pearson Natural Area

Description

This data set includes ring-width measurements for ponderosa pine (Pinus ponderosa) increment cores collected at the Gus Pearson Natural Area (GPNA) in northern Arizona, USA. There are 58 series from 29 trees (2 cores per tree). Data are further described by Biondi and Qeadan (2008) and references therein.

Usage

data(gp.rwl)

Format

A data.frame containing 58 ring-width series in columns and 421 years in rows.

Source

DendroLab, University of Nevada Reno, USA. https://dendrolaborg.wordpress.com/

References

Biondi, F. and Qeadan, F. (2008) A theory-driven approach to tree-ring standardization: Defining the biological trend from expected basal area increment. Tree-Ring Research, 64(2), 81–96.

Hanning Filter

Description

Applies a Hanning filter of length n to x.

Usage

hanning(x, n = 7)

Arguments

Details

This applies a low frequency Hanning (a.k.a. Hann) filter to x with weight set to n.

Value

A filtered vector.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

References

Oppenheim, A. V., Schafer, R. W., and Buck, J. R. (1999) Discrete-Time Signal Processing. Prentice-Hall, 2nd edition. ISBN-13: 978-0-13-754920-7.

See Also

filter

Examples

library(graphics)
library(utils)
data(ca533)
yrs <- time(ca533)
y <- ca533[, 1]
not.na <- !is.na(y)
yrs <- yrs[not.na]
y <- y[not.na]
plot(yrs, y, xlab = "Years", ylab = "Series1 (mm)",
     type = "l", col = "grey")
lines(yrs, hanning(y, n = 9), col = "red", lwd = 2)
lines(yrs, hanning(y, n = 21), col = "blue", lwd = 2)
legend("topright", c("Series", "n=9", "n=21"),
       fill=c("grey", "red", "blue"))

Interactively Detrend Multiple Ring-Width Series

Description

Interactively detrend multiple tree-ring series by one of two methods, a smoothing spline or a statistical model. This is a wrapper for detrend.series.

Usage

i.detrend(rwl, y.name = names(rwl), nyrs = NULL, f = 0.5,
          pos.slope = FALSE)

Arguments

Details

This function allows a user to choose detrending curves based on plots that are produced by detrend.series for which it is essentially a wrapper. The user enters their choice of detrended method via keyboard at a prompt for each ring width series in rwl. See detrend.series for examples and details on the detrending methods.

A series with no values at all is dropped, with a message naming it, before any series is shown, and its name is dropped from y.name too.

Value

An object of class c("rwi", "data.frame") (see as.rwi) containing each detrended series according to the method used as columns and rownames set to colnames(y). These are typically years. The method chosen for each series is recorded, named by series, in attr(x, "dplR.detrend")$method. Plots are also produced as the user chooses the detrending methods through keyboard input.

Author(s)

Andy Bunn

See Also

detrend.series

Interactively Detrend a Ring-Width Series

Description

Interactively detrend a tree-ring series by one of three methods, a smoothing spline, a linear model, or the mean. This is a wrapper for detrend.series.

Usage

i.detrend.series(y, y.name = NULL, nyrs = NULL, f = 0.5,
                 pos.slope = FALSE)

Arguments

Details

This function allows a user to choose a detrending method based on a plot that is produced by detrend.series for which it is essentially a wrapper. The user enters their choice of detrended method via keyboard at a prompt. See detrend.series for examples and details on the detrending methods.

Value

A vector containing the detrended series (y) according to the method used with names set to colnames(y). These are typically years. The name of the method chosen is attached as attr(x, "method"). A plot is also produced and the user chooses a method through keyboard input.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

detrend.series

Edit a Ring-Width Series

Description

Insert or delete rings from a ring-width series

Usage

insert.ring(rw.vec,rw.vec.yrs=as.numeric(names(rw.vec)),
            year,ring.value=mean(rw.vec,na.rm=TRUE),
            fix.last=TRUE,fix.length=TRUE)
delete.ring(rw.vec,rw.vec.yrs=as.numeric(names(rw.vec)),
            year,fix.last=TRUE,fix.length=TRUE)

Arguments

Details

Simple editing of ring widths.

Value

A named vector.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

dplR

Examples

library(utils)
data(gp.rwl)
dat <- gp.rwl
# insert a value of zero for the year 1950 in series 50A
# fix the last year of growth and maintain the length of the series
tmp <- insert.ring(rw.vec=dat$"50A",rw.vec.yrs = time(dat),
                   year=1950,ring.value=0,fix.length = TRUE)
# with fix.length=TRUE this can be merged back into the rwl object:
data.frame(dat$"50A",tmp)
dat$"50A" <- tmp

# note that if fix.last = FALSE and fix.length = FALSE inserting a ring causes the
# ending year of the series to be pushed forward and the length of the output to
# be longer than the original series.
tmp <- insert.ring(rw.vec=dat$"50A",rw.vec.yrs = time(dat),
                   year=1950,ring.value=0, fix.last = FALSE, 
                   fix.length = FALSE)
# with fix.length=FALSE this can't be merged back into the rwl object the
# same way as above
tail(tmp)
length(tmp)
nrow(dat)

# the same logic applies to deleting rings.
dat <- gp.rwl
# delete the year 1950 in series 50A
# fix the last year of growth and maintain the length of the series
tmp <- delete.ring(rw.vec=dat$"50A",rw.vec.yrs = time(dat),
                   year=1950,fix.last = TRUE, fix.length = TRUE)
# with fix.length=TRUE this can be merged back into the rwl object:
data.frame(dat$"50A",tmp)
dat$"50A" <- tmp

# note that if fix.last = FALSE and fix.length = FALSE inserting a ring causes the
# ending year of the series to be pushed forward and the length of the output to
# be longer than the original series.
tmp <- delete.ring(rw.vec=dat$"50A", rw.vec.yrs = time(dat),
                   year=1950, fix.last = FALSE, 
                   fix.length = FALSE)
# with fix.length=FALSE this can't be merged back into the rwl object the
# same way as above
tail(tmp)
length(tmp)
nrow(dat)

Individual Series Correlation Against a Master Chronology

Description

This function calculates the correlation between a series and a master chronology.

Usage

  interseries.cor(rwl, n = NULL, nyrs = NULL, prewhiten = TRUE,
    ar.order.max = NULL, biweight = TRUE,
    method = c("spearman", "pearson", "kendall"))

Arguments

Details

This function calculates correlation serially between each tree-ring series and a master chronology built from all the other series in the rwl object (leave-one-out principle).

Each series in the rwl object is optionally detrended as the residuals from a hanning filter with weight n. The filter is not applied if n is NULL. Detrending can also be done via prewhitening where the residuals of an ar model are added to each series mean. This is the default. The master chronology is computed as the mean of the rwl object using tbrm if biweight is TRUE and rowMeans if not. Note that detrending can change the length of the series. E.g., a hanning filter will shorten the series on either end by floor(n/2). The prewhitening default will change the series length based on the ar model fit. The effects of detrending can be seen with series.rwl.plot.

This function produces the same output of the overall portion of corr.rwl.seg. The mean correlation value given is sometimes referred to as the “overall interseries correlation” or the “COFECHA interseries correlation”. This output differs from the rbar statistics given by rwi.stats in that rbar is the average pairwise correlation between series where this is the correlation between a series and a master chronology.

Value

a data.frame with correlation values and p-values given from cor.test

Author(s)

Andy Bunn, patched and improved by Mikko Korpela

See Also

rwl.stats, rwi.stats

Examples

library(utils)
data(gp.rwl)
foo <- interseries.cor(gp.rwl)
# compare to: 
# corr.rwl.seg(rwl=gp.rwl,make.plot=FALSE)$overall
# using pearson's r
foo <- interseries.cor(gp.rwl,method="pearson")

# two measures of interseries correlation
# compare interseries.cor to rbar from rwi.stats
gp.ids <- read.ids(gp.rwl, stc = c(0, 2, 1))
bar <- rwi.stats(gp.rwl, gp.ids, prewhiten=TRUE)
bar$rbar.eff
mean(foo[,1])

Date Conversion to Character in LaTeX Format

Description

This is a simple convenience function that returns a date in the format used by ‘⁠\today⁠’ in LaTeX. A possible use case is fixing the date shown in a vignette at weaving time.

Usage

latexDate(x = Sys.Date(), ...)

Arguments

Value

A character vector

Author(s)

Mikko Korpela

Examples

latexDate()                              # today
latexDate(Sys.Date() + 5)                # today + 5 days
latexDate(c("2013-12-06", "2014-09-19")) # fixed dates
## [1] "December 6, 2013"   "September 19, 2014"
latexDate(5*60*60*24, origin=Sys.Date()) # today + 5 days

Convert Character Strings for Use with LaTeX

Description

Some characters cannot be entered directly into a LaTeX document. This function converts the input character vector to a form suitable for inclusion in a LaTeX document in text mode. It can be used together with ‘⁠\Sexpr⁠’ in vignettes.

Usage

latexify(x, doublebackslash = TRUE, dashdash = TRUE,
         quotes = c("straight", "curved"),
         packages = c("fontenc", "textcomp"))

Arguments

Details

The function is intended for use with unformatted inline text. Newlines, tabs and other whitespace characters ("[:space:]" in regex) are converted to spaces. Control characters ("[:cntrl:]") that are not whitespace are removed. Other more or less special characters in the ASCII set are ‘{’, ‘}’, ‘\’, ‘#’, ‘$’, ‘%’, ‘^’, ‘&’, ‘_’, ‘~’, double quote, ‘/’, single quote, ‘<’, ‘>’, ‘|’, grave and ‘-’. They are converted to the corresponding LaTeX commands. Some of the conversions are affected by user options, e.g. dashdash.

Before applying the substitutions described above, input elements with Encoding set to "bytes" are printed and the output is stored using captureOutput. The result of this intermediate stage is ASCII text where some characters are shown as their byte codes using a hexadecimal pair prefixed with "\x". This set includes tabs, newlines and control characters. The substitutions are then applied to the intermediate result.

The quoting functions sQuote and dQuote may use non-ASCII quote characters, depending on the locale. Also these quotes are converted to LaTeX commands. This means that the quoting functions are safe to use with any LaTeX input encoding. Similarly, some other non-ASCII characters, e.g. letters, currency symbols, punctuation marks and diacritics, are converted to commands.

Adding "eurosym" to packages enables the use of the euro sign as provided by the "eurosym" package (‘⁠\euro⁠’).

The result is converted to UTF-8 encoding, Normalization Form C (NFC). Note that this function will not add any non-ASCII characters that were not already present in the input. On the contrary, some non-ASCII characters, e.g. all characters in the "latin1" (ISO-8859-1) Encoding (character set), are removed when converted to LaTeX commands. Any remaining non-ASCII character has a good chance of working when the document is processed with XeTeX or LuaTeX, but the Unicode support available with pdfTeX is limited.

Assuming that ‘⁠pdflatex⁠’ is used for compilation, suggested package loading commands in the document preamble are:

\usepackage[T1]{fontenc}    % no '"' in OT1 font encoding
\usepackage{textcomp}       % some symbols e.g. straight single quote
\usepackage[utf8]{inputenx} % UTF-8 input encoding
\input{ix-utf8enc.dfu}      % more supported characters

Value

A character vector

Author(s)

Mikko Korpela

References

INRIA. Tralics: a LaTeX to XML translator, HTML documentation of all TeX commands. http://www-sop.inria.fr/marelle/tralics/.

Levitt, N., Persch, C., and Unicode, Inc. (2013) GNOME Character Map, application version 3.10.1.

Mittelbach, F., Goossens, M., Braams, J., Carlisle, D., and Rowley, C. (2004) The LaTeX Companion. Addison-Wesley, second edition. ISBN-13: 978-0-201-36299-2.

Pakin, S. (2009) The Comprehensive LaTeX Symbol List. https://www.ctan.org/tex-archive/info/symbols/comprehensive.

The Unicode Consortium. The Unicode Standard. https://home.unicode.org/.

Examples

x1 <- "clich\xe9\nma\xf1ana"
Encoding(x1) <- "latin1"
x1
x2 <- x1
Encoding(x2) <- "bytes"
x2
x3 <- enc2utf8(x1)
testStrings <-
    c("different     kinds\nof\tspace",
      "control\a characters \ftoo",
      "{braces} and \\backslash",
      '#various$ %other^ &characters_ ~escaped"/coded',
      x1,
      x2,
      x3)
latexStrings <- latexify(testStrings, doublebackslash = FALSE)
## All should be "unknown"
Encoding(latexStrings)
cat(latexStrings, sep="\n")
## Input encoding does not matter
identical(latexStrings[5], latexStrings[7])

Perform a Continuous Morlet Wavelet Transform

Description

This function performs a continuous wavelet transform on a time series.

Usage

morlet(y1, x1 = seq_along(y1), p2 = NULL, dj = 0.25, siglvl = 0.95)

Arguments

Details

This performs a continuous wavelet transform of a time series. This function is typically invoked with wavelet.plot.

Value

A list containing:

Note

This is a port of Torrence’s IDL code, which can be accessed through the Internet Archive Wayback Machine.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

References

Torrence, C. and Compo, G. P. (1998) A practical guide to wavelet analysis. Bulletin of the American Meteorological Society, 79(1), 61–78.

See Also

wavelet.plot

Examples

library(utils)
data(ca533)
ca533.rwi <- detrend(rwl = ca533, method = "ModNegExp")
ca533.crn <- chron(ca533.rwi, prewhiten = FALSE)
Years <- time(ca533.crn)
CAMstd <- ca533.crn[, 1]
out.wave <- morlet(y1 = CAMstd, x1 = Years, dj = 0.1, siglvl = 0.99)

Calculate NET

Description

Computes the \mathit{NET} parameter for a set of tree-ring records or other time-series data.

Usage

net(x, weights = c(v = 1, g = 1))

Arguments

Details

This function computes the \mathit{NET} parameter (Esper et al., 2001). The overall \mathit{NET} is an average of all (non-NA) yearly values \mathit{NET_j}, which are computed as follows:

\mathit{NET_j}=v_j+(1-G_j)

The yearly variation v_j is the standard deviation of the measurements of a single year divided by their mean. Gegenläufigkeit 1-G_j is based on one definition of Gleichläufigkeit G_j, similar to but not the same as what glk computes. Particularly, in the formula used by this function (Esper et al., 2001), simultaneous zero differences in two series are not counted as a synchronous change.

The weights of v_j and 1-G_j in the sum can be adjusted with the argument weights (see above). As a rather extreme example, it is possible to isolate variation or Gegenläufigkeit by setting one of the weights to zero (see ‘Examples’).

Value

A list with the following components, in the same order as described here:

Author(s)

Mikko Korpela

References

Esper, J., Neuwirth, B., and Treydte, K. (2001) A new parameter to evaluate temporal signal strength of tree-ring chronologies. Dendrochronologia, 19(1), 93–102.

Examples

library(utils)
data(ca533)
ca533.rwi <- detrend(rwl = ca533, method = "ModNegExp")
ca533.net <- net(ca533.rwi)
tail(ca533.net$all)
ca533.net$average
## Not run: 
## Isolate the components of NET
ca533.v <- net(ca533.rwi, weights=c(v=1,0))
ca533.g <- net(ca533.rwi, weights=c(g=1,0))

## End(Not run)

Los Alamos Tree Ring Widths

Description

This data set gives the raw ring widths for Douglas fir Pseudotsuga menziesii near Los Alamos New Mexico, USA. There are 8 series. Data set was created using read.rwl and saved to an .rda file using save.

Usage

data(nm046)

Format

A data.frame containing 8 tree-ring series in columns and 289 years in rows.

Source

International tree-ring data bank, Accessed on 20-April-2021 at https://www.ncei.noaa.gov/pub/data/paleo/treering/measurements/northamerica/usa/nm046.rwl

References

O'Brien, D (2002) Los Alamos Data Set. IGBP PAGES/World Data Center for Paleoclimatology Data Contribution Series 2002-NM046.RWL. NOAA/NCDC Paleoclimatology Program, Boulder, Colorado, USA.

Examples

library(utils)
data(nm046)

## Where the data came from. read.tucson() records what it saw and attaches
## it to what it returns; this data set carries that record.
prov <- attr(nm046, "dplR.provenance")
prov$file          # the archived file it was read from
prov$precision     # the precision each series was measured at
prov$gaps          # interior gaps, if any -- this file has none
prov$events        # problems found while reading -- none here either

## rwl.report() puts the file and precision at the head of its report
rwl.report(nm046)

Low-pass, high-pass, band-pass, and stop-pass filtering

Description

Applies low-pass, high-pass, band-pass, or stop-pass filtering to y with frequencies (or periods) supplied by the user.

Usage

pass.filt(y, W, type = c("low", "high", "stop", "pass"),
          method = c("Butterworth", "ChebyshevI"),
          n = 4, Rp = 1)

Arguments

Details

This function applies either a Butterworth or a Chebyshev type I filter of order n to a signal and is nothing more than a wrapper for functions in the signal package. The filters are designed via butter and cheby1. The filter is applied via filtfilt.

The input data (y) has the mean value subtracted and is then padded via reflection at the start and the end to a distance of twice the maximum period. The padded data and the filter are passed to filtfilt after which the data are unpadded and returned afer the mean is added back.

The argumement W can be given in either frequency between 0 and 0.5 or, for convenience, period (minimum value of 2). For low-pass and high-pass filters, W must have a length of one. For low-pass and high-pass filters W must be a two-element vector (c(low, high)) specifying the lower and upper boundaries of the filter.

Because this is just a wrapper for casual use with tree-ring data the frequencies and periods assume a sampling frequency of one. Users are encouraged to build their own filters using the signal package.

Value

A filtered vector.

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

hanning, detrend

Examples

data("co021")
x <- na.omit(co021[, 1])

# 20-year low-pass filter -- note freq is passed in
bSm <- pass.filt(x, W=0.05, type="low", method="Butterworth")
cSm <- pass.filt(x, W=0.05, type="low", method="ChebyshevI")
plot(x, type="l", col="grey")
lines(bSm, col="red")
lines(cSm, col="blue")

# 20-year high-pass filter -- note period is passed in
bSm <- pass.filt(x, W=20, type="high")
plot(x, type="l", col="grey")
lines(bSm, col="red")

# 20 to 100-year band-pass filter -- note freqs are passed in
bSm <- pass.filt(x, W=c(0.01, 0.05), type="pass")
cSm <- pass.filt(x, W=c(0.01, 0.05), type="pass", method="ChebyshevI")
plot(x, type="l", col="grey")
lines(bSm, col="red")
lines(cSm, col="blue")

# 20 to 100-year stop-pass filter -- note periods are passed in
cSm <- pass.filt(x, W=c(20, 100), type="stop", method="ChebyshevI")
plot(x, type="l", col="grey")
lines(cSm, col="red")

Plot a Tree-Ring Chronology

Description

This function makes a default plot of a tree-ring chronology from a data.frame of the type produced by chron, chron.ars, chron.stabilized, ssf.

Usage

## S3 method for class 'crn'
plot(x, add.spline = FALSE, nyrs = NULL, ...)

Arguments

Details

This makes a crude plot of one or more tree-ring chronologies.

Value

None. Invoked for side effect (plot).

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

chron

Examples

library(graphics)
data(wa082)
# Truncate the RW data to a sample depth at least 5
wa082Trunc <- wa082[rowSums(!is.na(wa082))>4,]
# Detrend with age-dependent spline
wa082RWI <- detrend(wa082Trunc,method = "AgeDep")
# make several chronologies
wa082CRN1 <- chron(wa082RWI)
wa082CRN2 <- chron.stabilized(wa082RWI,
                              winLength=51,
                              biweight = TRUE,
                              running.rbar = TRUE)
wa082CRN3 <- chron.ars(wa082RWI)
wa082CRN4 <- ssf(wa082Trunc)

# and plot
plot.crn(wa082CRN1,add.spline = TRUE,nyrs=20)
plot.crn(wa082CRN2,add.spline = TRUE,nyrs=20)
plot(wa082CRN3,add.spline = TRUE,nyrs=20)
plot(wa082CRN4,add.spline = TRUE,nyrs=20)

# a custom crn
foo <- data.frame(wa082CRN1,sfc=wa082CRN4$sfc)
foo <- foo[,c(1,3,2)]
class(foo) <- c("crn","data.frame")
plot.crn(foo,add.spline = TRUE,nyrs=20)

Plotting crs Objects

Description

Plots objects returned from corr.rwl.seg.

Usage

## S3 method for class 'crs'
plot(x, ...)

Arguments

Details

Each segment is drawn in blue if it correlates above the critical value, red if it does not, and green where it does not overlap the bin completely. If corr.rwl.seg was run with lag.max greater than 0, segments that correlate better at another lag (COFECHA’s B flag) are drawn in purple instead, whatever their correlation as dated.

Value

None. A plot is produced.

Author(s)

Andy Bunn

See Also

corr.rwl.seg

Examples

library(graphics)
data(co021)
foo <- corr.rwl.seg(co021, make.plot = FALSE)
plot(foo)

Plotting Rwl Objects

Description

Plots rwl objects.

Usage

## S3 method for class 'rwl'
plot(x, plot.type=c("seg","spag"), ...)

Arguments

Value

None. A plot is produced.

Author(s)

Andy Bunn

See Also

read.rwl

Examples

library(graphics)
library(utils)
data(co021)
plot(co021, plot.type="seg")
plot(co021, plot.type="spag")
plot(co021, plot.type="spag", zfac=2)

Convert Pith Offset to Wood Completeness

Description

This function creates a partial wood completeness data structure based on pith offset data.

Usage

po.to.wc(po)

Arguments

Details

Uses pith.offset - 1 as the number of missing heartwood rings.

Value

A data.frame containing one variable of wood completeness data: n.missing.heartwood (integer type). This can be used as input to write.tridas.

Author(s)

Mikko Korpela

See Also

wc.to.po, rcs, write.tridas

Examples

## Not run: 
library(utils)
data(gp.po)
all(wc.to.po(po.to.wc(gp.po)) == gp.po)

## End(Not run)

Calculates Pointer Years from a Group of Ring-Width Series

Description

This function calculates pointer years on a data.frame of ring-width series using the Becker algorithm. The pointer years are computed with adjustable thresholds of relative radial growth variation and number of series displaying similar growth pattern (i.e. positive or negative variations).

Usage

pointer(rwl, rgv.thresh = 10, nseries.thresh = 75, round.decimals = 2)

Arguments

Details

This calculates pointer years from ring-width series for each year t of the time period covered by the series using the Becker algorithm. This algorithm is based on, first, the calculation of the individual relative radial growth variation by comparison of ring-width of year t to that of year t-1 for each series, and second, the inter-series comparison of both sign and magnitude of these variations.

For example, if rgv.thresh and nseries.thresh are set at 10 and 75 respectively, pointer years will be defined as those years when at least 75% of the series present an absolute relative radial growth variation higher than 10%.

Users unfamiliar with the Becker algorithm should refer to Becker et al. (1994) and Mérian and Lebourgeois (2011) for further details.

Value

A data.frame containing the following columns (each row corresponds to one position of the window):

Author(s)

Pierre Mérian. Improved by Mikko Korpela and Andy Bunn.

References

Becker, M., Nieminen, T. M., and Gérémia, F. (1994) Short-term variations and long-term changes in oak productivity in northeastern France – the role of climate and atmospheric CO2. Annals of Forest Science, 51(5), 477–492.

Mérian, P. and Lebourgeois, F. (2011) Size-mediated climate–growth relationships in temperate forests: A multi-species analysis. Forest Ecology and Management, 261(8), 1382–1391.

See Also

skel.plot

Examples

## Pointer years calculation on ring-width series. Returns a data.frame.
library(utils)
data(gp.rwl)
py <- pointer(rwl=gp.rwl, rgv.thresh=10, nseries.thresh=75,
              round.decimals=2)
tail(py)

Power Transformation of Tree-Ring Data

Description

Power transformation of tree-ring width.

Usage

powt(rwl, method = "universal", rescale = FALSE,
     return.power=FALSE)

Arguments

Details

In dendrochronology, ring width series are sometimes power transformed to address heteroscedasticity.

The classic procedure used by method="cook" is a variance stabilization technique implemented after Cook & Peters (1997): for each series a linear model is fitted on the logs of level and spread, where level is defined as the local mean M_t = \left(R_t + R_{t-1}\right)/2 with ring widths R, and spread S is the local standard deviation defined as S_t = \left|R_t - R_{t-1}\right|. The regression coefficient b from a linear model \log S = k + b \log M is then used for the power transform \star{R}_t = R_t^{1-b}.

The procedure above is modified with method="universal" where all samples are used simultaneously in a linear mixed-effects model with time (year) as a random effect: lmer(log S ~ log M + (1|year). This "universal" or "signal free" approach accounts for the common year effect across all of the series in rwl and should address that not every year has the same change in environmental conditions to the previous year.

The rescale argument will return the series with a mean and standard deviation that matches the input data. While this is a common convention, users should note that this can produce negative values which can be confusing if thought of as "ring widths."

Value

Either an object of class c("rwl", "data.frame") containing the power transformed ring width series with the series in columns and the years as rows or in the case of a single series, a possibly named vector of the same. With class rwl, the series IDs are the column names and the years are the row names.

If return.power=TRUE the returned object is a list containing the power transformed data and a numeric with the power estimate(s) used to transform the data.

Author(s)

Christian Zang implemented the Cook and Peters method. Stefan Klesse conceived and wrote the universal method. Patched and improved by Mikko Korpela and Andy Bunn.

References

Cook, E. R. and Peters, K. (1997) Calculating unbiased tree-ring indices for the study of climatic and environmental change. The Holocene, 7(3), 361–370.

Examples

library(utils)
data(zof.rwl)
powtUniversal <- powt(zof.rwl, method = "universal")
powtCook <- powt(zof.rwl, method = "cook")

op <- par(no.readonly = TRUE)
par(mfcol = c(1, 3))
hist(summary(zof.rwl)$skew,
     breaks = seq(-2.25,2.25,by=0.25),
     main="Raw Data",xlab="Skew")
hist(summary(powtUniversal)$skew,
     breaks = seq(-2.25,2.25,by=0.25),
     main="Universal POWT",xlab="Skew")
hist(summary(powtCook)$skew,
     breaks = seq(-2.25,2.25,by=0.25),
     main="Cook POWT",xlab="Skew")
par(op) # restore graphical parameters

Printing Redfit Results

Description

Print information contained in or derived from a redfit object.

Usage

## S3 method for class 'redfit'
print(x, digits = NULL, csv.out = FALSE, do.table = FALSE,
      prefix = "", row.names = FALSE, file = "", ...)

Arguments

Value

Invisibly returns x.

Author(s)

Mikko Korpela

References

This function is based on the Fortran program REDFIT, which is in the public domain.

Schulz, M. and Mudelsee, M. (2002) REDFIT: estimating red-noise spectra directly from unevenly spaced paleoclimatic time series. Computers & Geosciences, 28(3), 421–426.

See Also

redfit

Examples

library(utils)
data(ca533)
tm <- time(ca533)
x <- ca533[[1]]
idx <- which(!is.na(x))
redf <- redfit(x[idx], tm[idx], "time",
               nsim = 100, iwin = 0, ofac = 1, n50 = 1)
print(redf)
fname <- tempfile(fileext=".csv")
print(fname) # tempfile used for output
print(redf, csv.out = TRUE, file = fname)
redftable <- read.csv(fname)
unlink(fname) # remove the file

Do some reporting on a RWL object

Description

This function prints the results of rwl.report

Usage

## S3 method for class 'rwl.report'
print(x, ...)

Arguments

Details

This function formats the list from rwl.report for the user to have a summary report of the number of series, the mean length of all the series, the first year, last year, the mean first-order autocorrelation (via summary.rwl), the mean interseries correlation (via interseries.cor), the years where a series has a missing ring (zero), internal NA, or a very small ring (<0.005).

Value

Invisible

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

See Also

rwl.report, summary.rwl, interseries.cor

Examples

data("gp.rwl")
rwl.report(gp.rwl)
foo <- gp.rwl
foo[177,1] <- NA 
foo[177:180,3] <- NA 
foo[185,4] <- 0.001 
rwl.report(foo)

Add Raster Elements to Plot

Description

This function takes plotting commands and uses a temporary bitmap graphics device to capture their output. The resulting raster image is drawn in the plot or figure region of the active high-level plot. A new plot is started if one does not exist.

Usage

rasterPlot(expr, res = 150, region = c("plot", "figure"), antialias,
           bg = "transparent", interpolate = TRUE, draw = TRUE,
           Cairo = FALSE, ...)

Arguments

Details

The appropriate graphical parameters of the current graphics device are copied to the temporary bitmap device. Therefore the appearance of the raster contents should be almost the same as when directly drawn.

The call or expression expr is evaluated in the environment of the caller.

It is possible that the raster contents will maintain a constant size when the graphics device is resized. If resizing works, however, the image may become distorted. For example, circle symbols will turn into ellipses if the width to height ratio is not maintained (see ‘Examples’). This is in contrast to a standard plot in a display graphics device, e.g. x11, where text and symbols maintain their size when the device is resized.

Value

If draw is TRUE, there is no return value. The function is used for the side effects.

If draw is FALSE, an object of class "nativeRaster" is returned. The object can be used as input for rasterImage or grid.raster. See readPNG. If no bitmap device is available (see ‘Note’), NULL is returned.

Note

• The graphics device used for the output must have support for including raster images. See "rasterImage" in dev.capabilities.

• The R build must have a functional png device, which requires one of the following capabilities: "png", "aqua" or "cairo". Alternatively, a Cairo device from package Cairo must be available with Cairo.capabilities "raster" or "png".

If either of these requirements is not met, at least one message is generated and the function reverts to regular plotting. The bg argument is then handled by drawing a filled rectangle. Also region is honored, but the other settings do not apply.

Author(s)

Mikko Korpela

Examples

library(graphics)
library(stats)

## Picture with various graphical elements
x <- 1:100
y0 <- quote(sin(pi * x / 20) + x / 100 + rnorm(100, 0, 0.2))
y <- eval(y0)
ylab <- deparse(y0)
spl <- smooth.spline(y)
plot(x, y, type = "n", axes = FALSE, ylab = ylab)
usr <- par("usr")
xrange <- usr[2] - usr[1]
xsize <- xrange * 0.4
nsteps <- 8
xmar <- xsize / 20
yrange <- usr[4] - usr[3]
ysize <- yrange / 20
ymar <- 0.5 * ysize
X <- seq(usr[1] + xmar, by = xsize / nsteps, length.out = nsteps + 1)
xleft <- X[-(nsteps + 1)]
xright <- X[-1]
pin <- par("pin")
maxrad <- xsize / 3 * min(1, pin[2] / pin[1])
nrad <- 16
minrad <- maxrad / nrad
Rad <- seq(maxrad, by = (minrad - maxrad) / (nrad - 1), length.out=nrad)
xmar2 <- xmar + maxrad
ymar2 <- (xmar2 / xrange) * pin[1] / pin[2] * yrange
expr <- quote({
    rect(xleft, usr[4] - 1.5 * ysize, xright, usr[4] - ymar,
         col = rainbow(8), border = NA)
    symbols(rep(usr[2] - xmar2, nrad), rep(usr[3] + ymar2, nrad),
            circles = Rad, inches = FALSE, add = TRUE, fg = NA,
            bg = gray.colors(nrad + 1, 1, 0)[-1])
    points(y)
    lines(spl)
})
rasterPlot(expr, res = 50)
box()
axis(1)
axis(2)

## The same picture with higher resolution but no antialiasing
plot(y, type = "n", axes = FALSE, ann = FALSE)
## No content in margin, but region = "figure" and bg = "white"
## paints margin white
rasterPlot(expr, antialias = "none", interpolate = FALSE,
           region = "figure", bg = "white")
## Draw box, axes, labels
parnew <- par(new = TRUE)
plot(x, y, type = "n", ylab = ylab)
par(parnew)

## Draw plot(1:5) with adjusted margins and additional axes.  Some parts
## are drawn with rasterPlot, others normally.  Resize to see stretching.
op <- par(no.readonly = TRUE)
par(mar = c(5.1, 4.1, 2.1, 2.1))
plot(1:5, type = "n", axes = FALSE, ann = FALSE)
expr2 <- quote({
    points(c(2, 4), c(2, 4))
    axis(2)
    axis(3)
})
rasterPlot(expr2, region = "figure", bg = "white")
points(c(1, 3, 5), c(1, 3, 5))
box()
axis(1)
axis(4)
title(xlab = "Index", ylab = "1:5")
par(op)

Regional Curve Standardization

Description

Detrend multiple ring-width series simultaneously using a regional curve.

Usage

rcs(rwl, po = NULL, nyrs = NULL, f = 0.5, biweight = TRUE, ratios = TRUE,
    rc.out = FALSE, make.plot = TRUE, method = c("caps", "ads"),
    min.n = NULL, pos.slope = TRUE, ...)

Arguments

Details

This method detrends and standardizes tree-ring series by calculating an age-related growth curve specific to the rwl. The detrending is the estimation and removal of the tree’s natural biological growth trend. The standardization is done by either dividing each series by the growth trend or subtracting the growth trend from each series to produce units in the dimensionless ring-width index (RWI). The option to produce indices by subtraction is intended to be used on series that have been subject to variance stabilization (e.g., using powt).

Two smoothers are available for fitting the regional curve. The default (method = "caps") uses a cubic smoothing spline where the frequency response is 0.50 at a wavelength of 10 percent of the maximum cambial age unless specified differently via nyrs and f (see caps). The alternative (method = "ads") uses an age-dependent spline whose stiffness increases with cambial age, which may better preserve low-frequency variability in the biological trend (see ads).

The min.n argument can be used to truncate the regional curve where sample depth becomes thin, preventing poorly replicated tail ages from influencing the fit.

This attempts to remove the low frequency variability that is due to biological or stand effects. See the references below for further details on detrending in general, and Biondi and Qeadan (2008) for an explanation of RCS.

A series with no values at all is dropped, with a message naming it, and the rest are standardized as if it had never been there.

Value

An object of class c("rwi", "data.frame") (see as.rwi) containing the dimensionless and detrended ring-width indices with column names, row names and dimensions of rwl. If rc.out is TRUE then a list will be returned with a data.frame containing the detrended ring widths as above and a vector containing the regional curve.

Author(s)

Code provided by DendroLab based on programming by F. Qeadan and F. Biondi, University of Nevada Reno, USA and adapted for dplR by Andy Bunn. Patched and improved by Mikko Korpela. Extended with method, min.n, and pos.slope arguments by Andy Bunn.

References

Biondi, F. and Qeadan, F. (2008) A theory-driven approach to tree-ring standardization: Defining the biological trend from expected basal area increment. Tree-Ring Research, 64(2), 81–96.

Cook, E. R. and Kairiukstis, L. A., editors (1990) Methods of Dendrochronology: Applications in the Environmental Sciences. Springer. ISBN-13: 978-0-7923-0586-6.

Fritts, H. C. (2001) Tree Rings and Climate. Blackburn. ISBN-13: 978-1-930665-39-2.

See Also

detrend, chron, cms, caps, ads

Examples

library(utils)
data(gp.rwl)
data(gp.po)

## Basic use: caps smoother (default), return indices only
gp.rwi <- rcs(rwl = gp.rwl, po = gp.po, biweight = TRUE,
              make.plot = FALSE)
str(gp.rwi)

## Return the regional curve alongside the indices
gp.out <- rcs(rwl = gp.rwl, po = gp.po, biweight = TRUE,
              rc.out = TRUE, make.plot = FALSE)
str(gp.out)

## Plot the regional curve with a title
rcs(rwl = gp.rwl, po = gp.po, biweight = TRUE,
    rc.out = FALSE, make.plot = TRUE, main = "Regional Curve")

## Use the age-dependent spline smoother instead of caps
gp.rwi.ads <- rcs(rwl = gp.rwl, po = gp.po, biweight = TRUE,
                  method = "ads", rc.out = FALSE, make.plot = FALSE)

## Truncate the regional curve where fewer than 5 series contribute.
## Indices for rings beyond the truncated curve are returned as NA.
gp.out.trunc <- rcs(rwl = gp.rwl, po = gp.po, biweight = TRUE,
                    method = "ads", min.n = 5,
                    rc.out = TRUE, make.plot = FALSE)

## Compare the two smoothers visually
op <- par(no.readonly = TRUE)
par(mfrow = c(1, 2))
rcs(rwl = gp.rwl, po = gp.po, biweight = TRUE,
    method = "caps", make.plot = TRUE, main = "caps")
rcs(rwl = gp.rwl, po = gp.po, biweight = TRUE,
    method = "ads", make.plot = TRUE, main = "ads")
par(op)

## Compare caps and ads indices for one series
gp.caps <- rcs(rwl = gp.rwl, po = gp.po, biweight = TRUE,
               method = "caps", rc.out = FALSE, make.plot = FALSE)
gp.ads  <- rcs(rwl = gp.rwl, po = gp.po, biweight = TRUE,
               method = "ads",  rc.out = FALSE, make.plot = FALSE)
cor(gp.caps[[1]], gp.ads[[1]], use = "complete.obs")

## Use subtraction rather than division (e.g. after variance stabilization)
gp.rwi.sub <- rcs(rwl = gp.rwl, po = gp.po, biweight = TRUE,
                  ratios = FALSE, make.plot = FALSE)

Read DPL Compact Format Ring Width File

Description

This function reads in a DPL compact format file of ring widths.

Usage

read.compact(fname)

Arguments

Details

This function should be able to read files written by the Dendrochronology Program Library (DPL) in its compact format.

Value

An object of class c("rwl", "data.frame") with the series in columns and the years as rows. The series IDs are the column names and the years are the row names.

Author(s)

Mikko Korpela

See Also

read.rwl, read.tucson, read.tridas, read.fh, write.compact

Examples

## Not run: 
data(co021)
fname <- write.compact(rwl.df = co021,
                       fname = tempfile(fileext=".rwl"),
                       append = FALSE, prec = 0.001)
foo <- read.compact(fname)
str(foo)
str(co021)
all.equal(foo,co021)
unlink(fname)

## End(Not run)

Read Tucson Format Chronology File

Description

This function reads in a Tucson (decadal) format file of tree-ring chronologies (.crn).

Usage

read.crn(fname, header = NULL, encoding = getOption("encoding"),
         long = TRUE)

Arguments

Details

This reads in a standard crn file as defined according to the standards of the ITRDB at https://www.ncei.noaa.gov/pub/data/paleo/treering/treeinfo.txt. Despite the standards at the ITRDB, this occasionally fails due to formatting problems.

Value

A data.frame with each chronology in columns and the years as rows. The chronology IDs are the column names and the years are the row names. If the file includes sample depth that is included as the last column (samp.depth). The output class is class "crn" and "data.frame"

Author(s)

Andy Bunn. Patched and improved by Mikko Korpela.

Read Heidelberg Format Ring Width File

Description

This function reads in a Heidelberg (block or column) format file of ring widths (.fh).

Usage

read.fh(fname, BC_correction = FALSE, encoding = NULL)

Arguments

Details

This reads in a fh-file with ring widths in blocks (decadal format) or in columns (e.g., as with comment flags) as used by TSAP program. Chronologies or half-chronos in fh-format are not supported.

Value

An object of class c("rwl", "data.frame") with the series in columns and the years as rows. The keycodes are the column names and the years are the row names. Depending on metadata available in the input file, the following attributes may be present in the data.frame:

Author(s)

Christian Zang. New features and patches by Mikko Korpela and Ronald Visser.

References

Rinn, F. (2003) TSAP-Win User Reference Manual. Rinntech, Heidelberg. https://rinntech.info/products/tsap-win/.

See Also

read.rwl

Read Site-Tree-Core IDs

Description

These functions try to read site, tree, and core IDs from a rwl data.frame.

Usage

read.ids(rwl, stc = c(3, 2, 3), ignore.site.case = FALSE,
         ignore.case = FALSE, fix.typos = FALSE, typo.ratio = 5,
         use.cor = TRUE)

autoread.ids(rwl, ignore.site.case = TRUE, ignore.case = "auto",
             fix.typos = TRUE, typo.ratio = 5, use.cor = TRUE)

Arguments

The following parameters affect the handling of suspected typing errors. Some have different default values in read.ids and autoread.ids.

Details

Because dendrochronologists often take more than one core per tree, it is occasionally useful to calculate within vs. between tree variance. The International Tree Ring Data Bank (ITRDB) allows the first eight characters in an rwl file for series IDs but these are often shorter. Typically the creators of rwl files use a logical labeling method that can allow the user to determine the tree and core ID from the label.

Argument stc tells how each series separate into site, tree, and core IDs. For instance a series code might be "ABC011" indicating site "ABC", tree 1, core 1. If this format is consistent then the stc mask would be c(3, 2, 3) allowing up to three characters for the core ID (i.e., pad to the right). If it is not possible to define the scheme (and often it is not possible to machine read IDs), then the output data.frame can be built manually. See Value for format.

The function autoread.ids is a wrapper to read.ids with stc="auto", i.e. automatic detection of the site / tree / core scheme, and different default values of some parameters. In automatic mode, the names in the same rwl can even follow different site / tree / core schemes. As there are numerous possible encoding schemes for naming measurement series, the function cannot always produce the correct result.

With stc="auto", the site part can be one of the following.

• In names mostly consisting of numbers, the longest common prefix is the site part

• Alphanumeric site part ending with alphabet, when followed by numbers and alphabets

• Alphabetic site part (quite complicated actual definition). Setting ignore.case to "auto" allows the function to try to guess when a case change in the middle of a sequence of alphabets signifies a boundary between the site part and the tree part.

• The characters before the first sequence of space / punctuation characters in a name that contains at least two such sequences

These descriptions are somewhat general, and the details can be found in regular expressions inside the function. If a name does not match any of the descriptions, it is matched against a previously found site part, starting from the longest.

The following ID schemes are detected and supported in the tree / core part. The detection is done per site.

• Numbers in tree part, core part starts with something else

• Alphabets in tree part, core part starts with something else

• Alphabets, either tree part all lower case and core part all upper case or vice versa. For this to work, ignore.case must be set to "auto" or FALSE.

• All digits. In this case, the number of characters belonging to the tree and core parts is detected with one of the following methods.

  • If numeric tree parts were found before, it is assumed that the core part is missing (one core per tree).

  • It the series are numbered continuously, one core per tree is assumed.

  • Otherwise, try to find a core part as the suffix so that the cores are numbered continuously.

  If none of the above fits, the tree / core split of the all-digit names will be decided with the methods described further down the list, or finally with the fallback mechanism.

• The combined tree / core part is empty or one character. In this case, the core part is assumed to be missing.

• Tree and core parts separated by a punctuation or white space character

If the split of a tree / core part cannot be found with any of the methods described above, the prefix of the string is matched against a previously found tree part, starting from the longest. The fallback mechanism for the still undecided tree / core parts is one of the following. The first one is used if use.cor is TRUE, number two if it is FALSE.

1. Pairwise correlation coefficients are computed between all remaining series. Pairs of series with above median correlation are flagged as similar, and the other pairs are flagged as dissimilar. Each possible number of characters (minimum 1) is considered for the share of the tree ID. The corresponding unique would-be tree IDs determine a set of clusterings where one cluster is formed by all the measurement series of a single tree. For each clustering (allocation of characters), an agreement score is computed. The agreement score is defined as the sum of the number of similar pairs with matching cluster number and the number of dissimilar pairs with non-matching cluster number. The number of characters with the maximum agreement is chosen.

2. If the majority of the names in the site use k characters for the tree part, that number is chosen. Otherwise, one core per tree is assumed. Parameter typo.ratio has a double meaning as it also defines what is meant by majority here: at least typo.ratio / (typo.ratio + 1) * n.tot, where n.tot is the number of names in the site.

In both fallback mechanisms, the number of characters allocated for the tree part will be increased until all trees have a non-zero ID or there are no more characters.

Suspected typing errors will be fixed by the function if fix.typos is TRUE. The parameter typo.ratio affects the eagerness to fix typos, i.e. the number of counterexamples required to declare a typo. The following main typo fixing mechanisms are implemented:

Site IDs.

If a rare site string resembles an at least typo.ratio times more frequent alternative, and if fixing it would not create any name collisions, make the fix. The alternative string must be unique, or if there is more than one alternative, it is enough if only one of them is a look-alike string. Any kind of substitution in one character place is allowed if the alternative string has the same length as the original string. The alternative string can be one character longer or one character shorter than the original string, but only if it involves interpreting one digit as the look-alike alphabet or vice versa. There are requirements to how long a site string must be in order to be eligible for replacement / typo fixing, i.e. cannot be shortened to zero length, cannot change the only character of a site string. The parameters ignore.case and ignore.site.case have some effect on this typo fixing mechanism.

Tree and core IDs.

If all tree / core parts of a site have the same length, each character position is inspected individually. If the characters in the i:th position are predominantly digits (alphabets), any alphabets (digits) are changed to the corresponding look-alike digit (alphabet) if there is one. The look-alike groups are {0, O, o}, {1, I, i}, {5, S, s} and {6, G}. The parameter typo.ratio determines the decision threshold of interpreting the type of each character position as alphabet (digit): the ratio of alphabets (digits) to the total number of characters must be at least typo.ratio / (typo.ratio + 1). If a name differs from the majority type in more than one character position, it is not fixed. Also, no fixes are performed if any of them would cause a possible monotonic order of numeric prefixes to break.

The function attempts to convert the tree and core substrings to integral values. When this succeeds, the converted values are copied to the output without modification. When non-integral substrings are observed, each unique tree is assigned a unique integral value. The same applies to cores within a tree, but there are some subtleties with respect to the handling of duplicates. Substrings are sorted before assigning the numeric IDs.

The order of columns in rwl, in most cases, does not affect the tree and core IDs assigned to each series.

Value

A data.frame with column one named "tree" giving an ID for each tree and column two named "core" giving an ID for each core. The original series IDs are copied from rwl as rownames. The order of the rows in the output matches the order of the series in rwl. If more than one site is detected, an additional third column named "site" will contain a site ID. All columns have integral valued numeric values.

Author(s)

Andy Bunn (original version) and Mikko Korpela (patches, stc="auto", fix.typos, etc.).

See Also

rwi.stats, read.rwl

Examples

library(utils)
data(ca533)
read.ids(ca533, stc = c(3, 2, 3))
autoread.ids(ca533)

Read Ring Width File

Description

This function reads in a file of ring widths (.rwl) in one of the available formats.

Usage

read.rwl(fname,
         format = c("auto", "tucson", "compact", "tridas", "heidelberg",
                    "sheet", "csv"),
         ...)

Arguments

Details

This is a simple wrapper to the functions actually implementing the read operation.

With format = "auto" the file's first lines are scanned and the most specific match wins: a DPL compact column specification or a HEADER: line, then a <tridas> element, then a NOAA/NCEI template header, then a comma separated sheet, and Tucson if nothing else matches. A sheet is recognised by every line splitting into the same number of comma separated fields, not by the presence of a comma. Only commas are looked for: a tab separated file may be a NOAA template, so tab separated sheets need format = "sheet" with a sep argument.

format = "sheet" is handled by read.sheet, which as of dplR 1.8.0 replaces the deprecated csv2rwl. Files that csv2rwl accepted may now be refused; see read.sheet.

"csv" is accepted as a synonym for "sheet" because it has been a valid value of format since before read.sheet existed. It is the less accurate of the two: read.sheet also reads tab, semicolon and pipe separated files, so format = "sheet" with a sep argument is an ordinary call. Prefer "sheet" in new code.

Value

If a "tucson", "compact", "heidelberg", "sheet" file is read (even through "auto"), returns an object of class c("rwl", "data.frame") with the series in columns and the years as rows. The series IDs are the column names and the years are the row names.

If a "tridas" file is read (even through "auto"), returns a list of results. See read.tridas for more information.

Author(s)

Mikko Korpela

See Also

read.tucson, read.tucson.legacy, read.compact, read.tridas, read.fh, read.sheet, write.rwl

Read Ring Widths from a Spreadsheet File

Description

Read a file of ring widths saved in “spreadsheet” layout: years down the rows, series across the columns, and the years in the first column.

Usage

read.sheet(fname, sep = NULL, dec = ".", layout = c("wide", "long"),
           transpose = FALSE, comment.char = "#", fill.internal.NA = NULL,
           fix.dup.char = "X", encoding = NULL, verbose = TRUE,
           strict = FALSE, ...)

csv2rwl(fname, ...)

Arguments

Details

The expected layout is the one a spreadsheet produces: the first column holds the years, the first row holds the series IDs, and each remaining cell holds one ring width.

An empty cell, NA, NaN, "." and "-" are all read as missing. Any other value that does not parse as a number is an error rather than a missing value, and the offending cells are named.

Unlike a Tucson file, a spreadsheet has no format to conform to, so this function checks the contents rather than assuming them. The years must be whole numbers, ascending, unique and consecutive: a year with no measurements is a row of NA, not a missing row. Every column after the first must be numeric. Series IDs are read verbatim, and anything that changes one – a stripped byte order mark, a resolved duplicate – is reported and recorded rather than done silently.

Separators

Comma, tab, semicolon and pipe separated files are all read. Left at NULL, sep is detected by looking for the candidate that gives every line the same number of fields, preferring the one that gives the most; quoted spans are ignored while counting, so a series ID may contain the separator. The detected separator is recorded in the provenance record as SEP_GUESS and printed when verbose is TRUE, but it is not a warning: detection is the default, and a defect it is not.

A semicolon separated file is usually a European export and carries decimal commas. When the separator is detected as a semicolon and a digit-comma-digit appears in the data, dec is set to "," and recorded as DEC_COMMA. Thousands separators are not handled: "1.234,5" and "1,234.5" are the same characters under two conventions and nothing in the file says which.

sep and dec may not be the same character.

Long format

With layout = "long" the file holds one row per observation rather than a rectangle: a series ID, a year and a value. Columns are taken by position, which is what write.sheet writes; if the header names all three of series, year and value, in any case and any order, the names are used instead.

The year span of the result runs from the earliest to the latest year in the file. A year with no observation in any series therefore comes back as a row of NA, and is reported, since years with no data break several dplR functions. One series may hold only one value per year.

read.sheet refuses a NOAA/NCEI template file. Those files are a commented metadata block followed by a tab separated data table, and reading the table alone would discard the coordinates, species, investigators and DOI above it.

Value

An object of class c("rwl", "data.frame") with the series in columns and the years as rows. The series IDs are the column names and the years are the row names.

The object carries an attribute "dplR.provenance", a list in the same shape as the one attached by read.tucson, recording what the reader found and what it changed. Its events element is a data.frame of findings, each with a stable id such as "ID_RENAMED" or "UNITS_SUSPECT".

Deprecated

csv2rwl is deprecated as of dplR 1.8.0 and now calls read.sheet. It will be made defunct in a later release.

This is a change of behaviour and not only a change of name. csv2rwl assigned the class "rwl" directly, which bypassed as.rwl and its requirement that the row names be consecutive integers. It therefore accepted files that read.sheet refuses:

• non-consecutive or duplicated years, which produced an object whose time was not a continuous span;

• columns holding text, which travelled inside the rwl until some later function failed on them;

• series IDs were passed through read.table(check.names = TRUE), so "1A" became "X1A" and "LF-2B" became "LF.2B", silently.

A file that csv2rwl read may therefore now stop with an error. That is the deprecation doing its job: the object it used to return was invalid. csv2rwl also stopped when a series ID appeared twice, where read.sheet renames the duplicate and records the change, as read.tucson does.

There is no csv2rwl.legacy. read.tucson.legacy exists because the behaviour of the old Tucson reader differs from the new one in ways that are occasionally wanted; that is not the case here.

Author(s)

Andy Bunn

See Also

write.sheet, read.rwl, read.tucson, rwl.check, as.rwl

Examples

library(utils)
data(ca533)
# write out a sheet that read.sheet will understand
tm <- time(ca533)
foo <- data.frame(tm, ca533, check.names = FALSE)
names(foo)[1] <- "Year"
tmpName <- tempfile(fileext = ".csv")
write.csv(foo, file = tmpName, row.names = FALSE)

# read it back in
bar <- read.sheet(tmpName, verbose = FALSE)
# the provenance record travels with the data
str(attr(bar, "dplR.provenance")$events)

unlink(tmpName)

Read Tree Ring Data Standard (TRiDaS) File

Description

This function reads in a TRiDaS format XML file. Measurements, derived series and various kinds of metadata are supported.

Usage

read.tridas(fname, ids.from.titles = FALSE,
            ids.from.identifiers = TRUE, combine.series = TRUE,
            trim.whitespace = TRUE, warn.units = TRUE)

Arguments

Details

The Tree Ring Data Standard (TRiDaS) is described in Jansma et. al (2010).

The parameters used for rearranging (ids.from.titles, ids.from.identifiers) and combining (combine.series) measurement series only affect the four lowest levels of document structure: element, sample, radius, measurementSeries. Series are not reorganized or combined at the upper structural levels (project, object).

Value

A list with a variable number of components according to the contents of the input file. The possible list components are:

Note

This is an early version of the function. Bugs are likely to exist, and parameters and return values are subject to change. Not all metadata defined in the TRiDaS specification is supported – unsupported elements are quietly ignored.

Author(s)

Mikko Korpela

References

Jansma, E., Brewer, P. W., and Zandhuis, I. (2010) TRiDaS 1.1: The tree-ring data standard. Dendrochronologia, 28(2), 99–130.

See Also

read.rwl, read.tucson, read.compact, read.fh, write.tridas

Read Tucson Format Ring Width File

Description

This function reads in a Tucson (decadal) format file of ring widths (.rwl).

Usage

read.tucson(fname, header = NULL, long = FALSE,
            encoding = getOption("encoding"),
            edge.zeros = TRUE, verbose = TRUE,
            comment.char = "#", fix.dup.char = "X",
            fill.internal.NA = NULL, strict = FALSE)

Arguments

Details

This reads in a standard rwl file as defined according to the standards of the ITRDB at https://www.ncei.noaa.gov/pub/data/paleo/treering/treeinfo.txt.

This function was replaced in dplR 1.8.0. The previous implementation is still available, unchanged, as read.tucson.legacy. The two differ in at least three important ways.

Interior gaps are no longer filled with zero. This is the change most likely to alter your numbers. Where a file records no measurement for a stretch of years inside a series, the old reader silently wrote zeros. A ring width of zero means a locally absent ring, which is a real observation about a tree in a year; a gap means the ring was not measured, which is a statement about the data. These are different claims and the reader no longer substitutes one for the other. Set fill.internal.NA = 0 to recover the old values exactly.

Note that several dplR functions, detrend among them, do not accept NA inside a series and will fail on such data. If you need to detrend a series containing interior gaps you must decide what those gaps mean and fill them deliberately, with fill.internal.NA or with fill.internal.NA = 0 here.

Malformed files are reported rather than guessed at. Where the ten fixed-width columns and a whitespace split of the same line disagree, the line does not conform to the format and nothing in it says which reading was meant, so it is returned as NA with both readings reported. Repeated series IDs are reported instead of being silently welded into one series, and are renamed when their records overlap. Tab characters, which have no defined width in a fixed-width format, are expanded to 8-column tab stops and reported.

The header and column layout are determined per line. This is why header and long are no longer used. header could only skip a fixed three lines and so could not describe the 4- and 6-line headers that real files contain. long is worse: years before -999 need five columns and so take column 8 away from the series ID, but that varies line by line, so a file mixing 8-character IDs with BC dates is read wrongly whichever way a per-file switch is set. Both are now decided from the content of each line. The arguments are still accepted, in their original positions, so that existing code and read.rwl's pass-through of ... do not fail with an “unused argument” error; passing either one warns.

Editorial: A closing word on the format itself. The Tucson decadal format dates from the era of the punch card. It is nominally fixed-width, yet files in the ITRDB arrive with tab characters sitting in the measurement fields, and a tab has no defined width at all. Column 8 belongs to the series ID unless the year is earlier than -999, in which case the year claims it, and nothing but a minus sign distinguishes the two readings. A negative number marks missing data, except when it marks the end of a record; 999 ends a record in a 0.01 mm file, so a genuine 9.99 mm ring cannot be told apart from a terminator. A single ID may carry several records written in any order at all: in ITRDB ‘ut542.rwl’ a series runs to 1841 and then starts over at 1353 further down the file. A decade line may resume off its boundary and quietly swallow the years it skipped, as ‘chil012.rwl’ does at 1242 for series AGU034. None of this is the fault of the people who made these data available, and the tree-ring community deserves credit for sharing its measurements openly long before most of science thought to do so. But we are still paying for a format built around eighty columns of cardboard. read.tucson tries to read what the file says, to report plainly what it did, and to refuse rather than guess when a line admits two readings. Even so, and despite our best efforts, this format can and does fail.

Value

An object of class c("rwl", "data.frame") with the series in columns and the years as rows. The series IDs are the column names and the years are the row names. Columns are returned in the order the series appear in the file.

The returned object carries a "dplR.provenance" attribute recording what the reader saw: the header lines, the precision each series was measured at and whether the file mixes them, any series the reader renamed (as old and new), the interior gaps and what the file held at each, and one row per recoverable problem found while parsing. Most of this cannot be recovered from the returned data or by reading the file again, and rwl.check reports on it. Subsetting an rwl object with [.rwl carries the record along and cuts it down to the series and years that are left, but functions that return a new object, such as detrend, drop it, so it describes the read rather than travelling with the data indefinitely.

Encoding

Ring widths are digits, but the header lines and occasionally the series IDs are not, and archived files carry accented site names and investigator names in whatever 8-bit encoding the contributor's machine used. read.tucson resolves this in three steps.

First, the file is tested against UTF-8. This is a decision and not a guess: UTF-8 is self-validating, and plain ASCII passes, so virtually every file takes this path and nothing is reported.

If the file is not valid UTF-8 and encoding was supplied, that encoding is used and an ENCODING_DECLARED event is recorded. If the declared encoding does not fit the file, read.tucson stops rather than transcode into replacement characters.

If the file is not valid UTF-8 and no encoding was supplied, it is read as "latin1" and an ENCODING_ASSUMED event is recorded, with a warning naming the line and quoting the text the assumption produced, so the guess can be judged at a glance. Under strict = TRUE this becomes an error.

The fallback guesses rather than running a charset detector because detectors do not work on this kind of file. They are byte-frequency language models and need a volume of non-ASCII text that a file of ring widths never has; on a rwl that is ASCII apart from one accented name, ICU's detector scores the right answer at 0.16 and ranks UTF-16 alongside it. "latin1" is used instead because it is total (every byte decodes, so it cannot fail) and lossless (the round trip is byte-identical, so a wrong guess discards nothing and the file can simply be re-read with encoding set).

Author(s)

Original reader by Hung Nguyen. Reworked and extended by Andy Bunn.

See Also

read.tucson.legacy, read.rwl, read.compact, read.tridas, read.fh, write.tucson, fill.internal.NA, [.rwl

Examples

## A short file with an interior gap: the 1910s are marked -999,
## i.e. not measured.
tf <- tempfile()
writeLines(c("TST01A  1900   123   134   145   156   167   178   189   150   141   132",
             "TST01A  1910  -999  -999  -999  -999  -999  -999  -999  -999  -999  -999",
             "TST01A  1920   211   222   233   244   255   266   277   288   299   300",
             "TST01A  1930   311   999"), tf)

## The default reports the gap and returns it as NA.
x <- read.tucson(tf)
x["1915", ]

## Ask for the old behaviour explicitly.
x0 <- read.tucson(tf, fill.internal.NA = 0)
x0["1915", ]

unlink(tf)

Read Tucson Format Ring Width File (Legacy Reader)

Description

This function reads in a Tucson (decadal) format file of ring widths (.rwl). It is the reader that dplR used through version 1.7.9, kept unchanged for backward compatibility. New code should use read.tucson.

Usage

read.tucson.legacy(fname, header = NULL, long = FALSE,
                   encoding = getOption("encoding"),
                   edge.zeros = TRUE, verbose = TRUE)

Arguments

Details

This reads in a standard rwl file as defined according to the standards of the ITRDB at https://www.ncei.noaa.gov/pub/data/paleo/treering/treeinfo.txt. Despite the standards at the ITRDB, this occasionally fails due to formatting problems.

This function was the default Tucson reader in dplR through version 1.7.9. As of 1.8.0 that role belongs to read.tucson, which reads files this one cannot, reports malformed lines rather than guessing at them, and by default returns interior gaps as NA instead of filling them with zero. This function is retained so that existing workflows continue to produce exactly the numbers they always have, and because it honours the encoding argument, which read.tucson does not. Use read.tucson(fname, fill.internal.NA = 0) to get the new reader's parsing with the old reader's treatment of gaps.

Value

An object of class c("rwl", "data.frame") with the series in columns and the years as rows. The series IDs are the column names and the years are the row names.

Author(s)

Andy Bunn. Patched and greatly improved by Mikko Korpela.

See Also

read.tucson, read.rwl, read.compact, read.tridas, read.fh, write.tucson

Red-Noise Spectra of Time-Series

Description

Estimate red-noise spectra from a possibly unevenly spaced time-series.

Usage

redfit(x, t, tType = c("time", "age"), nsim = 1000, mctest = TRUE,
       ofac = 4, hifac = 1, n50 = 3, rhopre = NULL,
       p = c(0.10, 0.05, 0.02), iwin = 2,
       txOrdered = FALSE, verbose = FALSE, seed = NULL,
       maxTime = 10, nLimit = 10000)

runcrit(n, p = c(0.10, 0.05, 0.02), maxTime = 10, nLimit = 10000)

Arguments

Details

Function redfit computes the spectrum of a possibly unevenly sampled time-series by using the Lomb-Scargle Fourier transform. The spectrum is bias-corrected using spectra computed from simulated AR1 series and the theoretical AR1 spectrum.

The function duplicates the functionality of program REDFIT by Schulz and Mudelsee. See the manual of that program for more information. The results of this function should be very close to REDFIT. However, some changes have been made:

• More precision is used in some constants and computations.

• All the data are used: the last segment always contains the last pair of (t, x). There may be small differences between redfit and REDFIT with respect to the number of points per segment and the overlap of consecutive segments.

• The critical values of the runs test (see the description of runcrit below) differ between redfit and REDFIT. The approximate equations in REDFIT produce values quite far off from the exact values when the number of frequencies is large.

• The user can select the significance levels of the runs test.

• Most of the window functions have been adjusted.

• 6 dB bandwidths have been computed for discrete-time windows.

Function runcrit computes the limits of the acceptance region of a number of runs test: assuming a sequence of n i.i.d. discrete random variables with two possible values a and b of equal probability (0.5), we are examining the distribution of the number of runs. A run is an uninterrupted sequence of only a or only b. The minimum number of runs is 1 (a sequence with only a or only b) while the maximum number is n (alternating a and b). See Bradley, p. 253–254, 259–263. The function is also called from redfit; see rcnt in ‘Value’ for the interpretation. In this case the arguments p, maxTime and nLimit are passed from redfit to runcrit, and n is the number of output frequencies.

The results of runcrit have been essentially precomputed for some values of p and n. If a precomputed result is not found and n is not too large (nLimit, maxTime), the exact results are computed on-demand. Otherwise, or if package "gmp" is not installed, the normal distribution is used for approximation.

Value

Function runcrit returns a list containing rcritlo, rcrithi and rcritexact (see below). Function redfit returns a list with the following elements: