#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.
####################################