In scoyoc/climateAnalyzeR: Import Data from ClimateAnalyzer.org into R
\raggedright
#-- Packages
# devtools::install_github("scoyoc/climateAnalyzeR")
library("climateAnalyzeR")
library("tidyr")
#-- Functions
# Number Contraction --
# Prints the appropriate number contraction for a given number.
number_contraction <- function(x){
if(is.na(x)){
y = FALSE
} else if(x == 11 | x == 12 | x == 13){
y = paste0(x, "th")
} else if(x - round(x, -1) == 1){
y = paste0(x, "st")
} else if(x - round(x, -1) == 2){
y = paste0(x, "nd")
} else if(x - round(x, -1) == 3){
y = paste0(x, "rd")
} else{
y = paste0(x, "th")
}
return(y)
}
# Glue Months --
# Returns a string of months with comma's between each month.
glue_mths <- function(mths){
if(length(mths) == 1){
mths
} else if(length(mths) == 2){
glue::glue("{mths[1]} and {mths[2]}")
} else{
glue::glue("{paste(mths[1:length(mths) - 1], collapse = ', ')}, and {mths[length(mths)]}")
}
}
# Round with Zeros --
# Rounds numbers and includes zeros for significant digits
round_with_zeros <- function(x, d = 2){
format(round(x, digits = d), nsmall = d)
}
#-- Load data frame of months
data(month_df)
month_df$month_name <- factor(month_df$month_name, month.abb)
#-- Retrieve station metadata
my_station <- stations(station_id = station_id) %>%
dplyr::distinct()
yrs_avail <- as.numeric(stringr::str_split(my_station$years_avail,
pattern = ";", simplify = T))
my_station <- my_station %>%
dplyr::mutate("label" = paste(min(yrs_avail, na.rm = T),
max(yrs_avail, na.rm = T), sep = "-"),
"water_yr" = my_year,
"name" = station_name)
#-- 30-year Normals
normals_30yr <- normals(station_id = station_id) %>%
tidyr::spread("element", "value")
annual_30yr <- normals_30yr %>%
dplyr::filter(month == "Annual") %>%
dplyr::mutate("TEMP_midpt" = ((TMAX - TMIN) / 2) + TMIN)
month_30yr <- normals_30yr %>%
dplyr::filter(month %in% month.abb) %>%
dplyr::mutate("month" = month_df$month_name) %>%
dplyr::rename("month_name" = month,
"prcp_30yr" = PRCP,
"tmax_30yr" = TMAX,
"tmin_30yr" = TMIN) %>%
dplyr::arrange(month_name) %>%
dplyr::mutate("prcp_accum_30yr" = cumsum(tidyr::replace_na(prcp_30yr, 0)),
"temp_midpt_30yr" =((tmax_30yr - tmin_30yr) / 2) + tmin_30yr)
#-- Annual Data
annual_dat <- import_annual(station_id = station_id,
start_year = min(yrs_avail, na.rm = T),
end_year = lubridate::year(lubridate::today()),
screen_blanks = 'true',
remove_missing = TRUE) %>%
dplyr::select(year | contains("cy")) %>%
dplyr::mutate("temp_midpt" = ifelse(is.na(tmax_cy) | is.na(tmin_cy), NA,
((tmax_cy - tmin_cy) / 2) + tmin_cy),
"prcp_driest" = dplyr::min_rank(prcp_cy),
"prcp_wettest" = dplyr::min_rank(-prcp_cy),
"temp_warmest" = dplyr::min_rank(-temp_midpt),
"temp_coldest" = dplyr::min_rank(temp_midpt),
"prcp_depart" = prcp_cy - annual_30yr$PRCP,
"temp_depart" = temp_midpt - annual_30yr$TEMP_midpt,
"tmax_depart" = tmax_cy - annual_30yr$TMAX,
"tmin_depart" = tmin_cy - annual_30yr$TMIN)
#-- Monthly Data
# Monthly weather data
mth_dat <- import_monthly(station_id = station_id,
start_year = min(yrs_avail, na.rm = T),
end_year = lubridate::year(lubridate::today())) %>%
dplyr::left_join(month_df) %>%
dplyr::arrange(year, month) %>%
dplyr::group_by(year) %>%
dplyr::mutate("prcp_accum" = cumsum(tidyr::replace_na(prcp, 0))) %>%
dplyr::mutate("temp_midpt" = ifelse(is.na(tmax) | is.na(tmin), NA,
((tmax - tmin) / 2) + tmin)) %>%
dplyr::group_by(water_month) %>%
dplyr::mutate("prcp_driest" = dplyr::min_rank(prcp),
"temp_warmest" = dplyr::min_rank(-temp_midpt))
# Monthly departure
mth_depart <- import_departure(station_id = station_id,
min(yrs_avail, na.rm = T),
end_year = lubridate::year(lubridate::today())) %>%
dplyr::left_join(month_df) %>%
dplyr::left_join(dplyr::left_join(mth_dat, month_30yr) %>%
dplyr::mutate("prcp_accum_depart" = prcp_accum -
prcp_accum_30yr) %>%
dplyr::ungroup() %>%
dplyr::select(year, month_name, prcp_accum_depart))
# Isolate month
my_mth <- dplyr::filter(mth_dat, year == my_year) %>%
dplyr::ungroup() %>%
dplyr::select(month_name, year) %>%
dplyr::mutate("month_name" = as.character(month_name))
my_mth <- my_mth[nrow(my_mth), ]
# Filter prcp and temp, then rank and reduce for referencing in text
mth_rank <- dplyr::filter(mth_dat, month_name %in% my_mth$month_name) %>%
dplyr::ungroup() %>%
dplyr::select("year", "prcp", "temp_midpt", "tmax", "tmin") %>%
tidyr::gather("element", "value", prcp:tmin) %>%
dplyr::group_by(element) %>%
dplyr::mutate("wet_warm" = dplyr::min_rank(-value),
"dry_cool" = dplyr::min_rank(value),
"n" = dplyr::n()) %>%
dplyr::filter(year == my_year) %>%
dplyr::left_join(
dplyr::filter(month_30yr, month_name %in% my_mth$month_name) %>%
dplyr::ungroup() %>%
dplyr::select(prcp_30yr, temp_midpt_30yr, tmax_30yr, tmin_30yr) %>%
dplyr::rename("prcp" = prcp_30yr,
"temp_midpt" = temp_midpt_30yr,
"tmax" = tmax_30yr,
"tmin" = tmin_30yr) %>%
tidyr::gather("element", "norms_30yr")
) %>%
dplyr::mutate("departure" = value - norms_30yr)
my_mth_prcp <- dplyr::filter(mth_rank, element == "prcp") %>%
dplyr::left_join(
dplyr::filter(mth_dat, month_name %in% my_mth$month_name &
year == my_year) %>%
dplyr::ungroup() %>%
dplyr::select(prcp_accum) %>%
dplyr::mutate("element" = "prcp")
) %>%
dplyr::left_join(
dplyr::filter(mth_depart, month_name %in% my_mth$month_name &
year == my_year) %>%
dplyr::ungroup() %>%
dplyr::select(prcp_accum_depart) %>%
dplyr::mutate("element" = "prcp") %>%
dplyr::rename("accum_depart" = prcp_accum_depart)
) %>%
dplyr::left_join(
dplyr::filter(month_30yr, month_name %in% my_mth$month_name) %>%
dplyr::ungroup() %>%
dplyr::select(prcp_accum_30yr) %>%
dplyr::mutate("element" = "prcp")
)
my_mth_temp <- dplyr::filter(mth_rank, element == "temp_midpt")
my_mth_tmax <- dplyr::filter(mth_rank, element == "tmax")
my_mth_tmin <- dplyr::filter(mth_rank, element == "tmin")
# Consolidate text references
text_refs <- my_station %>%
dplyr::bind_cols(dplyr::filter(annual_dat, year == my_year)) %>%
dplyr::mutate("timespan" = max(yrs_avail, na.rm = T) - min(yrs_avail, na.rm = T))
#-- Model annual trends
prcp_lm <- summary(lm(prcp_cy ~ year, data = annual_dat))$coefficients[2, 1]
temp_lm <- summary(lm(temp_midpt ~ year, data = annual_dat))$coefficients[2, 1]
#-- Subset months that exceed normals
tmax_exceed <- dplyr::filter(mth_depart, year == my_year &
tmax_depart > 0)$month_name
tmin_exceed <- dplyr::filter(mth_depart, year == my_year &
tmin_depart > 0)$month_name
prcp_above <- dplyr::filter(mth_depart, year == my_year &
prcp_pctavg > 100)$month_name
prcp_below <- dplyr::filter(mth_depart, year == my_year &
prcp_pctavg < 100)$month_name
# Error evaluation
error_eval <- dplyr::filter(mth_dat, year == my_year) %>%
dplyr::ungroup() %>%
dplyr::select(c("prcp", "temp_midpt")) %>%
dplyr::summarise_all(function(x)(sum(is.na(x))))
# Evaluate data for conditional text
# Calculate percentiles
prcp_quant <- quantile(dplyr::filter(mth_dat, month_name %in%
my_mth$month_name)$prcp,
c(0.25, 0.75), na.rm = TRUE)
temp_quant <- quantile(dplyr::filter(mth_dat, month_name %in%
my_mth$month_name)$temp_midpt,
c(0.25, 0.75), na.rm = TRUE)
# Assign condition
prcp_text <- if(error_eval$prcp >= 2 | is.na(my_mth_prcp$value)){
"error"
} else if(my_mth_prcp$value <= prcp_quant[1]){
"dry"
} else if(my_mth_prcp$value > prcp_quant[1] &
my_mth_prcp$value < prcp_quant[2]){
"avg"
} else if(my_mth_prcp$value >= prcp_quant[2]){
"wet"
} else(FALSE)
temp_text <- if(error_eval$temp_midpt >= 2 | is.na(my_mth_temp$value)){
'error'
} else if(my_mth_temp$value <= temp_quant[1]){
"cool"
} else if(my_mth_temp$value > temp_quant[1] &
my_mth_temp$value < temp_quant[2]){
"avg"
} else if(my_mth_temp$value >= temp_quant[2]){
"warm"
} else(FALSE)
This document summarizes temperature and precipitation anomalies relative to 30-year (1981-2010) averages for r my_station$name for calendar year r my_year through r month.name[match(my_mth$month_name, month.abb)].
The data used in these analyses are from the NOAA Global Historical Climatology Network (GHCN) Cooperative Observer Network (COOP) and were downloaded from ClimateAnalyzer.org on r format(as.Date(lubridate::today(), format = "%Y-%m-%d"), "%d %B, %Y") using R (version r paste(R.version$major, R.version$minor, sep = "."), R Core Team, r R.version$year) and the climateAnalyzeR package (ver r packageVersion("climateAnalyzeR"), Van Scoyoc, r lubridate::year(packageDate("climateAnalyzeR"))).
The data used for these analyses include the daily high temperature (TMAX), daily low temperature (TMIN), and daily precipitation totals (PRCP).
This is an automated summary, and all results are provisional.
# Identify location
extent <- dplyr::select(my_station, lon, lat)
# Expand to desired scale
extent[2, ] <- c(extent[1, 1] + 0.04, extent[1, 2] + 0.025)
extent[3, ] <- c(extent[1, 1] - 0.04, extent[1, 2] - 0.025)
# Create box to pull basemap from Google
my_box <- ggmap::make_bbox(lon = extent$lon, lat = extent$lat, f = 0.1)
# Pull basemap from Google
my_map <- ggmap::get_map(location = my_box, maptype = "terrain",
source = "stamen", zoom = 13, force = TRUE,
messaging = FALSE)
# Finish Map
my_map <- ggmap::ggmap(my_map) +
ggplot2::geom_point(data = my_station, ggplot2::aes(x = lon, y = lat),
color = "red", size = 4) +
shadowtext::geom_shadowtext(data = my_station,
mapping = ggplot2::aes(x = lon, y = lat,
label = name),
color = "black", bg.colour = "white", size = 3,
hjust = -0.1, vjust = -0.5) +
ggplot2::labs(x = "Longitude", y = "Latitude")
my_map
\newpage
Temperature
# Create text objects for temp discussion
temp_rank_text <- if(temp_text == "error"){
glue::glue("There are too many missing records to meaningfully summarize temperature data from this station through {month.name[match(my_mth$month_name, month.abb)]} {my_mth$year}.")
} else if(!is.na(text_refs$temp_coldest) & text_refs$temp_coldest == 1){
glue::glue("{month.name[match(my_mth$month_name, month.abb)]} {my_mth$year} was the coldest {month.name[match(my_mth$month_name, month.abb)]} on record at {my_station$name}.")
} else if(temp_text == "cool"){
glue::glue("{month.name[match(my_mth$month_name, month.abb)]} {my_mth$year} was the {number_contraction(my_mth_temp$dry_cool)} coolest {month.name[match(my_mth$month_name, month.abb)]} in the {text_refs$timespan}-year record at {my_station$name}.")
} else if(temp_text == "avg"){
glue::glue("Temperatures were normal in {month.name[match(my_mth$month_name, month.abb)]} {my_mth$year}. This month was the {number_contraction(my_mth_temp$wet_warm)} warmest {month.name[match(my_mth$month_name, month.abb)]} in the {text_refs$timespan}-year record at {my_station$name}.")
} else if(!is.na(text_refs$temp_warmest) & text_refs$temp_warmest == 1){
glue::glue("{month.name[match(my_mth$month_name, month.abb)]} {my_mth$year} was the hottest {month.name[match(my_mth$month_name, month.abb)]} on record at {my_station$name}.")
} else if(temp_text == "warm"){
glue::glue("{month.name[match(my_mth$month_name, month.abb)]} {my_mth$year} was the {number_contraction(my_mth_temp$wet_warm)} warmest {month.name[match(my_mth$month_name, month.abb)]} in the {text_refs$timespan}-year record at {my_station$name}.")
}
temp_mth_text <- glue::glue("The average temperature for {month.name[match(my_mth$month_name, month.abb)]} was {round_with_zeros(my_mth_temp$value)}°F which is {round_with_zeros(abs(my_mth_temp$departure))}°F {ifelse(my_mth_temp$departure > 0, 'above', 'below')} the 30-year average ({round_with_zeros(my_mth_temp$norms_30yr)}°F).")
temp_trend_text <- glue::glue("This summary suggests that temperatures are {ifelse(temp_lm > 0, 'increasing', 'decreasing')} {round_with_zeros(abs(temp_lm * 10))}°F per decade (Figure 1A).")
# Text about monthly exceedences
mths_so_far <- month_df[month_df$month_name == my_mth$month_name, ]$water_month
temp_tmax_text <- if(length(tmax_exceed) == 0){
glue::glue("Every monthly average TMAX has been below normal this water year (Figure 1B).")
} else if(length(tmax_exceed) == mths_so_far){
glue::glue("Every monthly average TMAX has been above normal this water year (Figure 1B).")
} else(
glue::glue("The monthly average TMAX has been {ifelse(length(tmax_exceed) >= (mths_so_far / 2), 'above', 'below')} normal most of the year and has exceeded the 30-year monthly average {length(tmax_exceed)} times ({glue_mths(tmax_exceed)}; Figure 1B).")
)
temp_tmin_text <- if(length(tmin_exceed) == 0){
glue::glue("Every monthly average TMIN has been below normal this water year (Figure 1B).")
} else if(length(tmin_exceed) == mths_so_far){
glue::glue("Every monthly average TMIN has been above normal this water year (Figure 1B).")
} else(
glue::glue("The monthly average TMIN has been {ifelse(length(tmin_exceed) >= (mths_so_far / 2), 'above', 'below')} normal most of the year and has exceeded the 30-year monthly average {length(tmin_exceed)} times ({glue_mths(tmin_exceed)}; Figure 1B).")
)
if(temp_text == "error"){
temp_error_trend_text <- glue::glue("However, this summary suggests that temperatures are {ifelse(temp_lm > 0, 'increasing', 'decreasing')} {round_with_zeros(abs(temp_lm * 10))}°F per decade (Figure 1A).")
}
# Text is printed in PDF if `temp_text %in% c("cool", "avg", "warm")` returns TRUE
cat(temp_rank_text,
temp_mth_text,
temp_trend_text,
temp_tmax_text,
temp_tmin_text,
sep = " ")
# Text for errors in monthly data set. Annual data set is good.
cat(temp_rank_text,
"Examine the monthly data set on ClimateAnalyzer.org for more information.",
temp_error_trend_text,
sep = " ")
# Text is printed if temp_text == FALSE
cat("Error: Temperature data could not be evaluated. Examine tabular data sets on ClimateAnalyzer.org.")
#-- Average Annual Temps
t1 <- annual_figure(x_var = annual_dat$year,
y_var = annual_dat$temp_midpt,
my_year = my_year,
normal = annual_30yr$TEMP_midpt,
reference_period = "1981-2010",
area_color = "orange",
line_color = "red3",
my_title = "A. Average Annual Temperature",
my_ylab = "Temperature (°F)")
#-- Average WY TMAX/TMIN
t2 <- monthly_figure(x_var = mth_depart$month_name,
y_var = mth_depart$tmax_depart,
year_group = mth_depart$year,
y_intercept = 0,
my_year = my_year,
line_color = "red3",
my_title = "B. Average Monthly TMAX",
my_ylab = "Departure (°F)") +
ggplot2::theme(axis.text.x = ggplot2::element_blank())
#-- TMAX/TMIN Departure
t3 <- monthly_figure(x_var = mth_depart$month_name,
y_var = mth_depart$tmin_depart,
year_group = mth_depart$year,
y_intercept = 0,
my_year = my_year,
line_color = "red3",
my_title = "C. Average Monthly TMIN",
my_ylab = "Departure (°F)")
cowplot::plot_grid(t1, t2, t3, labels = c(), ncol = 1, align = "v",
vjust = 1, scale = 1)
\newpage
Precipitation
# Create text objects for prcp discussion
prcp_rank_text <- if(prcp_text == "error"){
glue::glue("There are too many missing records to meaningfully summarize precipitation data from this station for {month.name[match(my_mth$month_name, month.abb)]} {my_mth$year}.")
} else if(prcp_text == "dry" & !is.na(text_refs$prcp_driest) &
text_refs$prcp_driest == 1){
glue::glue("{month.name[match(my_mth$month_name, month.abb)]} {my_mth$year} was the driest {month.name[match(my_mth$month_name, month.abb)]} on record at {my_station$name}.")
} else if(prcp_text == "dry"){
glue::glue("{month.name[match(my_mth$month_name, month.abb)]} {my_mth$year} was the {number_contraction(my_mth_prcp$dry_cool)} driest {month.name[match(my_mth$month_name, month.abb)]} in the {text_refs$timespan}-year record at {my_station$name}.")
} else if(prcp_text == "avg"){
glue::glue("Precipitation was normal in {month.name[match(my_mth$month_name, month.abb)]} {my_mth$year}. This month was the {number_contraction(my_mth_prcp$dry_cool)} driest water year in the {text_refs$timespan}-year record at {my_station$name}.")
} else if(prcp_text == "wet" & !is.na(text_refs$prcp_wettest) &
text_refs$prcp_wettest == 1){
glue::glue("{month.name[match(my_mth$month_name, month.abb)]} {my_mth$year} was the wettest {month.name[match(my_mth$month_name, month.abb)]} on record at {my_station$name}.")
} else if(prcp_text == "wet"){
glue::glue("{month.name[match(my_mth$month_name, month.abb)]} {my_mth$year} was the {number_contraction(my_mth_prcp$wet_warm)} wettest {month.name[match(my_mth$month_name, month.abb)]} in the {text_refs$timespan}-year record at {my_station$name}.")
}
prcp_annual_text <- glue::glue("Accumulated precipitation is {round_with_zeros(my_mth_prcp$prcp_accum)} inches so far this year which is {round_with_zeros(abs(my_mth_prcp$accum_depart))} inches {ifelse(my_mth_prcp$accum_depart < 0, 'below', 'above')} the 30-year average ({round_with_zeros(my_mth_prcp$prcp_accum_30yr)} inches).")
prcp_trend_text <- glue::glue("This summary suggests that precipitation is {ifelse(prcp_lm > 0, 'increasing', 'decreasing')} at a rate of {round_with_zeros(abs(prcp_lm * 10))} inches per decade (Figure 2A).")
prcp_condition_text <- if(prcp_text == "dry" & length(prcp_above) == 0){
glue::glue("Every month has received below average precipitation so far this year (Figures 2B). The accumulated precipitation is {ifelse(my_mth_prcp$accum_depart < 0, 'below', 'above')} the 30-year average up to this point (Figure 2C).")
} else if(prcp_text == "dry"){
glue::glue("{glue_mths(prcp_above)} received above average precipitation (Figure 2B) but have not been enought to compensate for the {length(prcp_below)} months that were below average ({glue_mths(prcp_below)}; Figure 2C).")
} else if(prcp_text == "avg"){
glue::glue("{glue_mths(prcp_above)} received above average precipitation and {glue_mths(prcp_below)} received below average precipitation (Figure 2B). {my_year} is {ifelse(my_mth_prcp$accum_depart < 0, 'below', 'above')} the 30-year average to this point (Figure 2C).")
} else if(prcp_text == "wet" & length(prcp_below) == 0){
glue::glue("Every month has received above average precipitation so far this year (Figures 2B). This station has received {round_with_zeros(my_mth_prcp$accum_depart)} inches of precipitation so far this year (Figure 2C).")
} else if(prcp_text == "wet"){
glue::glue("{glue_mths(prcp_above)} received above average precipitation (Figure 2B). The accumulated precipitation is {ifelse(my_mth_prcp$accum_depart < 0, 'below', 'above')} the 30-year average up to this point (Figure 2C).")
}
if(prcp_text == "error"){
prcp_error_trend_text <- glue::glue("However, this summary suggests that precipitation is {ifelse(prcp_lm > 0, 'increasing', 'decreasing')} by {round_with_zeros(abs(prcp_lm * 10))} inches per decade (Figure 1A).")
}
cat(
prcp_rank_text,
prcp_annual_text,
prcp_trend_text,
prcp_condition_text,
sep = " ")
# Text for errors in monthly and annual data sets.
cat(prcp_rank_text,
"Examine annual and monthly data sets on ClimateAnalyzer.org for more information.",
prcp_error_trend_text,
sep = " ")
# Text is printed if prcp_text == FALSE
cat("Error: precipitation data could not be evaluated. Examine tabular data set on ClimateAnalyzer.org.")
#-- Water year precipitation totals
p1 <- annual_figure(x_var = annual_dat$year,
y_var = annual_dat$prcp_cy,
my_year = my_year,
normal = annual_30yr$PRCP,
reference_period = "1981-2010",
area_color = "green3",
line_color = "darkgreen",
my_title = "A. Precipitation Totals",
my_ylab = "Precipitation (inches)")
#-- Monthly precipitation totals
p2 <- monthly_figure(x_var = mth_depart$month_name,
y_var = mth_depart$prcp_pctavg,
year_group = mth_depart$year,
y_intercept = 100,
my_year = my_year,
line_color = "darkgreen",
my_title = "B. Monthly Departure Realtive to the Historic Record",
my_ylab = "Percent Average (%)") +
ggplot2::theme(axis.text.x = ggplot2::element_blank())
#-- Water year accumulation
p3 <- monthly_figure(x_var = mth_depart$month_name,
y_var = mth_depart$prcp_accum_depart,
year_group = mth_depart$year,
y_intercept = 0,
my_year = my_year,
line_color = "darkgreen",
my_title = "C. Accumulated Departure Realtive to the Historic Record",
my_ylab = "Departure (inches)")
# Using cowplot to arrange figures
cowplot::plot_grid(p1, p2, p3, labels = c(""), ncol = 1, align = "v")
scoyoc/climateAnalyzeR documentation built on April 19, 2023, 9:57 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
\raggedright
#-- Packages # devtools::install_github("scoyoc/climateAnalyzeR") library("climateAnalyzeR") library("tidyr") #-- Functions # Number Contraction -- # Prints the appropriate number contraction for a given number. number_contraction <- function(x){ if(is.na(x)){ y = FALSE } else if(x == 11 | x == 12 | x == 13){ y = paste0(x, "th") } else if(x - round(x, -1) == 1){ y = paste0(x, "st") } else if(x - round(x, -1) == 2){ y = paste0(x, "nd") } else if(x - round(x, -1) == 3){ y = paste0(x, "rd") } else{ y = paste0(x, "th") } return(y) } # Glue Months -- # Returns a string of months with comma's between each month. glue_mths <- function(mths){ if(length(mths) == 1){ mths } else if(length(mths) == 2){ glue::glue("{mths[1]} and {mths[2]}") } else{ glue::glue("{paste(mths[1:length(mths) - 1], collapse = ', ')}, and {mths[length(mths)]}") } } # Round with Zeros -- # Rounds numbers and includes zeros for significant digits round_with_zeros <- function(x, d = 2){ format(round(x, digits = d), nsmall = d) } #-- Load data frame of months data(month_df) month_df$month_name <- factor(month_df$month_name, month.abb)
#-- Retrieve station metadata my_station <- stations(station_id = station_id) %>% dplyr::distinct() yrs_avail <- as.numeric(stringr::str_split(my_station$years_avail, pattern = ";", simplify = T)) my_station <- my_station %>% dplyr::mutate("label" = paste(min(yrs_avail, na.rm = T), max(yrs_avail, na.rm = T), sep = "-"), "water_yr" = my_year, "name" = station_name) #-- 30-year Normals normals_30yr <- normals(station_id = station_id) %>% tidyr::spread("element", "value") annual_30yr <- normals_30yr %>% dplyr::filter(month == "Annual") %>% dplyr::mutate("TEMP_midpt" = ((TMAX - TMIN) / 2) + TMIN) month_30yr <- normals_30yr %>% dplyr::filter(month %in% month.abb) %>% dplyr::mutate("month" = month_df$month_name) %>% dplyr::rename("month_name" = month, "prcp_30yr" = PRCP, "tmax_30yr" = TMAX, "tmin_30yr" = TMIN) %>% dplyr::arrange(month_name) %>% dplyr::mutate("prcp_accum_30yr" = cumsum(tidyr::replace_na(prcp_30yr, 0)), "temp_midpt_30yr" =((tmax_30yr - tmin_30yr) / 2) + tmin_30yr) #-- Annual Data annual_dat <- import_annual(station_id = station_id, start_year = min(yrs_avail, na.rm = T), end_year = lubridate::year(lubridate::today()), screen_blanks = 'true', remove_missing = TRUE) %>% dplyr::select(year | contains("cy")) %>% dplyr::mutate("temp_midpt" = ifelse(is.na(tmax_cy) | is.na(tmin_cy), NA, ((tmax_cy - tmin_cy) / 2) + tmin_cy), "prcp_driest" = dplyr::min_rank(prcp_cy), "prcp_wettest" = dplyr::min_rank(-prcp_cy), "temp_warmest" = dplyr::min_rank(-temp_midpt), "temp_coldest" = dplyr::min_rank(temp_midpt), "prcp_depart" = prcp_cy - annual_30yr$PRCP, "temp_depart" = temp_midpt - annual_30yr$TEMP_midpt, "tmax_depart" = tmax_cy - annual_30yr$TMAX, "tmin_depart" = tmin_cy - annual_30yr$TMIN) #-- Monthly Data # Monthly weather data mth_dat <- import_monthly(station_id = station_id, start_year = min(yrs_avail, na.rm = T), end_year = lubridate::year(lubridate::today())) %>% dplyr::left_join(month_df) %>% dplyr::arrange(year, month) %>% dplyr::group_by(year) %>% dplyr::mutate("prcp_accum" = cumsum(tidyr::replace_na(prcp, 0))) %>% dplyr::mutate("temp_midpt" = ifelse(is.na(tmax) | is.na(tmin), NA, ((tmax - tmin) / 2) + tmin)) %>% dplyr::group_by(water_month) %>% dplyr::mutate("prcp_driest" = dplyr::min_rank(prcp), "temp_warmest" = dplyr::min_rank(-temp_midpt)) # Monthly departure mth_depart <- import_departure(station_id = station_id, min(yrs_avail, na.rm = T), end_year = lubridate::year(lubridate::today())) %>% dplyr::left_join(month_df) %>% dplyr::left_join(dplyr::left_join(mth_dat, month_30yr) %>% dplyr::mutate("prcp_accum_depart" = prcp_accum - prcp_accum_30yr) %>% dplyr::ungroup() %>% dplyr::select(year, month_name, prcp_accum_depart))
# Isolate month my_mth <- dplyr::filter(mth_dat, year == my_year) %>% dplyr::ungroup() %>% dplyr::select(month_name, year) %>% dplyr::mutate("month_name" = as.character(month_name)) my_mth <- my_mth[nrow(my_mth), ] # Filter prcp and temp, then rank and reduce for referencing in text mth_rank <- dplyr::filter(mth_dat, month_name %in% my_mth$month_name) %>% dplyr::ungroup() %>% dplyr::select("year", "prcp", "temp_midpt", "tmax", "tmin") %>% tidyr::gather("element", "value", prcp:tmin) %>% dplyr::group_by(element) %>% dplyr::mutate("wet_warm" = dplyr::min_rank(-value), "dry_cool" = dplyr::min_rank(value), "n" = dplyr::n()) %>% dplyr::filter(year == my_year) %>% dplyr::left_join( dplyr::filter(month_30yr, month_name %in% my_mth$month_name) %>% dplyr::ungroup() %>% dplyr::select(prcp_30yr, temp_midpt_30yr, tmax_30yr, tmin_30yr) %>% dplyr::rename("prcp" = prcp_30yr, "temp_midpt" = temp_midpt_30yr, "tmax" = tmax_30yr, "tmin" = tmin_30yr) %>% tidyr::gather("element", "norms_30yr") ) %>% dplyr::mutate("departure" = value - norms_30yr) my_mth_prcp <- dplyr::filter(mth_rank, element == "prcp") %>% dplyr::left_join( dplyr::filter(mth_dat, month_name %in% my_mth$month_name & year == my_year) %>% dplyr::ungroup() %>% dplyr::select(prcp_accum) %>% dplyr::mutate("element" = "prcp") ) %>% dplyr::left_join( dplyr::filter(mth_depart, month_name %in% my_mth$month_name & year == my_year) %>% dplyr::ungroup() %>% dplyr::select(prcp_accum_depart) %>% dplyr::mutate("element" = "prcp") %>% dplyr::rename("accum_depart" = prcp_accum_depart) ) %>% dplyr::left_join( dplyr::filter(month_30yr, month_name %in% my_mth$month_name) %>% dplyr::ungroup() %>% dplyr::select(prcp_accum_30yr) %>% dplyr::mutate("element" = "prcp") ) my_mth_temp <- dplyr::filter(mth_rank, element == "temp_midpt") my_mth_tmax <- dplyr::filter(mth_rank, element == "tmax") my_mth_tmin <- dplyr::filter(mth_rank, element == "tmin") # Consolidate text references text_refs <- my_station %>% dplyr::bind_cols(dplyr::filter(annual_dat, year == my_year)) %>% dplyr::mutate("timespan" = max(yrs_avail, na.rm = T) - min(yrs_avail, na.rm = T)) #-- Model annual trends prcp_lm <- summary(lm(prcp_cy ~ year, data = annual_dat))$coefficients[2, 1] temp_lm <- summary(lm(temp_midpt ~ year, data = annual_dat))$coefficients[2, 1] #-- Subset months that exceed normals tmax_exceed <- dplyr::filter(mth_depart, year == my_year & tmax_depart > 0)$month_name tmin_exceed <- dplyr::filter(mth_depart, year == my_year & tmin_depart > 0)$month_name prcp_above <- dplyr::filter(mth_depart, year == my_year & prcp_pctavg > 100)$month_name prcp_below <- dplyr::filter(mth_depart, year == my_year & prcp_pctavg < 100)$month_name # Error evaluation error_eval <- dplyr::filter(mth_dat, year == my_year) %>% dplyr::ungroup() %>% dplyr::select(c("prcp", "temp_midpt")) %>% dplyr::summarise_all(function(x)(sum(is.na(x)))) # Evaluate data for conditional text # Calculate percentiles prcp_quant <- quantile(dplyr::filter(mth_dat, month_name %in% my_mth$month_name)$prcp, c(0.25, 0.75), na.rm = TRUE) temp_quant <- quantile(dplyr::filter(mth_dat, month_name %in% my_mth$month_name)$temp_midpt, c(0.25, 0.75), na.rm = TRUE) # Assign condition prcp_text <- if(error_eval$prcp >= 2 | is.na(my_mth_prcp$value)){ "error" } else if(my_mth_prcp$value <= prcp_quant[1]){ "dry" } else if(my_mth_prcp$value > prcp_quant[1] & my_mth_prcp$value < prcp_quant[2]){ "avg" } else if(my_mth_prcp$value >= prcp_quant[2]){ "wet" } else(FALSE) temp_text <- if(error_eval$temp_midpt >= 2 | is.na(my_mth_temp$value)){ 'error' } else if(my_mth_temp$value <= temp_quant[1]){ "cool" } else if(my_mth_temp$value > temp_quant[1] & my_mth_temp$value < temp_quant[2]){ "avg" } else if(my_mth_temp$value >= temp_quant[2]){ "warm" } else(FALSE)
This document summarizes temperature and precipitation anomalies relative to 30-year (1981-2010) averages for r my_station$name for calendar year r my_year through r month.name[match(my_mth$month_name, month.abb)].
The data used in these analyses are from the NOAA Global Historical Climatology Network (GHCN) Cooperative Observer Network (COOP) and were downloaded from ClimateAnalyzer.org on r format(as.Date(lubridate::today(), format = "%Y-%m-%d"), "%d %B, %Y") using R (version r paste(R.version$major, R.version$minor, sep = "."), R Core Team, r R.version$year) and the climateAnalyzeR package (ver r packageVersion("climateAnalyzeR"), Van Scoyoc, r lubridate::year(packageDate("climateAnalyzeR"))).
The data used for these analyses include the daily high temperature (TMAX), daily low temperature (TMIN), and daily precipitation totals (PRCP).
This is an automated summary, and all results are provisional.
# Identify location extent <- dplyr::select(my_station, lon, lat) # Expand to desired scale extent[2, ] <- c(extent[1, 1] + 0.04, extent[1, 2] + 0.025) extent[3, ] <- c(extent[1, 1] - 0.04, extent[1, 2] - 0.025) # Create box to pull basemap from Google my_box <- ggmap::make_bbox(lon = extent$lon, lat = extent$lat, f = 0.1) # Pull basemap from Google my_map <- ggmap::get_map(location = my_box, maptype = "terrain", source = "stamen", zoom = 13, force = TRUE, messaging = FALSE) # Finish Map my_map <- ggmap::ggmap(my_map) + ggplot2::geom_point(data = my_station, ggplot2::aes(x = lon, y = lat), color = "red", size = 4) + shadowtext::geom_shadowtext(data = my_station, mapping = ggplot2::aes(x = lon, y = lat, label = name), color = "black", bg.colour = "white", size = 3, hjust = -0.1, vjust = -0.5) + ggplot2::labs(x = "Longitude", y = "Latitude") my_map
\newpage
Temperature
# Create text objects for temp discussion temp_rank_text <- if(temp_text == "error"){ glue::glue("There are too many missing records to meaningfully summarize temperature data from this station through {month.name[match(my_mth$month_name, month.abb)]} {my_mth$year}.") } else if(!is.na(text_refs$temp_coldest) & text_refs$temp_coldest == 1){ glue::glue("{month.name[match(my_mth$month_name, month.abb)]} {my_mth$year} was the coldest {month.name[match(my_mth$month_name, month.abb)]} on record at {my_station$name}.") } else if(temp_text == "cool"){ glue::glue("{month.name[match(my_mth$month_name, month.abb)]} {my_mth$year} was the {number_contraction(my_mth_temp$dry_cool)} coolest {month.name[match(my_mth$month_name, month.abb)]} in the {text_refs$timespan}-year record at {my_station$name}.") } else if(temp_text == "avg"){ glue::glue("Temperatures were normal in {month.name[match(my_mth$month_name, month.abb)]} {my_mth$year}. This month was the {number_contraction(my_mth_temp$wet_warm)} warmest {month.name[match(my_mth$month_name, month.abb)]} in the {text_refs$timespan}-year record at {my_station$name}.") } else if(!is.na(text_refs$temp_warmest) & text_refs$temp_warmest == 1){ glue::glue("{month.name[match(my_mth$month_name, month.abb)]} {my_mth$year} was the hottest {month.name[match(my_mth$month_name, month.abb)]} on record at {my_station$name}.") } else if(temp_text == "warm"){ glue::glue("{month.name[match(my_mth$month_name, month.abb)]} {my_mth$year} was the {number_contraction(my_mth_temp$wet_warm)} warmest {month.name[match(my_mth$month_name, month.abb)]} in the {text_refs$timespan}-year record at {my_station$name}.") } temp_mth_text <- glue::glue("The average temperature for {month.name[match(my_mth$month_name, month.abb)]} was {round_with_zeros(my_mth_temp$value)}°F which is {round_with_zeros(abs(my_mth_temp$departure))}°F {ifelse(my_mth_temp$departure > 0, 'above', 'below')} the 30-year average ({round_with_zeros(my_mth_temp$norms_30yr)}°F).") temp_trend_text <- glue::glue("This summary suggests that temperatures are {ifelse(temp_lm > 0, 'increasing', 'decreasing')} {round_with_zeros(abs(temp_lm * 10))}°F per decade (Figure 1A).") # Text about monthly exceedences mths_so_far <- month_df[month_df$month_name == my_mth$month_name, ]$water_month temp_tmax_text <- if(length(tmax_exceed) == 0){ glue::glue("Every monthly average TMAX has been below normal this water year (Figure 1B).") } else if(length(tmax_exceed) == mths_so_far){ glue::glue("Every monthly average TMAX has been above normal this water year (Figure 1B).") } else( glue::glue("The monthly average TMAX has been {ifelse(length(tmax_exceed) >= (mths_so_far / 2), 'above', 'below')} normal most of the year and has exceeded the 30-year monthly average {length(tmax_exceed)} times ({glue_mths(tmax_exceed)}; Figure 1B).") ) temp_tmin_text <- if(length(tmin_exceed) == 0){ glue::glue("Every monthly average TMIN has been below normal this water year (Figure 1B).") } else if(length(tmin_exceed) == mths_so_far){ glue::glue("Every monthly average TMIN has been above normal this water year (Figure 1B).") } else( glue::glue("The monthly average TMIN has been {ifelse(length(tmin_exceed) >= (mths_so_far / 2), 'above', 'below')} normal most of the year and has exceeded the 30-year monthly average {length(tmin_exceed)} times ({glue_mths(tmin_exceed)}; Figure 1B).") ) if(temp_text == "error"){ temp_error_trend_text <- glue::glue("However, this summary suggests that temperatures are {ifelse(temp_lm > 0, 'increasing', 'decreasing')} {round_with_zeros(abs(temp_lm * 10))}°F per decade (Figure 1A).") }
# Text is printed in PDF if `temp_text %in% c("cool", "avg", "warm")` returns TRUE cat(temp_rank_text, temp_mth_text, temp_trend_text, temp_tmax_text, temp_tmin_text, sep = " ")
# Text for errors in monthly data set. Annual data set is good. cat(temp_rank_text, "Examine the monthly data set on ClimateAnalyzer.org for more information.", temp_error_trend_text, sep = " ")
# Text is printed if temp_text == FALSE cat("Error: Temperature data could not be evaluated. Examine tabular data sets on ClimateAnalyzer.org.")
#-- Average Annual Temps t1 <- annual_figure(x_var = annual_dat$year, y_var = annual_dat$temp_midpt, my_year = my_year, normal = annual_30yr$TEMP_midpt, reference_period = "1981-2010", area_color = "orange", line_color = "red3", my_title = "A. Average Annual Temperature", my_ylab = "Temperature (°F)") #-- Average WY TMAX/TMIN t2 <- monthly_figure(x_var = mth_depart$month_name, y_var = mth_depart$tmax_depart, year_group = mth_depart$year, y_intercept = 0, my_year = my_year, line_color = "red3", my_title = "B. Average Monthly TMAX", my_ylab = "Departure (°F)") + ggplot2::theme(axis.text.x = ggplot2::element_blank()) #-- TMAX/TMIN Departure t3 <- monthly_figure(x_var = mth_depart$month_name, y_var = mth_depart$tmin_depart, year_group = mth_depart$year, y_intercept = 0, my_year = my_year, line_color = "red3", my_title = "C. Average Monthly TMIN", my_ylab = "Departure (°F)") cowplot::plot_grid(t1, t2, t3, labels = c(), ncol = 1, align = "v", vjust = 1, scale = 1)
\newpage
Precipitation
# Create text objects for prcp discussion prcp_rank_text <- if(prcp_text == "error"){ glue::glue("There are too many missing records to meaningfully summarize precipitation data from this station for {month.name[match(my_mth$month_name, month.abb)]} {my_mth$year}.") } else if(prcp_text == "dry" & !is.na(text_refs$prcp_driest) & text_refs$prcp_driest == 1){ glue::glue("{month.name[match(my_mth$month_name, month.abb)]} {my_mth$year} was the driest {month.name[match(my_mth$month_name, month.abb)]} on record at {my_station$name}.") } else if(prcp_text == "dry"){ glue::glue("{month.name[match(my_mth$month_name, month.abb)]} {my_mth$year} was the {number_contraction(my_mth_prcp$dry_cool)} driest {month.name[match(my_mth$month_name, month.abb)]} in the {text_refs$timespan}-year record at {my_station$name}.") } else if(prcp_text == "avg"){ glue::glue("Precipitation was normal in {month.name[match(my_mth$month_name, month.abb)]} {my_mth$year}. This month was the {number_contraction(my_mth_prcp$dry_cool)} driest water year in the {text_refs$timespan}-year record at {my_station$name}.") } else if(prcp_text == "wet" & !is.na(text_refs$prcp_wettest) & text_refs$prcp_wettest == 1){ glue::glue("{month.name[match(my_mth$month_name, month.abb)]} {my_mth$year} was the wettest {month.name[match(my_mth$month_name, month.abb)]} on record at {my_station$name}.") } else if(prcp_text == "wet"){ glue::glue("{month.name[match(my_mth$month_name, month.abb)]} {my_mth$year} was the {number_contraction(my_mth_prcp$wet_warm)} wettest {month.name[match(my_mth$month_name, month.abb)]} in the {text_refs$timespan}-year record at {my_station$name}.") } prcp_annual_text <- glue::glue("Accumulated precipitation is {round_with_zeros(my_mth_prcp$prcp_accum)} inches so far this year which is {round_with_zeros(abs(my_mth_prcp$accum_depart))} inches {ifelse(my_mth_prcp$accum_depart < 0, 'below', 'above')} the 30-year average ({round_with_zeros(my_mth_prcp$prcp_accum_30yr)} inches).") prcp_trend_text <- glue::glue("This summary suggests that precipitation is {ifelse(prcp_lm > 0, 'increasing', 'decreasing')} at a rate of {round_with_zeros(abs(prcp_lm * 10))} inches per decade (Figure 2A).") prcp_condition_text <- if(prcp_text == "dry" & length(prcp_above) == 0){ glue::glue("Every month has received below average precipitation so far this year (Figures 2B). The accumulated precipitation is {ifelse(my_mth_prcp$accum_depart < 0, 'below', 'above')} the 30-year average up to this point (Figure 2C).") } else if(prcp_text == "dry"){ glue::glue("{glue_mths(prcp_above)} received above average precipitation (Figure 2B) but have not been enought to compensate for the {length(prcp_below)} months that were below average ({glue_mths(prcp_below)}; Figure 2C).") } else if(prcp_text == "avg"){ glue::glue("{glue_mths(prcp_above)} received above average precipitation and {glue_mths(prcp_below)} received below average precipitation (Figure 2B). {my_year} is {ifelse(my_mth_prcp$accum_depart < 0, 'below', 'above')} the 30-year average to this point (Figure 2C).") } else if(prcp_text == "wet" & length(prcp_below) == 0){ glue::glue("Every month has received above average precipitation so far this year (Figures 2B). This station has received {round_with_zeros(my_mth_prcp$accum_depart)} inches of precipitation so far this year (Figure 2C).") } else if(prcp_text == "wet"){ glue::glue("{glue_mths(prcp_above)} received above average precipitation (Figure 2B). The accumulated precipitation is {ifelse(my_mth_prcp$accum_depart < 0, 'below', 'above')} the 30-year average up to this point (Figure 2C).") } if(prcp_text == "error"){ prcp_error_trend_text <- glue::glue("However, this summary suggests that precipitation is {ifelse(prcp_lm > 0, 'increasing', 'decreasing')} by {round_with_zeros(abs(prcp_lm * 10))} inches per decade (Figure 1A).") }
cat( prcp_rank_text, prcp_annual_text, prcp_trend_text, prcp_condition_text, sep = " ")
# Text for errors in monthly and annual data sets. cat(prcp_rank_text, "Examine annual and monthly data sets on ClimateAnalyzer.org for more information.", prcp_error_trend_text, sep = " ")
# Text is printed if prcp_text == FALSE cat("Error: precipitation data could not be evaluated. Examine tabular data set on ClimateAnalyzer.org.")
#-- Water year precipitation totals p1 <- annual_figure(x_var = annual_dat$year, y_var = annual_dat$prcp_cy, my_year = my_year, normal = annual_30yr$PRCP, reference_period = "1981-2010", area_color = "green3", line_color = "darkgreen", my_title = "A. Precipitation Totals", my_ylab = "Precipitation (inches)") #-- Monthly precipitation totals p2 <- monthly_figure(x_var = mth_depart$month_name, y_var = mth_depart$prcp_pctavg, year_group = mth_depart$year, y_intercept = 100, my_year = my_year, line_color = "darkgreen", my_title = "B. Monthly Departure Realtive to the Historic Record", my_ylab = "Percent Average (%)") + ggplot2::theme(axis.text.x = ggplot2::element_blank()) #-- Water year accumulation p3 <- monthly_figure(x_var = mth_depart$month_name, y_var = mth_depart$prcp_accum_depart, year_group = mth_depart$year, y_intercept = 0, my_year = my_year, line_color = "darkgreen", my_title = "C. Accumulated Departure Realtive to the Historic Record", my_ylab = "Departure (inches)") # Using cowplot to arrange figures cowplot::plot_grid(p1, p2, p3, labels = c(""), ncol = 1, align = "v")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.