## The following R code provides everything needed to reproduce results (presented by Pederson et al. [2009]in the journal Climatic Change for single station daily data downloaded from the US Historic Climate Network. I have provided example code specifying file names and working directories based on the file structure and file names on my computer (or collaborators computers).  Do not expect this code to run on your computer without modifying it to reflect your file names and file structure.  This code is free from the authors and comes with ABSOLUTELY NO WARRENTY AND NO GAURENTEE OF SUPPORT. We only ask that if you use this code in support of your analysis you provide the following source citation: 

## Pederson GT, Graumlich LG, Fagre DF, Kipfer TK, and Muhlfeld CC (2009) A century of climate and ecosystem change in western Montana: What do temperature trends portend? Climatic Change XXX:XXX-XXX DOI:10.1007/s10584-009-9642-y)

## Daily HCN Temperature Data Analysis R Code (Run by pasting in R or preferably from a script window)
## Source data from http://cdiac.ornl.gov/epubs/ndp/ushcn/state_MT.html
## Clean up
	rm(list = ls())
## INPUT VARIABLES TO ENTER FOR EACH STATION AND ANALYSIS
	StationName <- c("Bozeman MSU")
	DataInputFile <- c("BMT241044daily_2702.csv")
	StartYearCutOff <- 1893
	# Calculate adjusted year using July 1 as the first day of a "year"
	# This starts a new Julian year with July 1 = 1 and June 31 = 365
		sday <- 1
		smonth <- 7
	# Data and output files in a directory called HCNtempAnalysis
		#Windoze setwd where you have your data file and want outputs
			setwd("C:/Documents and Settings/gpederson/My Documents/Montana Climate/Bozeman MSU/")
		#Mac setwd
			#setwd("/Volumes/Data/Users/todd/Rwork/HCNtempAnalysis")
	# View working directory and files in directory
		ls()
		dir()

## Read HCN data from file
## Assumes data are a .csv file and look like this
	# Source: CN Williams Jr., MJ Menne, RS Vose, DR Easterling, NOAA, National Climatic Data Center, Asheville, NC
	# Day  ,JD  ,Month  ,State_id  ,Year  ,TAVE (F),TMAX (F),TMIN (F) 
	# 1,60,3,'243139',1906,28,45,11
	# 2,61,3,'243139',1906,34,46,21
	# 3,62,3,'243139',1906,26,39,13
	# 4,63,3,'243139',1906,34,41,26
	# 5,64,3,'243139',1906,34,47,21
	# ...
	# Note1a - "." is no data in these files.
	# Note1b - In the last column "." is followed by a space
	# Note1c - Use both "." and ". " for no data fields
	Sta.T <- read.table(DataInputFile, header=T, skip=1, sep=",", na.strings=c(".",". "))
	Sta.T[1:10,]
	Sta.T <- cbind(Sta.T[,5],Sta.T[,3],Sta.T[,1:2], Sta.T[,6:8])
	Sta.T[1:10,]
	names(Sta.T) <- c("Year", "Month", "Day", "JD", "TAVE", "TMAX", "TMIN")
	# Clean-up:  cutoff data for selected start year
	#	Sta.T <- subset(Sta.T, Year >= StartYearCutOff)
	#	row.names(Sta.T) <- seq(1,nrow(Sta.T),by=1)
	# View the data input file
		Sta.T[1:10,]
		summary(Sta.T)

## Create output dataframe that contains sJD and sYear = rescaled year based on start year
	Sta.T.s <- data.frame(matrix(NA, nrow = nrow(Sta.T), ncol = 10))
	names(Sta.T.s) = c("Date", "Month", "Day", "Year", "JD", "sJD", "sYear","TAVE", "TMAX", "TMIN")
	Sta.T.s[1:10,]

## Bring over original data into new dataframe
	Sta.T.s[,2] <- Sta.T$Month
	Sta.T.s[,3] <- Sta.T$Day
	Sta.T.s[,4] <- Sta.T$Year
	Sta.T.s[,5] <- Sta.T$JD
	#Column 6 sJD calculated below
	#Column 7 sYear calculated below
	Sta.T.s[,8] <- Sta.T$TAVE
	Sta.T.s[,9] <- Sta.T$TMAX
	Sta.T.s[,10] <- Sta.T$TMIN
	#Display some rows to visually check data
	Sta.T.s[1:13,]


