## 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.  This version contains modifications to the original script made by Ali Macalady at the University of Arizona.  Modification here include new absolute temperature thresholds relevant to the southwestern U.S. and an additional analysis that track yearly extreme high and low temperatures. 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)
## 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("Gage 4ESE, NM")
	DataInputFile <- c("NM293368_Gage4ESE.csv")
	StartYearCutOff <- 1948
	# 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("/Users/alisonmacalady/Desktop/ClassesSpring2008/ClimateandWater_GEOG552/project/Daily_data")
	# 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',1956,28,45,11
	# 2,61,3,'243139',1956,34,46,21
	# 3,62,3,'243139',1956,26,39,13
	# 4,63,3,'243139',1956,34,41,26
	# 5,64,3,'243139',1956,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 95d 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, firstmin15, lastmin15, countmin15, firstmax32, lastmax32 countmax32, yearlymin,Year95, firstmax95, lastmax95, countmax95, yearlymax)
## 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 firstmin15 = 	first day of syear where the MIN temp is <= 15
## 6 lastmin15 = 	last day of syear where the MIN temp is <= 15
## 7 countmin15 = 	number of days in year where the MIN temp is <= 15
## 8 firstmax32 = 	MAX temp in the year
## 9 lastmax32 = 	Min temp in sYear
## 10 countmax32 = 	number of days in syear where MAX temp is <= 32
## 11 yearlymin = Min temp in sYear
## 12 Year95 = 		Calendar year for Max temp >= 95 calculations
## 13 firstmax95 = 	first day of CALENDAR year where the MAX temp is >= 95 -- not constrained by a season (e.g., spring)
## 14 lastmax95 = 	last day of CALENDAR year where the MAX temp is >= 95 -- not constrained by a season (e.g., autumnn)
## 15 countmax95 = number of days in CALENDAR year where the MAX temp is >= 95
## 16 firstmax100 = 	first day of CALENDAR year where the MAX temp is >= 100 -- not constrained by a season (e.g., spring)
## 17 lastmax100 = 	last day of CALENDAR year where the MAX temp is >= 100 -- not constrained by a season (e.g., autumnn)
## 18 countmax100 = number of days in CALENDAR year where the MAX temp is >= 100
## 19 yearlymax = Max temp in year

