# Functions to find the Parameters used for the functional
# model of Snow. Used to estimate the peak, max, and duration.
find_first <- function(x) {
snow_days <- which(x > 0.2)
first_day <- min(snow_days)
return(first_day)
}
find_peak_index <- function(x, first_day = 97, func = max) {
after_first <- round(x[97:length(x)], digits = 0)
per_max <- func(after_first, na.rm = TRUE)
peak_index <- which(after_first == per_max)[1]
if (is.na(peak_index)) {
return(NA)
} else {
return(peak_index + first_day)
}
}
find_val <- function(x, first_day = 97, func = mean) {
after_first <- round(x[97:length(x)], digits = 0)
per_max <- func(after_first, na.rm = TRUE)
}
find_duration <- function(x) {
beginning <- find_first(x)
peak_index <- which(x == max(x))[1]
zeros <- which(x == 0)
after_cycle <- zeros[which(zeros > peak_index)]
end <- min(after_cycle)
return(end - beginning)
}
find_period_start <- function(dataframe, period) {
index_of_period <- which(dataframe$period == period)
return(min(index_of_period))
}
find_period_peak <- function(dataframe, per, peak_var, func = max) {
period_start <- find_period_start(dataframe, per)
per_df <- dataframe[dataframe$period == per & !is.na(dataframe$period), ]
peak_var_vector <- per_df[, peak_var]
peak_var_index <- find_peak_index(peak_var_vector, first_day = 97, func)
df_index <-period_start + peak_var_index
df_value <- dataframe[df_index, peak_var]
return(list(index = df_index, value = df_value))
}
indices_window <- function(index, window = 5) {
ten_day_index <- seq(from = index - window, to = index + window, by = 1)
return(ten_day_index)
}
# Find the lambda corresponding to the glmnet fit, resticting
# the maximum number of variables allowed in the model.
upper_bound_lambda <- function(model, max_params = 15) {
# Make sure the modle has been built with
# glmnet and has cross validation
stopifnot(class(model) == "cv.glmnet")
# Restrict the acceptable models to the ones which
# have less than the maximum number of term
num_params <- as.numeric(model$nzero)
sparse_enough <- which(num_params < max_params)
mse <- model$cvm
lambda <- model$lambda
# Choose the model with the minimum cross validation
# score out of the remaining models.
min_cv <- min(mse[sparse_enough])
min_index <- which(mse == min_cv)
min_lambda <- lambda[min_index]
return(min_lambda)
}
# Establish a Period variable to seperate the dataset into.
# Then add this variable onto the data-set
add_period <- function(dataframe) {
years <- as.numeric(format.Date(dataframe$date, "%Y"))
year_index <- unique(years)
period <- rep(NA, length = nrow(dataframe))
count <- 0
for (year in year_index) {
# Index the start and end data of the dataset based on a
# data wen there is never swe, JULY 15th!
start_index <- which(dataframe$date == paste0(year, "-07-01"))
end_index <- which(dataframe$date == paste0(year + 1, "-06-30"))
# If the year doesn't have a complete year, then don't make the plot.
if (length(start_index) == 0 | length(end_index) == 0) {
next
}
# Update the cout and assign it to the period vector
count <- count + 1
period[start_index:end_index] <- count
}
dataframe$period <- period
# Remove everything not in one of the periods
dataframe <- dataframe[!is.na(dataframe$period), ]
}
generate_modeling_df <- function(dataframe) {
# Load in the packages for easy summary statistics.
library(tidyr)
library(dplyr)
library(magrittr)
# Use these monthly and yearly variables to obtain the monthly
# averages for our covariates of interest.
monthly_means <- select(dataframe, date, month, period, precip_inc, temp_avg, swe) %>%
group_by(month, period) %>%
summarise(cov_temp = mean(temp_avg, na.rm = TRUE),
cov_prec = mean(precip_inc, na.rm = TRUE))
period_parameters <- select(dataframe, date, year, month, swe, period, precip_accum, precip_inc, temp_avg) %>%
group_by(period) %>%
summarise(peaks = max(swe),
first = find_first(swe))
# Combine the parameter data-frame with the covariate dataframe. Then
# Create a new modeling df which makes it easy to fir linear models
modeling_df <- merge(period_parameters, monthly_means, by = "period")
lm_df <- modeling_df %>%
gather(Var, val, starts_with("cov")) %>%
unite(Var1, Var, month) %>%
spread(Var1, val)
colnames(lm_df) <- gsub("cov_", "", colnames(lm_df))
return(lm_df)
}
compute_monthly_covs <- function(temp_vec, precip_vec, date_vec = NULL) {
# Make sure that the temperature and precipitation vectors
# are the same
stopifnot(length(temp_vec) == length(precip_vec))
# Create a month vector based on the length of the temp and
# precipitation vectors (if it's a leap year)
if (!is.null(date_vec)) {
month_vec <- format.Date(date_vec, "%m")
} else if (length(temp_vec) == 365) {
dates <- seq(as.Date("1981-07-01"), to = as.Date("1982-06-30"), by = 1)
month_vec <- format.Date(dates, "%m")
} else if (length(temp_vec == 366)) {
dates <- seq(as.Date("1979-07-01"), to = as.Date("1980-06-30"), by = 1)
month_vec <- format.Date(dates, "%m")
} else {
stop("Vectors must be either 365 or 366 days!")
}
# Use this monthly vector to obtain the monthly averages for
# both temperature and precipication
month_vec <- as.numeric(month_vec)
dataframe <- as.data.frame(cbind(month_vec, temp_vec, precip_vec))
monthly_temps <- aggregate(formula = temp_vec ~ month_vec, data = dataframe, FUN = mean)
monthly_precip <- aggregate(formula = precip_vec ~ month_vec, data = dataframe, FUN = mean)
# Combine these monthly averages into a one row data-frame
# that can be used for making predictions on a model with the
# same covariate names
newvals <- c(monthly_temps[, 2], monthly_precip[, 2])
df <- data.frame(matrix(newvals, nrow = 1, ncol = length(newvals)))
cov_names <- c(paste0("temp_", 1:12), paste0("prec_", 1:12))
colnames(df) <- cov_names
return(df)
}
predict_swe_peak <- function(temp_vec, precip_vec, date_vec = NULL, model) {
# Compute the monthly covariates using this particular
# temp and precip
predictors <- compute_monthly_covs(temp_vec, precip_vec, date_vec)
# Use these predictors to make a prediction for the
# peak snowfall in this period
peak_pred <- predict(model, newdata = predictors)
return(peak_pred)
}
predict_swe_start <- function(temp_vec, precip_vec, date_vec = NULL, model) {
# Compute the monthly covariates using this particular
# temp and recip vector
predictors <- compute_monthly_covs(temp_vec, precip_vec, date_vec)
# Use these predictors to make a prediction for the
# peak snowfall in this period
first_day_pred <- predict(model, newdata = predictors)
return(first_day_pred)
}
# Function which takes in a precipitation vector, temperature
# vector, and compute what SWE would look like over the course
# of the year.
predict_swe_period <- function(temp_vec, precip_vec, date_vec = NULL, decline_rate = .5,
percent_snow = .9, first_day = 95, peak_model,
start_model, interval = FALSE) {
# Estimate the parameters (first day and peak) for this particular
# period using the two linear models and temp + precip vectors.
first_day <- floor(predict_swe_start(temp_vec, precip_vec, date_vec, start_model))
peak_snow <- predict_swe_peak(temp_vec, precip_vec, date_vec, peak_model)
# Calculate Precip Accumulation after the first day has
# been estimated
zero_prec <- rep(0, first_day - 1)
prec_accum <- cumsum(prec[first_day:length(precip_vec)])
precip_vec <- c(zero_prec, prec_accum)
# Loop through each one of the days and predict the SWE
# at this given day
num_days <- length(temp_vec)
swe_results <- rep(NA, num_days)
current_swe <- 0
peaked <- "no"
for (day in 1:num_days) {
if (current_swe < peak_snow & day < first_day) {
current_swe = 0
swe_results[day] <- current_swe
} else if (current_swe <= peak_snow & day >= first_day & peaked == "no") {
current_swe <- (percent_snow * precip_vec[day])
swe_results[day] <- current_swe
} else if (current_swe > peak_snow & peaked == "no") {
peaked <- "yes"
swe_results[day] <- current_swe
} else if (peaked == "yes" & current_swe > 0) {
current_swe <- current_swe - (decline_rate + rnorm(n = 1, mean = 0, sd = .2))
if (current_swe < 0) {
current_swe = 0
}
swe_results[day] <- current_swe
} else {
current_swe = 0
swe_results[day] <- current_swe
}
}
# Return the vector of SWE for this time period
return(swe_results)
}
convert_names <- function(names) {
# Create a lookup table which maps the differences between
# narcaap and ncep names
narcaap_names <- c("div", "hus", "huss", "rhum", "rhums",
"ta", "tas", "ua", "uas", "va", "vas",
"vort", "zg", "ugeo", "vgeo")
ncep_names <- c("div", "shum", "shum.2m", "rhum", "rhum.2m",
"air", "air.2m", "uwnd", "uwnd.10m", "vwnd",
"vwnd.10m", "vort", "hgt", "ugeo", "vgeo")
lookup_table <- data.frame(narcaap = narcaap_names,
ncep = ncep_names)
# Loop through each one of the column names and substitute
# in the varname from the correct format
new_names <- rep("", times = length(names))
for (name in seq_along(names)) {
cur_name <- names[name]
split_name <- unlist(strsplit(cur_name, "_"))
conv_name_index <- which(lookup_table[, "narcaap"] == split_name[1])
conv_name <- as.character(lookup_table[conv_name_index, "ncep"])
if (exists(split_name[3])) {
new_names[name] <- paste0(conv_name, "_", split_name[2], "_", split_name[3])
} else {
new_names[name] <- paste0(conv_name, "_", split_name[2])
}
}
return(new_names)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.