## Calculate sJD based on sday and smonth set above
## These selections appear more complex than they are.  These change data based on selection sets.
	# Calc sYear for days in calendar year prior to sday and smonth (sYear = previous year)
		Sta.T.s$sYear[Sta.T.s$Month < smonth | (Sta.T.s$Month == smonth & Sta.T.s$Day < sday)] <- Sta.T.s$Year[Sta.T.s$Month < smonth | 		(Sta.T.s$Month == smonth & Sta.T.s$Day < sday)] - 1
	# Calc sYear for days in calendar year equal to or greater than sday and smonth (sYear = current year)
		Sta.T.s$sYear[(Sta.T.s$Month == smonth & Sta.T.s$Day >= sday) | Sta.T.s$Month > smonth] <- Sta.T.s$Year[(Sta.T.s$Month 		== 	smonth & Sta.T.s$Day >= sday) | Sta.T.s$Month > smonth]
	# Create a date field by concatenating year, month, and day.  This helps for calculating new julian days.
		Sta.T.s[,1] <- as.Date(paste(Sta.T.s$Year,"-",Sta.T.s$Month,"-",Sta.T.s$Day,sep=""))
	# Display some rows from the middle to visually check data
		Sta.T.s[150:400,]
	## Calculate the new julian start and end years from the data.  Remember, we are using the new sYear (column 7)
		startyear <- Sta.T.s[1,7]
		endyear <- Sta.T.s[nrow(Sta.T.s),7]
		numyears <- endyear - startyear + 1
	## Calculate the calendar start year and end years for the 90d F data analysis (column 4)
		startyear2 <- Sta.T.s[1,4]
		endyear2 <- Sta.T.s[nrow(Sta.T.s),4]
	# Now we go through the data row by row to calculate the new julian day, sday, for the newly computed julian year, syear 
		for (y in startyear:endyear) {
    		Sta.T.s$sJD[Sta.T.s$sYear == y] <- julian(Sta.T.s$Date[Sta.T.s$sYear == y], origin = as.Date(paste(y,"-",smonth,"-",sday,sep=""))) + 1
		}
	# Display some rows to visually check data
		Sta.T.s[150:400,]
	# Look at summary of the data for sJD to see if this makes sense.
		summary(Sta.T.s$sJD)

## LOOK HERE FOR DEFINITIONS OF OUTPUT DATA FIELDS -- IMPORTANT FOR INTERPRETATION
## Sta.T.summary = Create new output dataframe 
##(sYear, firstmin32, lastmin32, countmin32, firstmin0, lastmin0, countmin0, firstmax32, lastmax32, countmax32, Year90, firstmax90, lastmax90, countmax90)
## 1 sYear = 		new julian year based on smonth and sday (calculated above)
## 2 firstmin32 = 	first day of syear where the MIN temp is <= 32
## 3 lastmin32 = 	last day of syear where the MIN temp is <= 32
## 4 countmin32 = 	number of days in syear where the MIN temp is <= 32
## 5 firstmin0 = 	first day of syear where the MIN temp is <= 0
## 6 lastmin0 = 	last day of syear where the MIN temp is <= 0
## 7 countmin0 = 	number of days in year where the MIN temp is <= 0
## 8 firstmax32 = 	first day of syear where MAX temp is <= 32
## 9 lastmax32 = 	last day of syear where MAX temp is <= 32
## 10 countmax32 = 	number of days in syear where MAX temp is <= 32
## 11 Year90 = 		Calendar year for Max temp >= 90 calculations
## 12 firstmax90 = 	first day of CALENDAR year where the MAX temp is >= 90 -- not constrained by a season (e.g., spring)
## 13 lastmax90 = 	last day of CALENDAR year where the MAX temp is >= 90 -- not constrained by a season (e.g., autumnn)
## 14 countmax90 = number of days in CALENDAR year where the MAX temp is >= 90

## Create the output dataframe
	Sta.T.summary <- data.frame(matrix(NA, nrow = numyears, ncol = 14))
	names(Sta.T.summary) = c("sYear", "firstmin32", "lastmin32", "countmin32", "firstmin0", "lastmin0", "countmin0", "firstmax32","lastmax32","countmax32", "Year90", "firstmax90", "lastmax90", "countmax90")
	Sta.T.summary$sYear <- seq(startyear, endyear, by=1)
	Sta.T.summary$Year90 <- seq(startyear2, endyear2, by=1)
	Sta.T.summary[1:10,]

## CALCULATE OUTPUT DATA

