#Description of the included functions is provided in Appendix B.

#Below we present two methods of reading in the data: Path and file specification, or
#Changing the working directory.  After that the commands are identical.

###################################
#Path and file specification.
#In our case, we were working off of a USB key.  Notice the path specification
#and filename designation. Please be sure to reference your appropriate download locationand
#note that the / is required at the end of the path.

path <- "K:/Weather/DLY/"
filename <- "CentralPark.dat";

#Create data field widths;
#This was obtained from the GHCN web site but this format is consistent for all HCN data sets;

fieldWidths <- c(11, 4, 2, 4, rep(c(5, 1, 1, 1), 31));

#Data input

centralPark <- read.fwf(file=paste(path, filename, sep=""), widths=fieldWidths)

###################################
#Change the working directory.
#Alternatively, the working directory can be changed to the appropriate download location by selecting
#File - Change dir... in the R interface.
#This approach may be used if the path for the data files is long (e.g. My Documents or Desktop in Windows).  
#To proceed in this manner, use the following block of code after changing the working directory.  Note that
#this assumes the data file is still named centralPark.dat in the appropriate directory location.  This can
#be changed by substituting the appropriate file name for "centralPark.dat" below:

fieldWidths <- c(11, 4, 2, 4, rep(c(5, 1, 1, 1), 31));

centralPark <- read.fwf("centralPark.dat", widths=fieldWidths)

###################################
#Once the data has been entered, we may check it using the following code:

dim(centralPark)

#The dimensions of the provided data set should be 10,354 rows x 128 columns.

head(centralPark)

#This command should list the first five rows of the data.  The location code is in the first column, followed by the year and month in columns
#two and three and the weather code (column 4).  The remaining 124 (=31 days x 4) columns repeat the Weather value and Quality Control Characters for 
#each possible day.
###################################

#Create the locationCleaning Function

locationCleaning <- function(data, weather="SNWD", month = 12, firstDay = 17, lastDay = 31)
{

#Basic Error Checking
if (firstDay > lastDay)
	stop("Sorry, firstDay must be less than lastDay.")

if (! ((weather == "SNWD") | (weather == "SNOW") | (weather == "TMAX") | (weather == "TMIN") | (weather == "PRCP") ) )
	stop("Sorry, you have entered an invalid weather condition.")

if ((month < 1) || (month > 12))
	stop("Sorry, month must be an integer between 1 and 12.")

if (! ( (is.integer(month)) | (!is.integer(firstDay)) | (!is.integer(lastDay))) )
	stop("Sorry, month, firstDay and lastDay must be integers.")


#Drop the station ID column since it is not needed;
data <- data[,-1]

#Remove all data that is not the required month and weather condition.
keep <- c();

for (i in 1:nrow(data))
{

if (data[i,2] == month)  
{
if (data[i,3] == weather) 
{
keep <- append(keep, i)
}
}

}

dataSNOW <- data[keep,] 

#Keep only the data that is in the appropriate day range.
SNOW <- dataSNOW[,1];

for (i in firstDay:lastDay)
{
SNOW <- cbind(SNOW, as.numeric(as.character(dataSNOW[,(4 + 4*(i-1))])))
}

colnames(SNOW) <- c("year", firstDay:lastDay)

#Create an object of type "SNOW" to ensure functionality with other methods.
class(SNOW) <- "SNOW";
return(SNOW)
}

###############################

#Run locationCleaning and store results into a 'clean' object;

clean <- locationCleaning(data=centralPark, weather="SNWD", month = 12, firstDay = 17, lastDay = 31)
clean;

#The first and last five rows should correspond to those listed in the table on p. 13 of Appendix B. 
###############################


#makeBinary requires an object of class "SNOW" created from locationCleaning
#This function dichotomizes the records according to SNOWCrit

makeBinary <- function(SNOW, SNOWCrit=50)
{
#Basic Error Checking
if (SNOWCrit < 0)
	stop("Sorry, SNOWCrit must be non-negative.")

if (class(SNOW) != "SNOW")
	stop("Sorry, SNOW must be the results of the locationCleaning function.")

#Loop to dichotomize the SNOW data
for (i in 1:nrow(SNOW))
{
for (j in 2:ncol(SNOW))
{

if (!is.na(SNOW[i,j]))
{

#SNOW becomes zero if an insufficient amount of snow is observed.
if ( (SNOW[i,j] < SNOWCrit ) && (SNOW[i,j] >= 0 ) )
{
SNOW[i,j] <- 0;
}

#Otherwise SNOW becomes one.
if (SNOW[i,j] >= SNOWCrit )
{
SNOW[i,j] <- 1;
}

#Note that all missing values (-9999) are changed to NA.
if (SNOW[i,j] == -9999)
{
SNOW[i,j] <- NA;
}

}
}
}

#SNOW01 object is created for the next function.
class(SNOW) <- "SNOW01";
return(SNOW)
}

###############################

#Run makeBinary, everything is now dichotomized according to 50 mm of snow.

SNOWBinary <- makeBinary(SNOW=clean, SNOWCrit=50)
SNOWBinary;

#The first and last five rows should correspond to those listed in the table on p. 14 of Appendix B. 
###############################

#makeP requires an object of class "SNOW01" from makeBinary
#This function estimates the transition matrix P for model 2 and
#the overall independence probability in model 1;

makeP <- function(SNOW01)
{

#Basic Error Checking
if (class(SNOW01) != "SNOW01")
	stop("Sorry, this function requires the results of the makeBinary function.")

#Initialize return objects
X <- NULL;
X$PSNOW <- matrix(0, nrow=2,ncol=2)

#Perform maximum likelihood estimation for Model 2
for (i in 1:nrow(SNOW01))
{
for (j in 2:(ncol(SNOW01) -1))
if ( (!is.na(SNOW01[i,j])) && (!is.na(SNOW01[i,(j + 1)])) )
{

if ( (SNOW01[i,j] == 0) && (SNOW01[i,j+1] == 0) )
{
X$PSNOW[1,1] <- X$PSNOW[1,1] + 1;
}

if ( (SNOW01[i,j] == 0) && (SNOW01[i,j+1] == 1) )
{
X$PSNOW[1,2] <- X$PSNOW[1,2] + 1;
}

if ( (SNOW01[i,j] == 1) && (SNOW01[i,j+1] == 0) )
{
X$PSNOW[2,1] <- X$PSNOW[2,1] + 1;
}

if ( (SNOW01[i,j] == 1) && (SNOW01[i,j+1] == 1) )
{
X$PSNOW[2,2] <- X$PSNOW[2,2] + 1;
}

}

}

#Normalize X$PSNOW
for (i in 1:2)
{
X$PSNOW[i,] <- X$PSNOW[i,]/sum(X$PSNOW[i,])
}


#The overall probability seeing a significant snowdepth
#on any day, should it be of interest.

X$PSNOWInd <- mean(SNOW01[,-1], na.rm=TRUE)

#PSNOW object is created for the prediction methods.
class(X) <- "PSNOW"
return(X);
}

###################################

#Estimate P and store results in Probs;

Probs <- makeP(SNOW01=SNOWBinary)
Probs;

#These values should correspond with the value for P on p. 14.
####################################