## Create the output dataframe
	Sta.T.summary <- data.frame(matrix(NA, nrow = (numyears), ncol = 19))
	names(Sta.T.summary) = c("sYear", "firstmin32", "lastmin32", "countmin32", "firstmin15", "lastmin15", "countmin15", "firstmax32","lastmax32","countmax32","yearlymin", "Year95", "firstmax95", "lastmax95", "countmax95", "firstmax100", "lastmax100", "countmax100", "yearlymax")
	Sta.T.summary$sYear <- seq(startyear, (endyear), by=1)
	Sta.T.summary$Year95 <- seq((startyear2-1), (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 >95)
	endyear2 # last calendar year (for >95)

## Loop through data row by row to calulate output data for 32 and 15 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	
    	## 15F MIN
    		# if then else statement checks that there are actually days TMIN <= 15
    			if (min(Sta.T.s$TMIN[Sta.T.s$sYear == y], na.rm=T) <= 15) {
        			# temp15.year is a temporary holder for TMIN <= 15F days in a year
        				temp15.year <- subset(Sta.T.s, sYear == y & TMIN <= 15)
        			# Create new row names from 1 to number of rows
        				row.names(temp15.year) <- seq(1, nrow(temp15.year), by=1)
        			# First day of syear <= 15F
        				Sta.T.summary$firstmin15[Sta.T.summary$sYear == y] <- temp15.year[1,6]
        			# Last day of syear <= 15F
        				Sta.T.summary$lastmin15[Sta.T.summary$sYear == y] <- temp15.year[nrow(temp15.year),6]
        			# Count number of days <= 15F
        				Sta.T.summary$countmin15[Sta.T.summary$sYear == y] <- nrow(temp15.year)
    			}  #close if
    			# if not days <= 15 add NA for start and end and 0 for count
    			else {
        			Sta.T.summary$firstmin15[Sta.T.summary$sYear == y] <- NA
        			Sta.T.summary$lastmin15[Sta.T.summary$sYear == y] <- NA
        			Sta.T.summary$countmin15[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 >= 32 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
##Calculate yearly min temp
 for (y in startyear:endyear) {
 Sta.T.summary$yearlymin[Sta.T.summary$sYear == y] <- min(Sta.T$TMIN[Sta.T.s$sYear == y], na.rm=T)
    }
    
## Loop through calendar years for calculating max95 data
    for (y in startyear2:endyear2) {
    	## 95F MAX
    			if (max(Sta.T.s$TMAX[Sta.T.s$Year == y], na.rm=T) >= 95) {
        				temp95max.year <- subset(Sta.T.s, Year == y & TMAX >= 95)
        				row.names(temp95max.year) <- seq(1, nrow(temp95max.year), by=1)
        				Sta.T.summary$firstmax95[Sta.T.summary$Year95 == y] <- temp95max.year[1,5]
        				Sta.T.summary$lastmax95[Sta.T.summary$Year95 == y] <- temp95max.year[nrow(temp95max.year),5]
        				Sta.T.summary$countmax95[Sta.T.summary$Year95 == y] <- nrow(temp95max.year)
    			} # close if
    			# if not days >= 95 add NA for start and end and 0 for count
    			else {
        			Sta.T.summary$firstmax95[Sta.T.summary$Year95 == y] <- NA
        			Sta.T.summary$lastmax95[Sta.T.summary$Year95 == y] <- NA
        			Sta.T.summary$countmax95[Sta.T.summary$Year95 == y] <- 0 
    			} # close else  			
	} #closes for loop
## Loop through calendar years for calculating max95 data
    for (y in startyear2:endyear2) {
    	## 100F MAX
    			if (max(Sta.T.s$TMAX[Sta.T.s$Year == y], na.rm=T) >= 100) {
        				temp100max.year <- subset(Sta.T.s, Year == y & TMAX >= 100)
        				row.names(temp100max.year) <- seq(1, nrow(temp100max.year), by=1)
        				Sta.T.summary$firstmax100[Sta.T.summary$Year95 == y] <- temp100max.year[1,5]
        				Sta.T.summary$lastmax100[Sta.T.summary$Year95 == y] <- temp100max.year[nrow(temp100max.year),5]
        				Sta.T.summary$countmax100[Sta.T.summary$Year95 == y] <- nrow(temp100max.year)
    			} # close if
    			# if not days >= 100 add NA for start and end and 0 for count
    			else {
        			Sta.T.summary$firstmax100[Sta.T.summary$Year95 == y] <- NA
        			Sta.T.summary$lastmax100[Sta.T.summary$Year95 == y] <- NA
        			Sta.T.summary$countmax100[Sta.T.summary$Year95 == y] <- 0 
    			} # close else  			
	} #closes for loop
	
##Calculate yearly max temp
 for (y in startyear2:endyear2) {
 Sta.T.summary$yearlymax[Sta.T.summary$Year95 == y] <- max(Sta.T$TMAX[Sta.T.s$Year == y], na.rm=T)
    }
Sta.T.summary
Sta.T.summary <- subset(Sta.T.summary, sYear >= (StartYearCutOff+1) & Year95 <= 2004)
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 < 2005] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005])
	summary(F32minstart.lm)
	F32minend.lm <- lm(Sta.T.summary$lastmin32[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005])
	summary(F32minend.lm)
	F32mincount.lm <- lm(Sta.T.summary$countmin32[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005])
	summary(F32mincount.lm)

	#Min15F
	F15start.lm <- lm(Sta.T.summary$firstmin15[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005])
	summary(F15start.lm)
	F15end.lm <- lm(Sta.T.summary$lastmin15[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005])
	summary(F15end.lm)
	F15count.lm <- lm(Sta.T.summary$countmin15[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005])
	summary(F15count.lm)
	
	#Max32F
	F32maxstart.lm <- lm(Sta.T.summary$firstmax32[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005])
	summary(F32maxstart.lm)
	F32maxend.lm <- lm(Sta.T.summary$lastmax32[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005])
	summary(F32maxend.lm)
	F32maxcount.lm <- lm(Sta.T.summary$countmax32[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005] ~ Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2005])
	summary(F32maxcount.lm)

	#95F
	F95start.lm <- lm(Sta.T.summary$firstmax95[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2005] ~ Sta.T.summary$Year95[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2005])
	summary(F95start.lm)
	F95end.lm <- lm(Sta.T.summary$lastmax95[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2005] ~ Sta.T.summary$sYear[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2005])
	summary(F95end.lm)
	F95count.lm <- lm(Sta.T.summary$countmax95[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2005] ~ Sta.T.summary$Year95[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2005])
	summary(F95count.lm)
	
		#100F
	F100start.lm <- lm(Sta.T.summary$firstmax100[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2005] ~ Sta.T.summary$Year95[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2005])
	summary(F100start.lm)
	F100end.lm <- lm(Sta.T.summary$lastmax100[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2005] ~ Sta.T.summary$sYear[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2005])
	summary(F100end.lm)
	F100count.lm <- lm(Sta.T.summary$countmax100[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2005] ~ Sta.T.summary$Year95[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2005])
	summary(F100count.lm)
	
			#min&max temps
	Fyearlymin.lm <- lm(Sta.T.summary$yearlymin[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2005] ~ Sta.T.summary$Year95[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2005])
	summary(Fyearlymin.lm)
	
	Fyearlymax.lm <- lm(Sta.T.summary$yearlymax[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2005] ~ Sta.T.summary$sYear[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2005])
	summary(Fyearlymax.lm)
	