## Revisit some variables we want to use (all set previously)

	startyear # first sYear
	endyear # last sYear
	startyear2 # first calendar year (for >90)
	endyear2 # last calendar year (for >90)

## Loop through data row by row to calulate output data for 32 and 0 thresholds
	for (y in startyear:endyear) {
    	## 32F MIN (sJD in Sta.T.s is column 6)
    		# if then else statement checks that there are actually days TMIN <= 32
    	
    			if (min(Sta.T.s$TMIN[Sta.T.s$sYear == y], na.rm=T) <= 32) {
    				# temp32.year is a temporary holder for TMIN <= 32F days in a year
    					temp32.year <- subset(Sta.T.s, sYear == y & TMIN <= 32)
    				# Create new row names from 1 to number of rows
    					row.names(temp32.year) <- seq(1, nrow(temp32.year), by=1)
    				# First day of syear <= 32F
    					Sta.T.summary$firstmin32[Sta.T.summary$sYear == y] <- temp32.year[1,6]
    				# Last day of syear <= 32F
    					Sta.T.summary$lastmin32[Sta.T.summary$sYear == y] <- temp32.year[nrow(temp32.year),6]
    				# Count number of days <= 32F
    					Sta.T.summary$countmin32[Sta.T.summary$sYear == y] <- nrow(temp32.year)
    			} # close if
    			# if not days <= 32 add NA for start and end and 0 for count
    			else {
    				Sta.T.summary$firstmin32[Sta.T.summary$sYear == y] <- NA
        			Sta.T.summary$lastmin32[Sta.T.summary$sYear == y] <- NA
        			Sta.T.summary$countmin32[Sta.T.summary$sYear == y] <- 0 
    			} #close else	
    	## 0F MIN
    		# if then else statement checks that there are actually days TMIN <= 0
    			if (min(Sta.T.s$TMIN[Sta.T.s$sYear == y], na.rm=T) <= 0) {
        			# temp0.year is a temporary holder for TMIN <= 0F days in a year
        				temp0.year <- subset(Sta.T.s, sYear == y & TMIN <= 0)
        			# Create new row names from 1 to number of rows
        				row.names(temp0.year) <- seq(1, nrow(temp0.year), by=1)
        			# First day of syear <= 0F
        				Sta.T.summary$firstmin0[Sta.T.summary$sYear == y] <- temp0.year[1,6]
        			# Last day of syear <= 0F
        				Sta.T.summary$lastmin0[Sta.T.summary$sYear == y] <- temp0.year[nrow(temp0.year),6]
        			# Count number of days <= 0F
        				Sta.T.summary$countmin0[Sta.T.summary$sYear == y] <- nrow(temp0.year)
    			}  #close if
    			# if not days <= 0 add NA for start and end and 0 for count
    			else {
        			Sta.T.summary$firstmin0[Sta.T.summary$sYear == y] <- NA
        			Sta.T.summary$lastmin0[Sta.T.summary$sYear == y] <- NA
        			Sta.T.summary$countmin0[Sta.T.summary$sYear == y] <- 0 
    			} #close else
    	## 32F MAX
    		# if then else statement checks that there are actually days TMAX <= 32
    			if (min(Sta.T.s$TMAX[Sta.T.s$sYear == y], na.rm=T) <= 32) {
        				temp32max.year <- subset(Sta.T.s, sYear == y & TMAX <= 32)
        				row.names(temp32max.year) <- seq(1, nrow(temp32max.year), by=1)
        				Sta.T.summary$firstmax32[Sta.T.summary$sYear == y] <- temp32max.year[1,6]
        				Sta.T.summary$lastmax32[Sta.T.summary$sYear == y] <- temp32max.year[nrow(temp32max.year),6]
        				Sta.T.summary$countmax32[Sta.T.summary$sYear == y] <- nrow(temp32max.year)
    			} # close if
    			# if not days >= 90 add NA for start and end and 0 for count
    			else {
        			Sta.T.summary$firstmax32[Sta.T.summary$sYear == y] <- NA
        			Sta.T.summary$lastmax32[Sta.T.summary$sYear == y] <- NA
        			Sta.T.summary$countmax32[Sta.T.summary$sYear == y] <- 0 
    			} # close else  	
    } #closes for loop

