scoyoc/climateAnalyzeR
Import Data from ClimateAnalyzer.org into R

Package index
Search the scoyoc/climateAnalyzeR package
Vignettes
Functions
45
Source code
16
Man pages
23
Home
/
GitHub
/
scoyoc/climateAnalyzeR
/

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)

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

# 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)))) %>%
  dplyr::mutate(year = my_year) %>%
  dplyr::left_join(tibble::tibble(annual_dat[, c("year", "prcp_cy", 
                                                 "temp_midpt")]) %>%
                     dplyr::filter(year == my_year) %>%
                     dplyr::rename('ann_temp' = 'temp_midpt'))

prcp_error <- if(is.na(error_eval$prcp_cy) & error_eval$prcp >= 2){
  "both"
  } else if(is.na(error_eval$prcp_cy) & error_eval$prcp < 2){
    "annual"
    } else if(!is.na(error_eval$prcp_cy) & error_eval$prcp >= 2){
      "month"
      } else("no error")

temp_error <- if(is.na(error_eval$ann_temp) & error_eval$temp_midpt >= 2){
  "both"
  } else if(is.na(error_eval$ann_temp) & error_eval$temp_midpt < 2){
    "annual"
    } else if(!is.na(error_eval$ann_temp) & error_eval$temp_midpt >= 2){
      "month"
      } else("no error")

# Evaluate data for conditional text
# Calculate percentiles
prcp_quant <- quantile(annual_dat$prcp_depart, c(0.25, 0.75), na.rm = TRUE)
temp_quant <- quantile(annual_dat$temp_midpt, c(0.25, 0.75), na.rm = TRUE)

# Assign condition
prcp_text <- if(prcp_error %in% c("both", "annual", "month")){
  "error"
  } else if(text_refs$prcp_depart <= prcp_quant[1]){
    "dry"
    } else if(text_refs$prcp_depart > prcp_quant[1] & 
               text_refs$prcp_depart < prcp_quant[2]){
      "avg"
      } else if(text_refs$prcp_depart >= prcp_quant[2]){
        "wet"
        } else(FALSE)

temp_text <- if(temp_error %in% c("both", "annual", "month")){
  'error'
  } else if(text_refs$temp_midpt <= temp_quant[1]){
    "cool"
    } else if(text_refs$temp_midpt > temp_quant[1] & 
               text_refs$temp_midpt < temp_quant[2]){
      "avg"
      } else if(text_refs$temp_midpt >= 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. 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_error == "annual"){
  glue::glue("There are too many missing records to meaningfully summarize the annual data from this station for {my_year}.")
  } else if(!is.na(text_refs$temp_coldest) & text_refs$temp_coldest == 1){
    glue::glue("{my_year} was the coldest year on record at {my_station$name}.")
    } else if(temp_text == "cool"){
      glue::glue("{my_year} was the {number_contraction(text_refs$temp_coldest)} coolest year in the {text_refs$timespan}-year record for {my_station$name}.")
      } else if(temp_text == "avg"){
        glue::glue("Temperatures were normal for {my_year}. This year was the {number_contraction(text_refs$temp_warmest)} warmest year in the {text_refs$timespan}-year record for {my_station$name}.")
        } else if(!is.na(text_refs$temp_warmest) & text_refs$temp_warmest == 1){
          glue::glue("{my_year} was the warmest year on record for {my_station$name}.")
          } else if(temp_text == "warm"){
            glue::glue("{my_year} was the {number_contraction(text_refs$temp_warmest)} warmest year in the {text_refs$timespan}-year record for {my_station$name}.")
            }
temp_annual_text <- glue::glue("The average annual temperature was {round_with_zeros(text_refs$temp_midpt)}°F which is {round_with_zeros(abs(text_refs$temp_depart))}°F {ifelse(text_refs$temp_depart > 0, 'above', 'below')} the 30-year average ({round_with_zeros(annual_30yr$TEMP_midpt)}°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).")

temp_tmax_text <- if(length(tmax_exceed) == 0){
    glue::glue("Every monthly average TMAX was below normal in {my_year} (Figure 1B).")
    } else if(length(tmax_exceed) == 6){
      glue::glue("Montly average TMAX was {ifelse(temp_text == 'cool', 'below', 'above')} normal for half of {my_year} (Figure 1B).")
      } else if(length(tmax_exceed) == 12){
        glue::glue("Every monthly average TMAX was above normal in {my_year} (Figure 1B).")
        } else(
          glue::glue("The monthly average TMAX was {ifelse(length(tmax_exceed) >= 6, 'above', 'below')} normal most of the year and 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 was below normal in {my_year} (Figure 1B).")
    } else if(length(tmin_exceed) == 6){
      glue::glue("Montly average TMIN was {ifelse(temp_text == 'cool', 'below', 'above')} normal for half of {my_year} (Figure 1B).")
      } else if(length(tmin_exceed) == 12){
        glue::glue("Every monthly average TMIN was above normal in {my_year} (Figure 1B).")
          } else(
            glue::glue("The monthly average TMIN was {ifelse(length(tmin_exceed) >= 6, 'above', 'below')} normal most of the year and exceeded the 30-year monthly average {length(tmin_exceed)} times ({glue_mths(tmin_exceed)}; Figure 1B).")
            )