##############################################
###                                        ###
##                 Plot fest!               ##
###                                        ###
##############################################

pdf("Gage 4ESE, NM _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)
#F15 Regression summary plots
opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(F15start.lm, las = 1, main=StationName)
plot(F15end.lm, las = 1, main=StationName)
plot(F15count.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)
#F95 Regression summary plots
opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(F95start.lm, las = 1, main=StationName)
plot(F95end.lm, las = 1, main=StationName)
plot(F95count.lm, las = 1, main=StationName)
par(opar)
#F100 Regression summary plots
opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(F100start.lm, las = 1, main=StationName)
plot(F100end.lm, las = 1, main=StationName)
plot(F100count.lm, las = 1, main=StationName)
par(opar)
pdf("Gage 4ESE, NM _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="Gage 4ESE, NM  (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 < 15 Cold Days       ###
##                                          ##
#########                          ###########

plot(Sta.T.summary$sYear, Sta.T.summary$firstmin15, type="l", col="black", lwd=0.5, ylim=c(120,300), xlim=c(StartYearCutOff,2006), xlab="Year", ylab="Winter Year Day", main="Gage 4ESE, NM  (HCN data) Start and End Days TMIN <= 15F")
lines(Sta.T.summary$sYear, Sta.T.summary$lastmin15, col="black")
abline(lsfit(Sta.T.summary$sYear, Sta.T.summary$firstmin15), col="blue")
abline(lsfit(Sta.T.summary$sYear, Sta.T.summary$lastmin15), 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="Gage 4ESE, NM  (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 >95 Frickin' Hot Days   ###
##                                          ##
#########                          ###########

plot(Sta.T.summary$Year95, Sta.T.summary$firstmax95, type="l", col="black", ylim=c(1,365), xlim=c(StartYearCutOff,2006), xlab="Year", ylab="Calendar Year Day", main="Gage 4ESE, NM  (HCN data) Start and End Days TMAX >= 95F")
lines(Sta.T.summary$Year95, Sta.T.summary$lastmax95, col="black")
abline(lsfit(Sta.T.summary$sYear, Sta.T.summary$firstmax95), col="green")
abline(lsfit(Sta.T.summary$sYear, Sta.T.summary$lastmax95), col="blue")
#########                          ###########
##                                          ##
### Plot 5. Number >100 Frickin' Hot Days   ###
##                                          ##
#########                          ###########

plot(Sta.T.summary$Year95, Sta.T.summary$firstmax100, type="l", col="black", ylim=c(1,365), xlim=c(StartYearCutOff,2006), xlab="Year", ylab="Calendar Year Day", main="Gage 4ESE, NM  (HCN data) Start and End Days TMAX >= 100F")
lines(Sta.T.summary$Year95, Sta.T.summary$lastmax100, col="black")
abline(lsfit(Sta.T.summary$sYear, Sta.T.summary$firstmax100), col="green")
abline(lsfit(Sta.T.summary$sYear, Sta.T.summary$lastmax100), col="blue")

#########                          ###########
##                                          ##
### Plot 6/7. Yearly Min/Max temperatures  ###
##                                          ##
#########                          ###########

plot(Sta.T.summary$sYear, Sta.T.summary$yearlymin, type="l", col="black", ylim=c(-20,50), xlim=c(StartYearCutOff,2006), xlab="Winter Year", ylab="Degrees F", main="Gage 4ESE, NM  (HCN data) Yearly Min Temps")
lines(Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006], Sta.T.summary$yearlymin[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$yearlymin[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006], f=5/100), col="blue", lwd=1.5)
abline(lsfit(Sta.T.summary$sYear[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006],Sta.T.summary$yearlymin[Sta.T.summary$sYear > StartYearCutOff & Sta.T.summary$sYear < 2006]), col="blue", lwd=1)


plot(Sta.T.summary$Year95, Sta.T.summary$yearlymax, type="l", col="black", ylim=c(85,120), xlim=c(StartYearCutOff,2006), xlab="Year", ylab="Degrees F", main="Gage 4ESE, NM  (HCN data) Yearly Max Temps")
lines(Sta.T.summary$sYear[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$sYear < 2006], Sta.T.summary$yearlymin[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2006], col="black")
lines(lowess(Sta.T.summary$Year95[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2006],Sta.T.summary$yearlymax[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2006], f=5/100), col="blue", lwd=1.5)
abline(lsfit(Sta.T.summary$Year95[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2006],Sta.T.summary$yearlymax[Sta.T.summary$Year95 > StartYearCutOff & Sta.T.summary$Year95 < 2006]), col="blue", lwd=1)




#########                          ###########
##                                          ##
### Plot 8. 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="Gage 4ESE, NM  (HCN data) Days Tmin <= 32F")
barplot(Sta.T.summary$countmin15,names.arg=Sta.T.summary$sYear, width=1, space=0.75, xlab="Year", ylab="Number of days", main="Gage 4ESE, NM  (HCN data) Days Tmin <= 15F")
barplot(Sta.T.summary$countmax32,names.arg=Sta.T.summary$sYear, width=1, space=0.75, xlab="Year", ylab="Number of days", main="Gage 4ESE, NM  (HCN data) Days Tmax <= 32F")
barplot(Sta.T.summary$countmax95,names.arg=Sta.T.summary$Year95, width=1, space=0.75, xlab="Year", ylab="Number of days", main="Gage 4ESE, NM  (HCN data) Days Tmax >= 95F")
barplot(Sta.T.summary$countmax100,names.arg=Sta.T.summary$Year95, width=1, space=0.75, xlab="Year", ylab="Number of days", main="Gage 4ESE, NM  (HCN data) Days Tmax >= 100F")