## 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("Fortine")
    DataInputFile <- c("MT243139_2605.csv")
    StartYearCutOff <- 1907
    # 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("G:/Rwork/HCNtempAnalysis")
        #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,]
        dim(Sta.T.s)
    # Look at summary of the data for sJD to see if this makes sense.
        summary(Sta.T.s$sJD)

## Output data for each analysis category in matrices of year and day

Sta.T.TMAX <- data.frame(matrix(0, nrow = numyears, ncol = 365))
names(Sta.T.TMAX) <- c(seq(1,365, by=1), "sYear")
Sta.T.TMAX$sYear <- seq(startyear, endyear, by=1)

Sta.T.TMIN <- data.frame(matrix(0, nrow = numyears, ncol = 365))
names(Sta.T.TMIN) <- c(seq(1,365, by=1), "sYear")
Sta.T.TMIN$sYear <- seq(startyear, endyear, by=1)


    for (y in 1:numyears) {
        for(d in 1:365) {
            if ( length(Sta.T.s$JD[Sta.T.s$Year == (y + startyear -1) & Sta.T.s$JD == d]) > 0) { 
                Sta.T.TMAX[y,d] <- Sta.T.s$TMAX[Sta.T.s$Year == (y + startyear -1) & Sta.T.s$JD == d]
            }
            else {
                Sta.T.TMAX[y,d] <- NA
                
            }
        }
    }
    

    for (y in 1:numyears) {
        for(d in 1:365) {
            if ( length(Sta.T.s$sJD[Sta.T.s$sYear == (y + startyear -1) & Sta.T.s$sJD == d]) > 0) { 
                Sta.T.TMIN[y,d] <- Sta.T.s$TMIN[Sta.T.s$sYear == (y + startyear -1) & Sta.T.s$sJD == d]
            }
            else {
                Sta.T.TMIN[y,d] <- NA
                
            }
        }
    }
    
## PLOT Temperature Data as filled contours

my.colors <- function(n)c("blue",heat.colors(n-1))

#Plot TMIN  
x.year <- seq(startyear, endyear, by=1)
y.day <- seq(1,365,b=1)
TMIN.matrix <- as.matrix(Sta.T.TMIN[,1:365])
summary(TMIN.matrix)
filled.contour(x.year, y.day, TMIN.matrix, color = rainbow,
    plot.title = title(main = paste(StationName, "TMIN", sep=" "),
    xlab = "Year", ylab = "Julian Day, Winter Centered"),
    plot.axes = { axis(1, seq(startyear, endyear, by = 10))
                  axis(2, seq(10,360, by=20)) },
    key.title = title(main="Degrees Farenheit"),
    key.axes = axis(4, seq(-50, 70, by = 10)))
    
#Plot TMAX 
x.year <- seq(startyear, endyear, by=1)
y.day <- seq(1,365,b=1)
TMAX.matrix <- as.matrix(Sta.T.TMAX[,1:365])
summary(TMAX.matrix)
filled.contour(x.year, y.day, TMAX.matrix, color = heat.colors,
    plot.title = title(main = paste(StationName, "TMAX", sep=" "),
    xlab = "Year", ylab = "Julian Day, Calendar Year"),
    plot.axes = { axis(1, seq(startyear, endyear, by = 10))
                  axis(2, seq(10,360, by=10)) },
    key.title = title(main="Degrees Farenheit"),
    key.axes = axis(4, seq(0, 100, by = 10)))
#contour(x.year, y.day, TMAX.matrix, levels = pretty(32,10))