Objectives
- Create groups of trees by temperature and growth pattern for Mount Washington site
- Compare to results of Salzer et al. (2014) from White Mountains
library(dplR)
## Warning: package 'dplR' was built under R version 3.1.3
library(treeclim)
## Warning: package 'treeclim' was built under R version 3.1.3
## Loading required package: Rcpp
## Warning: package 'Rcpp' was built under R version 3.1.3
## This is treeclim version 1.0.11.
##
## Need help? See the wiki at https://github.com/cszang/treeclim/wiki
## (Use suppressPackageStartupMessages to eliminate package startup messages.)
library(rgdal)
## Warning: package 'rgdal' was built under R version 3.1.3
## Loading required package: sp
## Warning: package 'sp' was built under R version 3.1.3
## rgdal: version: 0.9-2, (SVN revision 526)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
## Path to GDAL shared files: C:/Users/Tyler/Documents/R/win-library/3.1/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: C:/Users/Tyler/Documents/R/win-library/3.1/rgdal/proj
library(lattice)
library(Hmisc)
## Warning: package 'Hmisc' was built under R version 3.1.3
## Loading required package: grid
## Loading required package: survival
## Loading required package: splines
## Loading required package: Formula
## Warning: package 'Formula' was built under R version 3.1.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.1.3
##
## Attaching package: 'Hmisc'
##
## The following objects are masked from 'package:base':
##
## format.pval, round.POSIXt, trunc.POSIXt, units
Reading in tree-ring data and calculating statistics
mwa <- read.tucson(fname = "MWA_308series.rwl",long=T)
mwa <- mwa[,-(1:13)]
mwa <- mwa[,-6]
mwa <- mwa[,-9]
mwa <- mwa[,-(12:13)]
mwa.stats <- summary(mwa)
Reading in temperature data for each tree
mwaTrees <- readOGR(".","mwaTreesTopoTmin")
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "mwaTreesTopoTmin"
## with 291 features
## It has 34 fields
trees2get <- colnames(mwa)
mwaTrees <- mwaTrees[mwaTrees$MSMT_ID2 %in% trees2get,]
mwaTmins <- data.frame(sample=mwaTrees@data$MSMT_ID2, elev=mwaTrees@data$elev, tminMay=mwaTrees@data$tminMay, tminJun=mwaTrees@data$tminJun,
tminJul=mwaTrees@data$tminJul, tminAug=mwaTrees@data$tminAug, tminSep=mwaTrees@data$tminSep)
# Write as a csv and then rearrange in descending order of Tmin in excel
mwaTmins <- read.csv('mwaTmins.csv')
mwaTmins.elev <- read.csv('mwaTmins.elev.csv')
rownames(mwaTmins) <- mwaTmins$sample
rownames(mwaTmins.elev) <- mwaTmins$sample
acf plot for mean chronology
mwaRWI <- detrend(mwa,method="Mean")
mwa.crn <- chron(mwaRWI)
acf(mwa.crn$xxxstd)
pacf(mwa.crn$xxxstd)
Plotting AR1 for all trees vs elevation and temperature
plot(mwaTmins$elev, mwa.stats$ar1, xlab='Elevation (m)', ylab='AR(1)')
plot(mwaTmins$tminGS, mwa.stats$ar1, xlab='Growing season minimum temperature (C)', ylab='AR(1)')
plot(mwaTmins$elev, mwaTmins$tminGS, xlab='Elevation (m)', ylab='Growing season minimum temperature (C)')
Create 10 elevation bins
# By elevation
bin10 <- mwaTmins.elev[1:13,] # bin 10 is warmest, bin 1 is coldest
bin9 <- mwaTmins.elev[14:26,]
bin8 <- mwaTmins.elev[27:39,]
bin7 <- mwaTmins.elev[40:52,]
bin6 <- mwaTmins.elev[53:65,]
bin5 <- mwaTmins.elev[66:78,]
bin4 <- mwaTmins.elev[79:91,]
bin3 <- mwaTmins.elev[92:104,]
bin2 <- mwaTmins.elev[105:117,]
bin1 <- mwaTmins.elev[118:129,]
rownames(mwaTrees@data) <- mwaTrees@data$MSMT_ID2
rownames(mwa.stats) <- mwa.stats$series
bin10.stats <- mwa.stats[match(rownames(bin10), rownames(mwa.stats)),]
bin9.stats <- mwa.stats[match(rownames(bin9), rownames(mwa.stats)),]
bin8.stats <- mwa.stats[match(rownames(bin8), rownames(mwa.stats)),]
bin7.stats <- mwa.stats[match(rownames(bin7), rownames(mwa.stats)),]
bin6.stats <- mwa.stats[match(rownames(bin6), rownames(mwa.stats)),]
bin5.stats <- mwa.stats[match(rownames(bin5), rownames(mwa.stats)),]
bin4.stats <- mwa.stats[match(rownames(bin4), rownames(mwa.stats)),]
bin3.stats <- mwa.stats[match(rownames(bin3), rownames(mwa.stats)),]
bin2.stats <- mwa.stats[match(rownames(bin2), rownames(mwa.stats)),]
bin1.stats <- mwa.stats[match(rownames(bin1), rownames(mwa.stats)),]
se <- function(x) sd(x)/sqrt(length(x))
bin.clust.elev <- data.frame(xbar = c(mean(bin10.stats$ar1), mean(bin9.stats$ar1), mean(bin8.stats$ar1),
mean(bin7.stats$ar1), mean(bin6.stats$ar1), mean(bin5.stats$ar1),
mean(bin4.stats$ar1), mean(bin3.stats$ar1), mean(bin2.stats$ar1),
mean(bin1.stats$ar1)), se = c(se(bin10.stats$ar1), se(bin9.stats$ar1),
se(bin8.stats$ar1), se(bin7.stats$ar1), se(bin6.stats$ar1),
se(bin5.stats$ar1), se(bin4.stats$ar1), se(bin3.stats$ar1),
se(bin2.stats$ar1), se(bin1.stats$ar1)), cluster = c('3505.8 (bin10)',
'3489.9 (bin9)','3481.2 (bin8)','3459.4 (bin7)','3434.6 (bin6)','3417.6 (bin5)',
'3413.2 (bin4)','3406.1 (bin3)','3401.4 (bin2)','3398.8 (bin1)'))
Dotplot(cluster~Cbind(xbar, xbar-se, xbar+se), bin.clust.elev,
xlab="AR(1)", ylab = "Bin mean elevation (m)",las =2)
Create 10 temperature bins
bin10 <- mwaTmins[1:13,] # bin 10 is warmest, bin 1 is coldest
bin9 <- mwaTmins[14:26,]
bin8 <- mwaTmins[27:39,]
bin7 <- mwaTmins[40:52,]
bin6 <- mwaTmins[53:65,]
bin5 <- mwaTmins[66:78,]
bin4 <- mwaTmins[79:91,]
bin3 <- mwaTmins[92:104,]
bin2 <- mwaTmins[105:117,]
bin1 <- mwaTmins[118:129,]
bin10.stats <- mwa.stats[match(rownames(bin10), rownames(mwa.stats)),]
bin9.stats <- mwa.stats[match(rownames(bin9), rownames(mwa.stats)),]
bin8.stats <- mwa.stats[match(rownames(bin8), rownames(mwa.stats)),]
bin7.stats <- mwa.stats[match(rownames(bin7), rownames(mwa.stats)),]
bin6.stats <- mwa.stats[match(rownames(bin6), rownames(mwa.stats)),]
bin5.stats <- mwa.stats[match(rownames(bin5), rownames(mwa.stats)),]
bin4.stats <- mwa.stats[match(rownames(bin4), rownames(mwa.stats)),]
bin3.stats <- mwa.stats[match(rownames(bin3), rownames(mwa.stats)),]
bin2.stats <- mwa.stats[match(rownames(bin2), rownames(mwa.stats)),]
bin1.stats <- mwa.stats[match(rownames(bin1), rownames(mwa.stats)),]
bin.clust <- data.frame(xbar = c(mean(bin10.stats$ar1), mean(bin9.stats$ar1), mean(bin8.stats$ar1),
mean(bin7.stats$ar1), mean(bin6.stats$ar1), mean(bin5.stats$ar1),
mean(bin4.stats$ar1), mean(bin3.stats$ar1), mean(bin2.stats$ar1),
mean(bin1.stats$ar1)), se = c(se(bin10.stats$ar1), se(bin9.stats$ar1),
se(bin8.stats$ar1), se(bin7.stats$ar1), se(bin6.stats$ar1),
se(bin5.stats$ar1), se(bin4.stats$ar1), se(bin3.stats$ar1),
se(bin2.stats$ar1), se(bin1.stats$ar1)), cluster = c('6.182 (bin10)',
'6.017 (bin9)','5.885 (bin8)','5.732 (bin7)','5.635 (bin6)','5.570 (bin5)','5.501 (bin4)','5.237 (bin3)','5.168 (bin2)','5.142 (bin1)'))
Dotplot(cluster~Cbind(xbar, xbar-se, xbar+se), bin.clust,
xlab="AR(1)", ylab = "Bin mean growing season min temp (C)",las =2)
Looking at ar1 by growth clusters
load("C:/Users/Tyler/Dropbox/treelineTMCs (1)/working/tyler/522analysis/mwaMonthlyPRISM.Rdata")
mwaPRISM <- mwaPRISM[,c(1:3,6)]
# Separate by clusters
yrs <- as.numeric(rownames(mwa))
cut.yr <- 1895 # start at 1895
rows2keep <- yrs >= cut.yr
mwa <- mwa[rows2keep,]
yrs <- as.numeric(rownames(mwa))
cut.yr <- 2002 # keep only series that go until 2002
rows2keep <- yrs <= cut.yr
mwa <- mwa[rows2keep,]
# get rid of series that don't fit completely in this time frame
cols2keep <- apply(mwa,2,function(x){ !any(is.na(x))})
mwa <- mwa[,cols2keep]
mwaRWI <- detrend(mwa,method="Mean")
d <- dist(t(mwaRWI), method = "euclidean") # distance matrix
fit <- hclust(d, method="ward.D")
plot(fit, labels=F, main='Growth Cluster Dendrogram', xlab='') # display dendogram
# four clusters?
ngrp <- 4
clusters.growth <- rect.hclust(fit, k=ngrp, border=c('cornflowerblue','darkolivegreen','mediumpurple3','mediumpurple3'))
clusts <- cutree(fit, k=ngrp)
cluster1.growth <- names(clusters.growth[[1]])
cluster2.growth <- names(clusters.growth[[2]])
cluster3.growth <- names(clusters.growth[[3]])
cluster4.growth <- names(clusters.growth[[4]])
mwa1 <- mwa.stats[match(cluster1.growth, rownames(mwa.stats)),]
mwa2 <- mwa.stats[match(cluster2.growth, rownames(mwa.stats)),]
mwa3 <- mwa.stats[match(cluster3.growth, rownames(mwa.stats)),]
mwa4 <- mwa.stats[match(cluster4.growth, rownames(mwa.stats)),]
# Clean up, get rid of NAs
mwa1 <- mwa1[-(1:4),]
mwa3 <- mwa3[-1,]
mwa4 <- mwa4[-1,]
bin.4clust <- data.frame(xbar = c(mean(mwa1$ar1), mean(mwa2$ar1), mean(mwa3$ar1),
mean(mwa4$ar1)), se = c(0.01887279, 0.02507388,
0.02579756, 0.02787271), cluster = c('cluster 1',
'cluster 2', 'cluster 3', 'cluster 4'))
bin.4clust$se <- as.numeric(bin.4clust$se)
Fourclusters <- c(mean(mwa1$ar1), mean(mwa2$ar1), mean(mwa3$ar1),
mean(mwa4$ar1))
se4clust <- c(se(mwa1$ar1), se(mwa2$ar1),
se(mwa3$ar1), se(mwa4$ar1))
dotchart(Fourclusters, labels = c('Cluster 1', 'Cluster 2', 'Cluster 3', 'Cluster 4'),
xlim = c(0.5, 0.73), xlab = 'AR(1) coefficient')
segments(Fourclusters-se4clust, 1:4, Fourclusters+se4clust, 1:4)
