2 min read

Cumulative rainfall graph code

This code builds the cumulative rainfall graph in this post.

library(tidyverse)
library(ggthemes)
library(readr)
library(gsheet)

#1938-2017
old <- gsheet2tbl("https://docs.google.com/spreadsheets/d/1536nQ_xtBdZuccG2tOKRCdWvlJvhf7D1xFFMq3_g4rc/edit?usp=sharing")
old <- old[old$year>1940,]
d2018 <- read_table(file = "ftp://ftp.ncdc.noaa.gov/pub/data/gsod/2018/726400-14839-2018.op.gz")
d2018 <- d2018[,c("STN---","YEARMODA","MAX","MIN","PRCP")]
d2018 <- separate(d2018, "YEARMODA", sep = c(4,6),
                  into = c("year","month","day"))
names(d2018) <- c("station","year","month","day","MaxTemp","MinTemp",
                  "Precip")
d2018$MaxTemp <- as.numeric(str_replace(d2018$MaxTemp, "[*]", ""))
d2018$MinTemp <- as.numeric(str_replace(d2018$MinTemp, "[*]", ""))
d2018$Precip <- as.numeric(str_replace(d2018$Precip, "G|I", ""))
d2018$Precip[is.na(d2018$Precip)] <- 0
d2018$dayofyear <- 1:length(d2018$station)
d2018$cumprecip <- cumsum(d2018$Precip)
allweather <- rbind(old, d2018)
allweather <- allweather[allweather$MaxTemp!=9999.9 & allweather$MinTemp!=9999.9 &
                           allweather$Precip!=99.99,]
allweather <- allweather[as.Date(paste(allweather$year, allweather$month, allweather$day,
                                       sep = "-"))!=Sys.Date(),]
last.day <- as.Date(paste(last(allweather$year), last(allweather$month),
                          last(allweather$day), sep="/"))


# Build the ggplot graph
precipitation <- ggplot(allweather, aes(allweather$dayofyear, allweather$cumprecip, group=year)) +
  geom_step(col="grey", size=.2) +
  geom_step(data = allweather[allweather$year==2017,], aes(dayofyear, cumprecip, group=year),
            col="black", size = 1) +
  geom_step(data = allweather[allweather$year==2018,], aes(dayofyear, cumprecip, group=year),
            col="blue", size = 1) +
  ggtitle("Cumulative precipitation (in.), 1941-2018") + 
  geom_vline(xintercept = 31, colour = "wheat4", linetype=3, size=.2) +
  geom_vline(xintercept = 59, colour = "wheat4", linetype=3, size=.2) +
  geom_vline(xintercept = 90, colour = "wheat4", linetype=3, size=.2) +
  geom_vline(xintercept = 120, colour = "wheat4", linetype=3, size=.2) +
  geom_vline(xintercept = 151, colour = "wheat4", linetype=3, size=.2) +
  geom_vline(xintercept = 181, colour = "wheat4", linetype=3, size=.2) +
  geom_vline(xintercept = 212, colour = "wheat4", linetype=3, size=.2) +
  geom_vline(xintercept = 243, colour = "wheat4", linetype=3, size=.2) +
  geom_vline(xintercept = 273, colour = "wheat4", linetype=3, size=.2) +
  geom_vline(xintercept = 304, colour = "wheat4", linetype=3, size=.2) +
  geom_vline(xintercept = 334, colour = "wheat4", linetype=3, size=.2) +
  geom_vline(xintercept = 365, colour = "wheat4", linetype=3, size=.2) +
  theme(plot.background = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.ticks = element_blank(),
        #axis.text = element_blank(),  
        axis.title = element_blank()) +
  scale_x_continuous(expand = c(0, 0), 
                     breaks = c(15,45,75,105,135,165,195,228,258,288,320,350),
                     labels = c("Jan", "Feb", "March", "April",
                                "May", "June", "July", "Aug", "Sept",
                                "Oct", "Nov", "Dec")) +
  theme(plot.title=element_text(face="bold",hjust=.012,vjust=.8,colour="#3C3C3C",size=20)) +
  annotate("text", x=10, y=49, label = "General Mitchell Airport", size=4, fontface="bold", hjust=0) +
  annotate("segment", x=15, xend=5, y=10, yend=0, color="blue") +
  annotate("text", x=15, y=11, label = "2018", color="blue") +
  annotate("segment", x=340, xend=365, y=46, yend=44.3) +
  annotate("text", x=340, y=47, label = "2000 & 2008") +
  annotate("segment", x=153, xend=159, y=26, yend=22) +
  annotate("text", x=153, y=27, label = "2008") +
  annotate("segment", x=80, xend=90, y=15, yend=8, color="black") +
  annotate("text", x=80, y=16, label = "2017", color="black")
print(precipitation)