# Text is printed in PDF if `temp_text %in% c("cool", "avg", "warm")` returns TRUE
cat(
  temp_rank_text, 
  temp_annual_text, 
  temp_trend_text, 
  temp_tmax_text, 
  temp_tmin_text, 
  sep = " ")

# Text for errors in temperature data sets.
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 for errors in annual and monthly data
if(temp_error == "both"){
  cat(
    glue::glue("There are too many missing records to meaningfully summarize temperature data from this station for {my_year}."), 
    "Examine annual and monthly data sets on ClimateAnalyzer.org for more information.", 
    temp_error_trend_text,
    sep = " ")
}

# Text for errors in annual data set. Monthly data set is good.
if(temp_error == "annual"){
  cat(
    temp_rank_text,  
    "Examine the annual data set on ClimateAnalyzer.org for more information.", 
    temp_error_trend_text,
    temp_tmax_text,
    temp_tmin_text,
    sep = " ")
}

# Text for errors in monthly data set. Annual data set is good.
if(temp_error == "month"){
  cat(
    glue::glue("There are too many missing records to meaningfully summarize the montly temperature data from this station for {my_year}."),
    "Examine the monthly data set on ClimateAnalyzer.org for more information.", 
    temp_rank_text, 
    temp_annual_text,
    temp_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_error == "annual"){
  glue::glue("There are too many missing records to meaningfully summarize the annual data from this station for {my_year}.")
  } else if(text_refs$prcp_driest == 1){
    glue::glue("{my_year} was the driest year on record for {my_station$name}.")
    } else if(prcp_text == "dry"){
      glue::glue("{my_year} was the {number_contraction(text_refs$prcp_driest)} driest on record for {station_name}.")
      } else if(prcp_text == "avg"){
        glue::glue("{station_name} received an average amount of precipitation during {my_year} and was the {number_contraction(text_refs$prcp_driest)} driest year on record.")
        } else if(text_refs$prcp_wettest == 1){
          glue::glue("{my_year} was the wettest year on record for {my_station$name}.")
          } else(
            glue::glue("{my_year} was the {number_contraction(text_refs$prcp_wettest)} wettest on record for {station_name}.")
            )

prcp_annual_text <- glue::glue("Total accumulated precipitation was {text_refs$prcp_cy} inches and was {round_with_zeros(abs(text_refs$prcp_cy - annual_30yr$PRCP))} inches {ifelse(text_refs$prcp_cy - annual_30yr$PRCP < 0, 'below', 'above')} the 30-year average ({round_with_zeros(annual_30yr$PRCP)} 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 received below average precipitation in {my_year} (Figures 2B and  2C).")
  } else if(prcp_text == "dry"){
    glue::glue("{glue_mths(prcp_above)} received above average precipitation (Figure 2B) but were not enough 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 (Figure 2B) and {glue_mths(prcp_below)} received below average precipitation. {my_year} finished {ifelse(text_refs$prcp_cy < 0, 'below', 'above')}  the 30-year average (Figure 2C).")
      } else if(prcp_text == "wet" & length(prcp_below) == 0){
        glue::glue("Every month received above average precipitation (Figures 2B 2C).")
        } else if(prcp_text == "wet"){
          glue::glue("{glue_mths(prcp_above)} received above average precipitation (Figure 2B) and {my_year} finished {ifelse(text_refs$prcp_cy < 0, 'below', 'above')} the 30-year average (Figure 2C).")
          }

cat(
  prcp_rank_text, 
  prcp_annual_text, 
  prcp_trend_text, 
  prcp_condition_text, 
  sep = " ")

# Text objects for errors in prcp data sets.
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).")

# Text for errors in monthly and annual data sets.
if(prcp_error == "both"){
  cat(
  glue::glue("There are too many missing records to meaningfully summarize precipitation data from this station for {my_year}."),
  "Examine annual and monthly data sets on ClimateAnalyzer.org for more information.",
 prcp_error_trend_text, 
 sep = " ")
  }

# Text for errors in annual data set. Monthly data set is good.
if(prcp_error == "annual"){
  cat(
  prcp_rank_text,
  "Examine the annual data set on ClimateAnalyzer.org for more information.",
  prcp_error_trend_text,
  prcp_condition_text, 
  sep = " ")
}

# Text for errors in monthly data set. Annual data set is good.
if(prcp_error == "month"){
  cat(
    glue::glue("There are too many missing records to meaningfully summarize the montly precipitation data from this station for {my_year}."),
  "Examine the montly data set on ClimateAnalyzer.org for more information.",
  prcp_rank_text, 
  prcp_annual_text, 
  prcp_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 Relative 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 Relative 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

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
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.

Close