Code for downloading and plotting deviation in average precipitation for a given weather station. Using R and ggplot

precip_deviation_by_year.R

library('rnoaa')
library("tidyverse")
library("lubridate")
library("ggrepel")

token = 'GET YOUR API KEY at: http://www.ncdc.noaa.gov/cdo-web/token'
locs <- ncdc_locs(locationcategoryid='CITY', sortfield='name', sortorder='desc', token = token, limit = 800)
loc_data <- locs$data
dplyr::filter(loc_data, grepl(", PA",loc_data$name))

station_locs <- ncdc_stations(datasetid='GHCND', 
                              locationid='CITY:US420015', 
              token = token, limit = 500, sortfield = 'name')
station_data <- station_locs$data
dplyr::filter(station_data, grepl("PHILADELPHIA",station_data$name))

years = seq(1948,2018,1)
months = str_pad(seq(1:12),width = 2,side = "left", pad = "0")
years_dat <- list()
indx = 1
for(i in seq_along(years)){
  for(j in seq_along(months)){
    start_date = paste0(years[i],"-",months[j],"-01")
    end_date = as.character(ceiling_date(date(start_date), "month")-1)
    cat("looking up year",years[i],"and month",months[j],"\n")
    year_dat <- ncdc(datasetid='GHCND', 
                stationid='GHCND:USW00013739', datatypeid='PRCP', 
                startdate = start_date, 
                enddate = end_date, 
                limit=800,
                token = token)$data
    
    Sys.sleep(0.5) # avoid rate limits
    year_dat2 <- year_dat %>%
      mutate(total = cumsum(value),
             id = seq(1:n()),
             year = year(date),
             month = month(date))
    
    years_dat[[indx]] <- year_dat2
    indx = indx + 1
    } # end month loop
}

year_plot <- dplyr::bind_rows(years_dat)  %>% 
  arrange(date) %>% 
  group_by(year) %>% 
  mutate(year_id = seq(1:n()),
         month_label = month(date, label = TRUE),
         cm  = value/100,
         year_cumsum = cumsum(cm),
         year_total = sum(cm),
         day = day(date),
         year_label = ifelse(year == 2018 & day == 13 & month == 12,
                             "2018", "")) %>% 
  ungroup() %>% 
  group_by(month, day) %>% 
  mutate(day_cumsum_med = median(year_cumsum)) %>% 
  ungroup() %>% 
  mutate(over_under = year_cumsum - day_cumsum_med)

# write_csv(year_plot, "./R/rain_totals/PHL_1948_to_midDec_2018.csv")

ggplot(data = year_plot, aes(x = year_id, y = over_under,
                             group = year)) +
  geom_hline(yintercept = 0, color = "gray30", alpha = 0.8) +
  geom_line(color = ifelse(year_plot$year == 2018, "red", "gray45"),
            alpha = ifelse(year_plot$year == 2018, 0.9, 0.35),
            size  = ifelse(year_plot$year == 2018, 1, 0.5)) +
  geom_text_repel(aes(label = year_label), 
                  nudge_x = 15, nudge_y = 8) +
  theme_bw() +
  labs(y = "Precipitation (cm)",
       x = "Month",
       title = "Yearly Deviation from Median Precipitation: 1948 to 2018",
       subtitle = "Philadelphia International Airport") +
  scale_x_continuous(breaks = seq(1,360,30), # not exact
                     labels = unique(year_plot$month_label)) +
  scale_color_viridis_d(option = "D", direction = -1) +
  theme(
    text=element_text(size=16,  family="Trebuchet MS"),
    legend.position = "none",
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.border = element_blank()
  )

ggsave("./R/rain_totals/PHL_1948_to_midDec_2018.png", width = 10, height = 6)