Introduction
Coding can be a major hurdle to conducting data analysis in
dendrochronology. This workbook is meant to get you up and running
quickly. This script will allow you to enter your data and start
analyzing. I’ve set up the code so that you can load in your own data
and run many common analyses. But first, let’s run through some key
packages and functions. If you are new to R, I recommend checking out Essential
R by my colleague Eric Nord. If you are new to dendrochronological
analysis in R, I recommend checking out OpenDendro by my
colleagues Andy Bunn, Kevin Anchunkaitis, Ed Cook, and Tyson Swetnam. If
you want step-by-step video help, check out my website The
Treeringist for videos on how to use common R packages, CooRecorder,
Cdendro, and links to dendro stuff.
Libraries Used in this Document
There are a number of useful libraries for manipulating and analyzing
tree-ring data in R. In this document we will use the following
tree-ring specific and general libraries. You can use
citation("x") to find the proper citation for each package
where
x is the name of the package that you used.
- dplR
- treeclim
- TRADER
- DendroSync
- dendroTools
- burnr
- defoliatR
- graphics
- ggplot2
- utils
- dplyr
- tidyr
- tidyverse
Installing and Loading Libraries
If you do not have one or many of these libraries installed you can
add them by using the install.packages("x") function where
x is the name of the package you desire to install. Be sure to
put the name in parentheses. Start by loading all of the
libraries:
library(dplR)
library(treeclim)
library(TRADER)
library(DendroSync)
library(dendroTools)
library(burnr)
library(dfoliatR)
library(graphics)
library(utils)
library(ggplot2)
library(dplyr)
library(tidyr)
library(tidyverse)
Crossdating in dplR
Using the script below, you will be able to read in raw ring width
files for crossdating or COFECHA-like analysis. To run your own
data, change the fname to match the data in your
grow.rwl <- read.rwl(fname = ... , format = "auto")
script. Alternatively you can use the RWL file included in this
project.
Note: The format of most raw ring width files is the Tucson format
that organizes measurements in a decadal format. You might have
measurements formatted in columns. If you do, try using
grow.rwl <- read.table(fname = ...). THe
head() function shows the top 10 lines of the data. I use
this just to make sure it has been read in correctly. “va013.rwl” is a
red spruce (Picea rubens) chronology collected by Ed Cook back
in the 1980s near Mountain Lake Biological Station in Virginia.
grow.rwl <- read.tucson(fname = "va013.rwl")
## There does not appear to be a header in the rwl file
## [1] 0
## There are 41 series
## 1 089011 1752 1982 0.01
## 2 089012 1725 1982 0.01
## 3 089031 1725 1982 0.01
## 4 089032 1747 1982 0.01
## 5 089041 1835 1982 0.01
## 6 089042 1813 1982 0.01
## 7 089051 1795 1982 0.01
## 8 089052 1790 1982 0.01
## 9 089061 1800 1982 0.01
## 10 089062 1800 1982 0.01
## 11 089071 1804 1982 0.01
## 12 089072 1816 1982 0.01
## 13 089112 1780 1982 0.01
## 14 089121 1792 1982 0.01
## 15 089131 1748 1982 0.01
## 16 089132 1728 1982 0.01
## 17 089141 1811 1982 0.01
## 18 089151 1694 1982 0.01
## 19 089152 1742 1982 0.01
## 20 089161 1707 1970 0.01
## 21 089162 1709 1982 0.01
## 22 089171 1741 1961 0.01
## 23 089172 1782 1982 0.01
## 24 089181 1826 1982 0.01
## 25 089191 1734 1982 0.01
## 26 089192 1778 1982 0.01
## 27 089201 1709 1982 0.01
## 28 089202 1750 1982 0.01
## 29 089211 1728 1982 0.01
## 30 089412 1726 1982 0.01
## 31 089421 1743 1982 0.01
## 32 089422 1749 1982 0.01
## 33 089431 1750 1982 0.01
## 34 089432 1747 1982 0.01
## 35 089441 1771 1982 0.01
## 36 089442 1766 1982 0.01
## 37 089451 1820 1982 0.01
## 38 089452 1801 1982 0.01
## 39 089461 1812 1982 0.01
## 40 089471 1777 1982 0.01
## 41 089472 1784 1982 0.01
head(grow.rwl)
You can examine summary statistics of your ring width file using the
rwl.stats(x) function.
rwl.stats(grow.rwl)
You can plot the series to examine sample depth and observe the raw
ring widths.
seg.plot(grow.rwl)
Notice the periodic growth releases seen in the raw ring width data.
This indicates some disturbance is present at the site. This is common
in closed canopy ecosystems.
spag.plot(grow.rwl)
To assess your crossdating, print out a report of the rwl file. The
mean series intercorrelation is a key metric to assess the quality of
crossdating.
rwl.report(grow.rwl)
## Number of dated series: 41
## Number of measurements: 8863
## Avg series length: 216.1707
## Range: 289
## Span: 1694 - 1982
## Mean (Std dev) series intercorrelation: 0.4898686 (0.06904016)
## Mean (Std dev) AR1: 0.8795122 (0.04203399)
## -------------
## Years with absent rings listed by series
## Series 089011 -- 1914
## Series 089162 -- 1915
## 2 absent rings (0.023%)
## -------------
## Years with internal NA values listed by series
## None
COFECHA
You can perform a number of crossdating functions in dplR. If you are
like me, you might prefer Andy Bunn’s Shiny app called XdateR. I
primarily use the app for crossdating but using the code below helps you
get under the hood a bit more. Or go old school and rock COFECHA in the
command prompt.
corr.rwl.seg(rwl = grow.rwl, seg.length = 50, bin.floor = 100, n = NULL, prewhiten = TRUE, pcrit = 0.05, biweight = TRUE, method = c("spearman"), make.plot = TRUE, label.cex = 1, floor.plus1 = FALSE,
master = NULL)
You can examine an individual core with this function. Just replace
the series name.
series.rwl.plot(grow.rwl, series = "089472", series.yrs = as.numeric(names(series)),
seg.length = 50, bin.floor = 100, n = NULL,
prewhiten = TRUE, biweight = TRUE, floor.plus1 = FALSE)
The series intercorrelation is a key statistic that you will need to
review to validate your crossdating of each sample. Values are dependent
on the species and location of sampling.
interseries.cor(grow.rwl, n = NULL, prewhiten = TRUE, biweight = TRUE, method = "spearman")
If you need assistance with visual crossdating, you might try loading
in a well-dated chronology and using the pointer function
to find narrow rings.
markers <- pointer(grow.rwl)#Calculate marker rings from raw ring width series
tail(markers)
Detrending and Standardization
ARSTAN
ARSTAN is the brain child of Dr. Ed Cook. It stands
for autoregressive standardization. Birthed from his 1985
dissertation, Ed laid the modern groundwork for detrending and
standardizing tree-ring (i.e., time series) data. Buy him an IPA if you
run into him.
The standalone ARSTAN program allows you to
interactively detrend and standardize a tree-ring series and develop
site chronologies. You can download ARSTAN here
if you want the full variety of options. But please be cautious. You can
do a lot of things to your data that might not be appropriate.
dplR and Detrending
Many ARSTAN options are also baked into
dplR and it is constantly being updated by some
dedicated colleagues. It is best to explore your individual tree-ring
series when you work at a new site and/or with a new species. To run
interactive detrending you can use the i.detrend()
function. This allows you to explore curve fits for each tree ring
series. The main idea of detrending is the removal of age-related growth
trends and disturbance signals that obscure the climate-related signal.
Andy Bunn’s iDetrend Shiny
App is a great way to use this code to explore detrending.
If you want to run interactive detrending in Rstudio, use this code:
grow.rwi.int <- i.detrend(rwl = grow.rwl, nyrs = NULL, f =
0.5,pos.slope = FALSE)
You will need to select a detrending method for each series. You
could use the “grow.rwi.int” data frame in the rest of the code below if
you’d like. I typically use only 1 or 2 detrending methods for a given
site and species. This is easier to justify and replicate.
If you want to utilize a singular detrending method for all of your
series after you used interactive detrending to determine the best fit
for your data, you can use the detrend() function instead
of i.detrend(). In this function, you will choose one of
the following detrending options with method = : “Spline”,
“ModNegExp”, “Mean”, “Ar”, “Friedman”, “ModHugershoff” or
“AgeDepSpline”.
Let’s detrend the red spruce raw ring widths using the Age Dependent
Spline that allow for some flexibility when the tree is young and gets
more stiff as the tree ages. This is a good first method to try. The
results are time series of indexed values with a mean of 1.
grow.rwi <- detrend(rwl = grow.rwl, method = c("AgeDepSpline"), nyrs = NULL, f = 0.5, pos.slope = FALSE)
tail(grow.rwi)
To view the spaghetti plot for the detrended series you can once
again use the spag.plot() function similar to above except
substituting the raw ring width data for the detrended data.
spag.plot(rwl = grow.rwi, zfac = 1, useRaster = FALSE, res = 300)
Notice some of the artifacts of the age dependent spline towards the
end of the time series. You might consider using a spline (e.g., 50 year
spline) that is more flexible at the tails.
grow.rwi <- detrend(rwl = grow.rwl, method = c("Spline"), nyrs = 50, f = 0.5, pos.slope = FALSE)
tail(grow.rwi)
spag.plot(rwl = grow.rwi, zfac = 1, useRaster = FALSE, res = 300)
To examine the statistics for the entire chronology you can use
rwi.stats() or rwi.stats.running() to use
running statistics in order to adjust the time periods. See
?rwi.stats.running for help on this function.
rwi.stats(grow.rwi)
I like to look at 50 year windows with a 49 year overlap to determine
when the Expressed Population Signal (EPS) falls below 0.85, indicating
that the common signal in the chronology is weakening. This is common as
fewer older trees are used to calculate the mean. I called the object
“stat_out” so that I can review the data frame later. Otherwise, it is a
large print out.
stat_out <- rwi.stats.running(grow.rwi, method = c("spearman"), prewhiten = FALSE, window.length = 50, window.overlap = 49)
tail(stat_out)
dplR and Generating Chronologies
The ARSTAN program generates a standardized chronology (maintains
temporal autocorrelation), residual chronology (autoregressive model
applied), or the ARSTAN chronology (retains a portion of the pooled
autocorrelation). In dplR you can use the chron()
function to build a mean value chronology from detrended ring widths
produced from a detrend() function.
Let’s look at the standardized chronology that retains temporal
autocorrelation, or the effect of the previous years’ growth on the
current year. The biweight function is a method to calculate the mean
that is common in dendrochronology. It attempts to reduce the effects of
outlier values in the calculation of the mean.
grow.crn <- chron(x = grow.rwi, prefix = "MTL", biweight = TRUE, prewhiten = FALSE)
And make a figure of the chronology with a smoothing spline
highlighting decadal trends.
plot.crn(x = grow.crn, add.spline = TRUE, nyrs = 20)
Now, we can apply an autoregressive model to remove the
autocorrelation from the individual indexed series. This is also known
as prewhitening.
grow.crn <- chron(x = grow.rwi, prefix = "MTL", biweight = TRUE, prewhiten = TRUE)
And make a figure to compare the standardized vs. residual
chronologies.
plot.crn(x = grow.crn, add.spline = TRUE, nyrs = 20)
Finally, we have the ARSTAN chronlogy that applies Ed’s special sauce
in the chronology building process. Essentially, the ARSTAN chronology
removes the autocorrelation with an autoregressive model, calculates the
pooled autocorrelation across all series, and applies a portion of the
common autocorrelation back in the time series. It has the benefit of
retaining some autocorrelation that might be related to climate but
removing some of that biological signal.
## Pooled AR Summary
## ACF
## ar(0) ar(1) ar(2) ar(3) ar(4) ar(5)
## 1.00000000 0.58342621 0.30073586 0.12027204 -0.04652745 -0.16403166
## ar(6) ar(7) ar(8) ar(9) ar(10)
## -0.23512777 -0.22318104 -0.21156920 -0.18212902 -0.14112659
## AR Coefs
## 1 2 3 4 5 6
## ar(1) -0.5834262 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## ar(2) -0.6184968 0.06011136 0.00000000 0.00000000 0.00000000 0.00000000
## ar(3) -0.6156924 0.03125682 0.04665269 0.00000000 0.00000000 0.00000000
## ar(4) -0.6097196 0.03525851 -0.03217207 0.12802620 0.00000000 0.00000000
## ar(5) -0.5965402 0.03194661 -0.02854243 0.06525953 0.10294350 0.00000000
## ar(6) -0.5864574 0.03833843 -0.03133801 0.06838853 0.04451560 0.09794461
## ar(7) -0.5846967 0.03913869 -0.03010859 0.06782517 0.04520481 0.08740186
## ar(8) -0.5833798 0.04554083 -0.02679736 0.07279333 0.04299937 0.09026875
## ar(9) -0.5798381 0.04433920 -0.02243267 0.07487244 0.04651909 0.08897304
## ar(10) -0.5775063 0.04651138 -0.02352497 0.07916323 0.04876250 0.09258382
## 7 8 9 10
## ar(1) 0.00000000 0.00000000 0.00000000 0.00000000
## ar(2) 0.00000000 0.00000000 0.00000000 0.00000000
## ar(3) 0.00000000 0.00000000 0.00000000 0.00000000
## ar(4) 0.00000000 0.00000000 0.00000000 0.00000000
## ar(5) 0.00000000 0.00000000 0.00000000 0.00000000
## ar(6) 0.00000000 0.00000000 0.00000000 0.00000000
## ar(7) 0.01797701 0.00000000 0.00000000 0.00000000
## ar(8) -0.02485176 0.07324956 0.00000000 0.00000000
## ar(9) -0.02264976 0.04504190 0.04835214 0.00000000
## ar(10) -0.02373159 0.04718019 0.02038903 0.04822572
## AIC
## ar(0) ar(1) ar(2) ar(3) ar(4) ar(5) ar(6) ar(7) ar(8) ar(9) ar(10)
## 2557 2439 2440 2441 2439 2438 2437 2439 2439 2440 2442
## Selected Order
## [1] 1
And make a figure to compare the standardized, residual chronologies,
and ARSTAN chronologies.
plot.crn(x = grow.crn.ars, add.spline = TRUE, nyrs = 20)
Growth-Climate Analysis in Treeclim
The Treeclim package allows for the assessment of
growth-climate relationships. The statistical analysis is somewhat more
complex than a standard response or correlation function. Treeclim
was created by Dr. Christian Zang. Be sure to cite the publication if
you use the code. If you want to learn more about the analysis, read Biondi
and Waikul 2004. The package requires a chronology and monthly
climate data for your location. If you are in the contiguous U.S., you
can download interpolated monthly and daily climate data from the PRISM Climate Group.
Review your chronology data with the summary()
function.
summary(grow.crn)
## std res samp.depth
## Min. :0.3256 Min. :0.5186 Min. : 1.00
## 1st Qu.:0.9012 1st Qu.:0.9205 1st Qu.:21.00
## Median :0.9759 Median :0.9929 Median :40.00
## Mean :0.9845 Mean :0.9930 Mean :30.67
## 3rd Qu.:1.0758 3rd Qu.:1.0601 3rd Qu.:41.00
## Max. :1.4461 Max. :1.4301 Max. :41.00
## NA's :5
Here, we will read in PRISM climate data and format for treeclim. To
run with your own data, change file name.
clim <- read.table(file = "MTL_PRISM_ppt_tmean.csv", skip = 10, header = TRUE, sep = ",") #skips reading the header
head(clim)
#Transform data with custom function; h/t V. Harris
PRISM.ym <- function(PRISM.data){
PRISM.data <- PRISM.data %>% separate(Date, sep="-", into=c("Year", "Month"))
PRISM.data$Year <- as.numeric(PRISM.data$Year)
PRISM.data$Month <- as.numeric(PRISM.data$Month)
return(PRISM.data)
}
clim3 <- PRISM.ym(clim)
colnames(clim3) <- (c("Year", "Month", "PPT", "TMEAN"))
head(clim3)
You will need to decide which version of the chronology you would
like to test against climate data. Here, I pull out the residual
chronology for analysis.
grow.crn.res <- grow.crn[,-2]
Main treeclim function
You are required to call up your chronology, climate data, identify
your month selection (negative indicates the previous year), method
(i.e., response or correlation), and window size. I have used some
common settings below.
The dynamic argument can take dynamic = “static”, “moving”,
“evolving”. Static calculates coefficients on the maximum period of
overlap between growth and climate. Moving use a window of user defined
length and shifts it incrementally through time to assess the time
stability of the growth-climate relationship. Similarly, the evolving
dynamic assesses stability but it initiates with a window of user
defined length and then growths the window length to assess the
influence of incorporating additional climate data on the
relationship.
Here is a static analysis examining the correlation between ring
width and precipitation and mean temperature from the previous May
through the current October. The analysis first calculates the
correlation with precipitation. Then, it removes that influence and
calculates the correlation with mean temperature. Oh, and there is
bootstrapping to calculate confidence intervals. That’s too much to
explain here.
resp <- dcc(chrono = grow.crn.res, climate = clim3, selection = -5:10,
method = "correlation", dynamic = "static", win_size = 35, win_offset = 1, start_last = FALSE, timespan = NULL, var_names = NULL, ci = 0.05, boot = "stationary", sb = FALSE)
## Running for timespan 1896 - 1982...
plot(resp) #call the object resp to show coefficients
Here is a moving analysis using a 35-year window starting towards the
pith but I have adjusted the months to the current May through the
current October. You have to be careful to not violate the assumptions
of parametric regression. You will use all of your degrees of freedom
with an 18-month window with two variables moving through time. The
function will tell you when you messed up.
resp <- dcc(chrono = grow.crn.res, climate = clim3, selection = -5:-10,
method = "correlation", dynamic = "moving", win_size = 35, win_offset = 1, start_last = FALSE, timespan = NULL, var_names = NULL, ci = 0.05, boot = "stationary", sb = FALSE)
plot(resp) #call the object resp to show coefficients
Here is an evolving analysis that starts with a 35-year window
anchored in 1895. The window then grows to 36 years, 37 years, etc.
resp <- dcc(chrono = grow.crn.res, climate = clim3, selection = -5:-10,
method = "correlation", dynamic = "evolving", win_size = 35, win_offset = 1, start_last = FALSE, timespan = NULL, var_names = NULL, ci = 0.05, boot = "stationary", sb = FALSE)
plot(resp) #call the object resp to show coefficients
Seasonal Correlations
The following code for seasonal correlation analysis was developed by
Dr. Dave Meko in Matlab. It is commonly used to assess seasons of
varying lengths.
seas <- seascorr(grow.crn.res, clim3, var_names = NULL, timespan = NULL, complete = 9, season_lengths = c(1, 3, 6), primary = "PPT", secondary = "TMEAN", ci = 0.05)
## treeclim tries to use the maximum overlap in timespan for chronology and climate data. The overlap starts in 1895, but to be able to use climate data from the previous year(s) (as you chose by setting 'selection' accordingly), the analysis starts in 1896.
## Running for timespan 1896 - 1982...
plot(seas)
More advanced growth-climate analysis
You might have need to evaluate the skill of a reconstruction model.
Here, I demonstrate how to use a split calibration period to assess
model viability. The function requires 2 climate variables.
calib <- dcc(chrono = grow.crn.res, climate = clim3, selection = .mean(-6:-9, "TMEAN") + .sum(-6:-7, "PPT"), #use a selection with recon variable of interest - modifiers like .mean or .sum can be used to average across months
method = "response")
## treeclim tries to use the maximum overlap in timespan for chronology and climate data. The overlap starts in 1895, but to be able to use climate data from the previous year(s) (as you chose by setting 'selection' accordingly), the analysis starts in 1896.
## Running for timespan 1896 - 1982...
calib #show model results
## Coefficients (significance flags correspond to p < 0.05):
## id varname month coef
## PPT.sum.prev.jun.prev.jul 1 PPT.sum jun...jul 0.285
## TMEAN.mean.prev.jun...ul.prev.aug.prev.sep 2 TMEAN.mean jun...sep -0.285
## significant ci_lower ci_upper
## PPT.sum.prev.jun.prev.jul FALSE -0.156 0.381
## TMEAN.mean.prev.jun...ul.prev.aug.prev.sep TRUE -0.382 -0.062
plot(calib)
This does the model evaluation.
skillz <- skills(object = calib, model = "ols", calibration = "50%", timespan = NULL)
## Warning in skills(object = calib, model = "ols", calibration = "50%", timespan
## = NULL): Reconstruction skills cannot be evaluated when using more than one
## independent variable. We use only the first variable (by alphabet:
## PPT.sum.prev.jun.prev.jul).
plot(skillz)
skillz #show model results
## Call:
## skills(object = calib, model = "ols", calibration = "50%", timespan = NULL,
## target = .mean(-6:-9, "TMEAN") + .sum(-6:-7, "PPT"))
##
## Using climate target PPT.sum.prev.jun.prev.jul and calibration model with 50 percent (= 44 years) of data starting at recent end as calibration period:
##
## r p-value
## 0.1776091 0.1089109
##
## Coefficients:
## intercept slope
## 154.76902 65.89314
##
## Verification statistics:
## RE CE prediction RMSE
## 0.085 0.019 14.414
##
## Durbin-Watson Test: 1.105 (p = 0.000592928)
##
## Model for whole period:
##
## r p-value
## 0.26247049 0.00990099
##
## Coefficients:
## intercept slope
## 124.1532 104.9594
Growth Release Detection in TRADER
Most of the analysis above concentrated on honing the climate signal
while reducing the effect of biological trends and ecological noise. But
if you are an ecologist, you want to keep that signal related to stand
dynamics and disturbance. One method for understanding forest growth
dynamics relies on detecting growth releases. There are many methods to
detect releases. TRADER
was created by Dr. Jan Altman. Use ?TRADER to learn about
the techniques included in the package. I will run through the common
growth-averaging technique from Nowacki
and Abrams 1997.
rwl.stats(grow.rwl)
The growth-averaging technique requires raw ring width data so you
will use your original rwl file. The growthAveragingALL
function looks at windows before and after each year to determine if the
percent growth change has exceeded a defined threshold. I have used
common settings in the example. The function creates output for each
series and summary figures in your working directory.
growthAveragingALL(grow.rwl, releases = NULL, m1 = 15, m2 = 15,buffer = 10, drawing = TRUE, criteria = 0.25, criteria2 = 0.50,prefix = "ga", gfun = mean, length = 5, storedev = jpeg)
## [1] "## Nowacki & Abrams analysis!"
## [1] "Criteria 0.25 Criteria2 0.5 m1 15 m2 15 Buffer 10 Length 5"
## [1] "Total number of releases >= 0.25 & < 0.5 is 55"
## inyears
## 1748 1762 1771 1777 1780 1781 1785 1789 1792 1799 1811 1813 1820 1821 1829 1830
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3
## 1832 1833 1847 1853 1859 1860 1867 1879 1883 1896 1898 1899 1905 1907 1916 1922
## 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1
## 1924 1927 1928 1936 1937 1938 1940 1946 1947 1955 1956 1957 1967
## 1 3 1 1 1 2 1 1 2 1 2 1 3
## [1] "Total number of releases >= 0.5 is 52"
## inyears
## 1728 1742 1752 1760 1785 1787 1788 1789 1807 1808 1830 1832 1834 1836 1861 1905
## 1 2 1 1 1 1 1 1 1 1 1 1 2 1 1 2
## 1906 1907 1911 1921 1926 1927 1931 1934 1935 1937 1938 1939 1940 1946 1947 1948
## 2 1 1 1 1 1 1 1 1 4 1 1 1 3 2 1
## 1949 1950 1956 1957 1961 1965
## 2 2 1 3 1 1
Assessing growth synchronony with dendrosync
Tree growth should ideally show some synchronous growth across a site
and a region because climate is a primary factor controlling tree
growth, enabling dendrochronologists to crossdate. We can quantify the
level of synchrony using the dendrosync package. Dendrosync
was created by Josu Alday. “It combines variance-covariance mixed
modelling with functions that quantify the degree to which tree-ring
chronologies contain a common signal over a fixed time period. It also
estimates temporal changes in synchrony using a moving window
algorithm.” This function calculates the synchronous growth changes
(sgc), semi synchronous growth changes (ssgc) and the length of the
compared overlap for a given set of tree-ring records. The package
creates a covariance matrix comparing all series in the collection.
These could be individual trees as in the example or the time series
could be chronologies in a region.
changes <- sgc(grow.rwl,overlap = 50, prob = TRUE)
mean(changes$sgc_mat, na.rm = TRUE)
## [1] 0.5866611
mean(changes$ssgc_mat, na.rm = TRUE)
## [1] 0.05492645
Basal area increment calculation in dplR
Let’s start with a summary of the raw ring width data.
rwl.stats(grow.rwl) #summary and stats of raw ring width file
Next, let’s run the BAI function that goes from pith to bark
and plot the results. Finally, you can write the data to a csv file for
later use.
basal <- bai.out(grow.rwl, diam = NULL)
spag.plot(basal, zfac = 1, useRaster = FALSE, res = 300)
write.csv(basal, "bai.csv")
Fire history graphing in burnr
Fire history analyis relies on a list of events or fire dates. This
type of graphing is/was commonly done in the standalone programs FHX and
FHAES. FHAES works quite
well as a standalone program with a point and click interface. The
format of the files is somewhat unique so I highly recommend using
?burnr when using this code. The example here is simply
plotting fire scar data from Zion (h/t to Chris Gentry and Peter
Brown).
Zion <- read_fhx('Zion.fhx')
Sites <- read.csv('ZionSiteIDs.csv')
facetplot <- plot_demograph(Zion, facet_group = Sites$SiteID, facet_id = Sites$series, plot_legend = TRUE)
print(facetplot)
rugplot <- plot_demograph(Zion, composite_rug = TRUE, plot_legend = TRUE)
compositerug <- rugplot + annotate('rect', xmin = 1721, xmax = 1723, ymin = 0, ymax = 21, alpha = 0.4) + annotate('rect', xmin = 1734, xmax = 1736, ymin = 0, ymax = 21, alpha = 0.4) + annotate('rect', xmin = 1748, xmax = 1750, ymin = 0, ymax = 21, alpha = 0.4) + annotate('rect', xmin = 1777, xmax = 1779, ymin = 0, ymax = 21, alpha = 0.4) + annotate('rect', xmin = 1793, xmax = 1795, ymin = 0, ymax = 21, alpha = 0.4) + scale_x_continuous(limits=c(1450, 2005), breaks = seq(1450,2005,25)) + theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust = 0.2), panel.grid.major = element_blank(), panel.grid.minor = element_blank())
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
print(compositerug)
Detecting insect outbraks with dfoliatR
The Outbreak program was originally programmed by Dr. Tom Swetnam.
The defoliatR package with generously coded by Dr. Chris
Guiterman. This package will allow you to detect insect outbreak
events using host and non-host tree-ring chronologies. I’ve used an
example from the Intro to doliatR
Let’s load in some data.
data(dmj_h) #host data is Douglas fir with western spruce budworm damage.
data(dmj_nh) #non-host data is ponderosa pine from a nearby control site
Here is the main function for analysis.
dmj_defol <- defoliate_trees(host_tree = dmj_h, nonhost_chron = dmj_nh,
duration_years = 8, max_reduction = -1.28, list_output = FALSE)
plot(dmj_defol)
defol_stats(dmj_defol)
dmj_obr <- outbreak(dmj_defol, filter_perc = 25, filter_min_series = 3)
plot_outbreak(dmj_obr)
dmj_obr_stats <- outbreak_stats(dmj_obr)
head(dmj_obr_stats)
dmj_interv <- diff(dmj_obr_stats$start)
Output all outbreak return intervals, mean interval, and median
interval respectivtly.
dmj_interv# All return intervals
## [1] 73 72 24 32 78 28
mean(dmj_interv)# Mean interval
## [1] 51.16667
median(dmj_interv)# Median interval
## [1] 52
That’s the end of our regularly scheduled programming. More to come
as I build examples out for dendroTools and other tree-ring
packages. You can email at rmaxwell2@radford.edu and I will attempt to answer
questions.