## Loop through calendar years for calculating max90 data
    for (y in startyear2:endyear2) {
    	## 90F MAX
    			if (max(Sta.T.s$TMAX[Sta.T.s$Year == y], na.rm=T) >= 90) {
        				temp90max.year <- subset(Sta.T.s, Year == y & TMAX >= 90)
        				row.names(temp90max.year) <- seq(1, nrow(temp90max.year), by=1)
        				Sta.T.summary$firstmax90[Sta.T.summary$sYear == y] <- temp90max.year[1,5]
        				Sta.T.summary$lastmax90[Sta.T.summary$sYear == y] <- temp90max.year[nrow(temp90max.year),5]
        				Sta.T.summary$countmax90[Sta.T.summary$sYear == y] <- nrow(temp90max.year)
    			} # close if
    			# if not days >= 90 add NA for start and end and 0 for count
    			else {
        			Sta.T.summary$firstmax90[Sta.T.summary$sYear == y] <- NA
        			Sta.T.summary$lastmax90[Sta.T.summary$sYear == y] <- NA
        			Sta.T.summary$countmax90[Sta.T.summary$sYear == y] <- 0 
    			} # close lese  			
	} #closes for loop

Sta.T.summary
Sta.T.summary <- subset(Sta.T.summary, sYear >= StartYearCutOff & sYear <= 2006)
Sta.T.summary
summary(Sta.T.summary)
write.table(Sta.T.summary, file=paste(StationName,"-","HCNday-analysis.csv", sep=""), quote=F, sep=",", row.names=F, col.names=T) 

##Linear Models of significance testing of trends in first, last, and total count of days
	#Min32F
	F32minstart.lm <- lm(Sta.T.summary$firstmin32[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006])
	summary(F32minstart.lm)
	F32minend.lm <- lm(Sta.T.summary$lastmin32[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006])
	summary(F32minend.lm)
	F32mincount.lm <- lm(Sta.T.summary$countmin32[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006])
	summary(F32mincount.lm)

	#Min0F
	F0start.lm <- lm(Sta.T.summary$firstmin0[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006])
	summary(F0start.lm)
	F0end.lm <- lm(Sta.T.summary$lastmin0[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006])
	summary(F0end.lm)
	F0count.lm <- lm(Sta.T.summary$countmin0[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006])
	summary(F0count.lm)
	
	#Max32F
	F32maxstart.lm <- lm(Sta.T.summary$firstmax32[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006])
	summary(F32maxstart.lm)
	F32maxend.lm <- lm(Sta.T.summary$lastmax32[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006])
	summary(F32maxend.lm)
	F32maxcount.lm <- lm(Sta.T.summary$countmax32[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006])
	summary(F32maxcount.lm)

	#90F
	F90start.lm <- lm(Sta.T.summary$firstmax90[Sta.T.summary$Year90 > StartYearCutOff & Sta.T.summary$Year90 < 2006] ~ Sta.T.summary$Year90[Sta.T.summary$Year90 > StartYearCutOff & Sta.T.summary$Year90 < 2006])
	summary(F90start.lm)
	F90end.lm <- lm(Sta.T.summary$lastmax90[Sta.T.summary$Year90 > StartYearCutOff & Sta.T.summary$Year90 < 2006] ~ Sta.T.summary$sYear[Sta.T.summary$Year90 > StartYearCutOff & Sta.T.summary$Year90 < 2006])
	summary(F90end.lm)
	F90count.lm <- lm(Sta.T.summary$countmax90[Sta.T.summary$Year90 > StartYearCutOff & Sta.T.summary$Year90 < 2006] ~ Sta.T.summary$Year90[Sta.T.summary$Year90 > StartYearCutOff & Sta.T.summary$Year90 < 2006])
	summary(F90count.lm)
	
##############################################
###                                        ###
##                 Plot fest!               ##
###                                        ###
##############################################

pdf("Bozeman MSU_RegressionSummary.pdf",8,8)
par(mar = rep(3,4)+0.1)

#F32min Regression summary plots
opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(F32minstart.lm, las = 1, main=StationName)
plot(F32minend.lm, las=1, main=StationName)
plot(F32mincount.lm, las=1, main=StationName)
par(opar)
#F0 Regression summary plots
opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(F0start.lm, las = 1, main=StationName)
plot(F0end.lm, las = 1, main=StationName)
plot(F0count.lm, las = 1, main=StationName)
par(opar)
#F32max Regression summary plots
opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(F32maxstart.lm, las = 1, main=StationName)
plot(F32maxend.lm, las=1, main=StationName)
plot(F32maxcount.lm, las=1, main=StationName)
par(opar)
#F90 Regression summary plots
opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(F90start.lm, las = 1, main=StationName)
plot(F0end.lm, las = 1, main=StationName)
plot(F90count.lm, las = 1, main=StationName)
par(opar)

pdf("Bozeman MSU_DaySummary.pdf",8,4)
par(mar = rep(3,4)+0.1)

#########                          ###########
##                                          ##
### Plot 1. First/Last Frost Days          ###
##                                          ##
#########                          ###########

plot(Sta.T.summary$sYear, Sta.T.summary$firstmin32, type="l", col="black", ylim=c(1,365), xlim=c(StartYearCutOff,2006), xlab="Year", ylab="Winter Year Day", main="Bozeman MSU (HCN data) Start and End Days TMIN <= 32F")
lines(Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006], Sta.T.summary$lastmin32[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006], col="black")
lines(lowess(Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006],Sta.T.summary$firstmin32[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006], f=5/100), col="blue", lwd=1.5)
lines(lowess(Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006],Sta.T.summary$lastmin32[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006], f=5/100), col="green", lwd=1.5)
abline(lsfit(Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006],Sta.T.summary$firstmin32[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006]), col="blue", lwd=1)
abline(lsfit(Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006],Sta.T.summary$lastmin32[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006]), col="green", lwd=1)

#########                          ###########
##                                          ##
### Plot 2. First/Last < 0 Cold Days       ###
##                                          ##
#########                          ###########

plot(Sta.T.summary$sYear, Sta.T.summary$firstmin0, type="l", col="black", lwd=0.5, ylim=c(120,300), xlim=c(StartYearCutOff,2006), xlab="Year", ylab="Winter Year Day", main="Bozeman MSU (HCN data) Start and End Days TMIN <= 0F")
lines(Sta.T.summary$sYear, Sta.T.summary$lastmin0, col="black")
abline(lsfit(Sta.T.summary$sYear, Sta.T.summary$firstmin0), col="blue")
abline(lsfit(Sta.T.summary$sYear, Sta.T.summary$lastmin0), col="green")

#########                          ###########
##                                          ##
### Plot 3. First/Last max < 32            ###
##                                          ##
#########                          ###########

plot(Sta.T.summary$sYear, Sta.T.summary$firstmax32, type="l", col="black", ylim=c(1,365), xlim=c(StartYearCutOff,2006), xlab="Year", ylab="Winter Year Day", main="Bozeman MSU (HCN data) Start and End Days TMAX <= 32F")
lines(Sta.T.summary$sYear, Sta.T.summary$lastmax32, col="black")
abline(lsfit(Sta.T.summary$sYear, Sta.T.summary$firstmax32), col="green")
abline(lsfit(Sta.T.summary$sYear, Sta.T.summary$lastmax32), col="blue")

#########                          ###########
##                                          ##
### Plot 4. Number >90 Frickin' Hot Days   ###
##                                          ##
#########                          ###########

plot(Sta.T.summary$Year90, Sta.T.summary$firstmax90, type="l", col="black", ylim=c(1,365), xlim=c(StartYearCutOff,2006), xlab="Year", ylab="Calendar Year Day", main="Bozeman MSU (HCN data) Start and End Days TMAX >= 90F")
lines(Sta.T.summary$Year90, Sta.T.summary$lastmax90, col="black")
abline(lsfit(Sta.T.summary$sYear, Sta.T.summary$firstmax90), col="green")
abline(lsfit(Sta.T.summary$sYear, Sta.T.summary$lastmax90), col="blue")

#########                          ###########
##                                          ##
### Plot 5. Day Counts for temp. thresholds###
##                                          ##
#########                          ###########

barplot(Sta.T.summary$countmin32,names.arg=Sta.T.summary$sYear, width=1, space=0.75, xlab="Year", ylab="Number of days", main="Bozeman MSU (HCN data) Days Tmin <= 32F")
barplot(Sta.T.summary$countmin0,names.arg=Sta.T.summary$sYear, width=1, space=0.75, xlab="Year", ylab="Number of days", main="Bozeman MSU (HCN data) Days Tmin <= 0F")
barplot(Sta.T.summary$countmax32,names.arg=Sta.T.summary$sYear, width=1, space=0.75, xlab="Year", ylab="Number of days", main="Bozeman MSU (HCN data) Days Tmax <= 32F")
barplot(Sta.T.summary$countmax90,names.arg=Sta.T.summary$Year90, width=1, space=0.75, xlab="Year", ylab="Number of days", main="Bozeman MSU (HCN data) Days Tmax >= 90F")