Adapted from a blog post by Sam Tyner

library(dplyr)
library(tidyr)
library(rnoaa)
library(ggplot2)
library(ggrepel)
library(patchwork)
theme_set(theme_minimal() + theme(text = element_text(size = 16)))

Getting the Data from NOAA

Climate data is available via a web service API from NOAA.

The rnoaa package provides an interface to this API. A tutorial is available from rOpenSci.

This code would be used to find the most suitable Iowa City station:

options("<key>") ## may not be needed
stat <- ghcnd_stations()
ic_stat <- filter(stat, grepl("IOWA CITY", name) & element == "SNOW")
select(ic_stat, name, id, ends_with("year"))

The best station seems to be

ic_id <- "USC00134101"

The current data for this station can then be obtained from NOAA:

ic_data_raw <- ghcnd(ic_id, refresh = TRUE)
## goes into a cache file as csv
## /home/luke/.cache/rnoaa/ghcnd/USC00134101.dly
write.csv(ic_data_raw, row.names = FALSE, file = "ic_noaa.csv")
## gzip and move to data directory

I have done these steps and stored the data locally; I will update the data as the course progresses.

Download the local data copy:

if (! file.exists("ic_noaa.csv.gz"))
    download.file("http://www.stat.uiowa.edu/~luke/data/ic_noaa.csv.gz",
                  "ic_noaa.csv.gz")
## make sure data is up to date
ic_data_raw <- as_tibble(read.csv("ic_noaa.csv.gz", head = TRUE))

Processing the Data

The raw data:

ic_data_raw
## # A tibble: 10,133 × 128
##    id        year month element VALUE1 MFLAG1 QFLAG1 SFLAG1 VALUE2 MFLAG2 QFLAG2
##    <chr>    <int> <int> <chr>    <int> <chr>  <chr>  <chr>   <int> <chr>  <chr> 
##  1 USC0013…  1893     1 TMAX       -17 " "    " "    6         -28 " "    " "   
##  2 USC0013…  1893     1 TMIN       -67 " "    " "    6        -161 " "    " "   
##  3 USC0013…  1893     1 PRCP         0 "P"    " "    6           0 "P"    " "   
##  4 USC0013…  1893     1 SNOW         0 " "    " "    6           0 " "    " "   
##  5 USC0013…  1893     2 TMAX        11 " "    " "    6         -56 " "    " "   
##  6 USC0013…  1893     2 TMIN      -233 " "    " "    6        -206 " "    " "   
##  7 USC0013…  1893     2 PRCP         0 "P"    " "    6         102 " "    " "   
##  8 USC0013…  1893     2 SNOW         0 " "    " "    6          20 " "    " "   
##  9 USC0013…  1893     3 TMAX        56 " "    " "    6          44 " "    " "   
## 10 USC0013…  1893     3 TMIN       -78 " "    " "    6         -67 " "    " "   
## # ℹ 10,123 more rows
## # ℹ 117 more variables: SFLAG2 <chr>, VALUE3 <int>, MFLAG3 <chr>, QFLAG3 <chr>,
## #   SFLAG3 <chr>, VALUE4 <int>, MFLAG4 <chr>, QFLAG4 <chr>, SFLAG4 <chr>,
## #   VALUE5 <int>, MFLAG5 <chr>, QFLAG5 <chr>, SFLAG5 <chr>, VALUE6 <int>,
## #   MFLAG6 <chr>, QFLAG6 <chr>, SFLAG6 <chr>, VALUE7 <int>, MFLAG7 <chr>,
## #   QFLAG7 <chr>, SFLAG7 <chr>, VALUE8 <int>, MFLAG8 <chr>, QFLAG8 <chr>,
## #   SFLAG8 <chr>, VALUE9 <int>, MFLAG9 <chr>, QFLAG9 <chr>, SFLAG9 <chr>, …

The documentation for the data is available at

https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/readme.txt

The available elements, or measurements, are:

unique(ic_data_raw$element)
##  [1] "TMAX" "TMIN" "PRCP" "SNOW" "WT16" "WT01" "TOBS" "SNWD" "WT18" "WT03"
## [11] "WT04" "WT14" "EVAP" "WDMV" "WT05" "DAEV" "MDEV" "DAWM" "MDWM" "DAPR"
## [21] "MDPR" "MNPN" "MXPN" "WESD"

The primary ones are:

PRCP = Precipitation (tenths of mm)
SNOW = Snowfall (mm)
SNWD = Snow depth (mm)
TMAX = Maximum temperature (tenths of degrees C)
TMIN = Minimum temperature (tenths of degrees C)

The values for day \(k\) are stored in the variable VALUE\(k\). A selection:

ic_feb2019 <- filter(ic_data_raw, year == 2019, month == 2)
select(ic_feb2019, element, VALUE1, VALUE28, VALUE31)
## # A tibble: 6 × 4
##   element VALUE1 VALUE28 VALUE31
##   <chr>    <int>   <int>   <int>
## 1 TMAX      -167     -72      NA
## 2 TMIN      -344    -139      NA
## 3 TOBS      -194    -128      NA
## 4 PRCP        18       0      NA
## 5 SNOW        51       0      NA
## 6 SNWD       330      51      NA

We need to:

Start by selecting the columns we need:

ic_data <- select(ic_data_raw, year, month, element, starts_with("VALUE"))
ic_data
## # A tibble: 10,133 × 34
##     year month element VALUE1 VALUE2 VALUE3 VALUE4 VALUE5 VALUE6 VALUE7 VALUE8
##    <int> <int> <chr>    <int>  <int>  <int>  <int>  <int>  <int>  <int>  <int>
##  1  1893     1 TMAX       -17    -28   -150    -44    -17   -106    -56    -67
##  2  1893     1 TMIN       -67   -161   -233   -156   -156   -206   -111   -211
##  3  1893     1 PRCP         0      0      0     64      0     38      0      0
##  4  1893     1 SNOW         0      0      0     64      0     38      0      0
##  5  1893     2 TMAX        11    -56   -106   -122     33     28   -161    -94
##  6  1893     2 TMIN      -233   -206   -217   -267   -122   -211   -256   -278
##  7  1893     2 PRCP         0    102      0      0      0    102      0      0
##  8  1893     2 SNOW         0     20      0      0      0    102      0      0
##  9  1893     3 TMAX        56     44     17    -28     17     61     83     72
## 10  1893     3 TMIN       -78    -67   -122   -156   -122   -106    -28     17
## # ℹ 10,123 more rows
## # ℹ 23 more variables: VALUE9 <int>, VALUE10 <int>, VALUE11 <int>,
## #   VALUE12 <int>, VALUE13 <int>, VALUE14 <int>, VALUE15 <int>, VALUE16 <int>,
## #   VALUE17 <int>, VALUE18 <int>, VALUE19 <int>, VALUE20 <int>, VALUE21 <int>,
## #   VALUE22 <int>, VALUE23 <int>, VALUE24 <int>, VALUE25 <int>, VALUE26 <int>,
## #   VALUE27 <int>, VALUE28 <int>, VALUE29 <int>, VALUE30 <int>, VALUE31 <int>

Reshape from (very) wide to (too) long:

ic_data <- pivot_longer(ic_data, VALUE1 : VALUE31,
                        names_to = "day", values_to = "value")
ic_data
## # A tibble: 314,123 × 5
##     year month element day     value
##    <int> <int> <chr>   <chr>   <int>
##  1  1893     1 TMAX    VALUE1    -17
##  2  1893     1 TMAX    VALUE2    -28
##  3  1893     1 TMAX    VALUE3   -150
##  4  1893     1 TMAX    VALUE4    -44
##  5  1893     1 TMAX    VALUE5    -17
##  6  1893     1 TMAX    VALUE6   -106
##  7  1893     1 TMAX    VALUE7    -56
##  8  1893     1 TMAX    VALUE8    -67
##  9  1893     1 TMAX    VALUE9     11
## 10  1893     1 TMAX    VALUE10  -150
## # ℹ 314,113 more rows

Extract the day as a number:

ic_data <- mutate(ic_data, day = readr::parse_number(day))
ic_data
## # A tibble: 314,123 × 5
##     year month element   day value
##    <int> <int> <chr>   <dbl> <int>
##  1  1893     1 TMAX        1   -17
##  2  1893     1 TMAX        2   -28
##  3  1893     1 TMAX        3  -150
##  4  1893     1 TMAX        4   -44
##  5  1893     1 TMAX        5   -17
##  6  1893     1 TMAX        6  -106
##  7  1893     1 TMAX        7   -56
##  8  1893     1 TMAX        8   -67
##  9  1893     1 TMAX        9    11
## 10  1893     1 TMAX       10  -150
## # ℹ 314,113 more rows

Reshape from too long to tidy with one row per day, keeping only the primary variables:

corevars <- c("TMAX", "TMIN", "PRCP", "SNOW", "SNWD")
ic_data <- filter(ic_data, element %in% corevars)
ic_data <- pivot_wider(ic_data, names_from = "element", values_from = "value")
ic_data
## # A tibble: 47,120 × 8
##     year month   day  TMAX  TMIN  PRCP  SNOW  SNWD
##    <int> <int> <dbl> <int> <int> <int> <int> <int>
##  1  1893     1     1   -17   -67     0     0    NA
##  2  1893     1     2   -28  -161     0     0    NA
##  3  1893     1     3  -150  -233     0     0    NA
##  4  1893     1     4   -44  -156    64    64    NA
##  5  1893     1     5   -17  -156     0     0    NA
##  6  1893     1     6  -106  -206    38    38    NA
##  7  1893     1     7   -56  -111     0     0    NA
##  8  1893     1     8   -67  -211     0     0    NA
##  9  1893     1     9    11  -144     0     0    NA
## 10  1893     1    10  -150  -250     0     0    NA
## # ℹ 47,110 more rows

Add a date variable for plotting and to help get rid of bogus days:

ic_data <- mutate(ic_data, date = lubridate::make_date(year, month, day))
filter(ic_data, year == 2019, month == 2, day >= 27)
## # A tibble: 5 × 9
##    year month   day  TMAX  TMIN  PRCP  SNOW  SNWD date      
##   <int> <int> <dbl> <int> <int> <int> <int> <int> <date>    
## 1  2019     2    27   -44  -128     0     0    51 2019-02-27
## 2  2019     2    28   -72  -139     0     0    51 2019-02-28
## 3  2019     2    29    NA    NA    NA    NA    NA NA        
## 4  2019     2    30    NA    NA    NA    NA    NA NA        
## 5  2019     2    31    NA    NA    NA    NA    NA NA

Get rid of the bogus days:

ic_data <- filter(ic_data, ! is.na(date))

Make units more standard (American):

mm2in <- function(x) x / 25.4
C2F <- function(x) 32 + 1.8 * x
ic_data <- transmute(ic_data, year, month, day, date,
                     PRCP = mm2in(PRCP / 10),
                     SNOW = mm2in(SNOW),
                     SNWD = mm2in(SNWD),
                     TMIN = C2F(TMIN / 10),
                     TMAX = C2F(TMAX / 10))

The Winters of 2018/2019 and 2020/2021

Daily Snowfall

Extract the data for days between November 1 and March 31:

w18 <- filter(ic_data, date >= "2018-11-01" & date <= "2019-03-31")
w20 <- filter(ic_data, date >= "2020-11-01" & date <= "2021-03-31")

A first set of plots:

p1 <- ggplot(w18, aes(x = date)) + geom_col(aes(y = SNOW))
p2 <- ggplot(w18, aes(x = date)) + geom_line(aes(y = SNWD))
p1 + p2

For daily snowfall it is probably reasonable to replace missing values by zeros.

w18 <- mutate(w18, SNOW = ifelse(is.na(SNOW), 0, SNOW))
w20 <- mutate(w20, SNOW = ifelse(is.na(SNOW), 0, SNOW))

For snow depth it may make sense to fill forward, replacing misusing values by the previous non-missing value. Using the fill function from tidyr is one way to do this.

w18 <- fill(w18, SNWD)
w20 <- fill(w20, SNWD)

The plot of snow depth looks a bit better with this change; the axis range for the plot of daily snowfall is also affected. Also, the warnings are eliminated.

p1 <- ggplot(w18, aes(x = date)) + geom_col(aes(y = SNOW))
p2 <- ggplot(w18, aes(x = date)) + geom_line(aes(y = SNWD))
(p1 + labs(title = "2018/19")) + p2

(p1 %+% w20 + labs(title = "2020/21")) + (p2 %+% w20)

Often snowfall is visualized in terms of the cumulative snowfall for the season up to each date:

p <- ggplot(w18, aes(x = date)) + geom_line(aes(y = cumsum(SNOW)))
p + labs(title = "2018/19")

p %+% w20 + labs(title = "2020/21")

Comparing to Long Term Averages

Compute the average snowfall for each day of the year and merge the results for the days in w18 into w18 with a left_join:

ic_avg <- group_by(ic_data, month, day) %>%
    summarize(av_snow = mean(SNOW, na.rm = TRUE),
              av_snwd = mean(SNWD, na.rm = TRUE)) %>%
    ungroup()
w18 <- left_join(w18, ic_avg, c("month", "day"))
w20 <- left_join(w20, ic_avg, c("month", "day"))

Daily snowfall and snow cover for the 2018/19 winter along with averages over all years in the record:

p1 <- ggplot(w18, aes(x = date)) +
    geom_col(aes(y = SNOW)) +
    geom_line(aes(y = av_snow), color = "blue")
p2 <- ggplot(w18, aes(x = date)) +
    geom_line(aes(y = SNWD)) +
    geom_line(aes(y = av_snwd), color = "blue")
(p1 + labs(title = "2018/19")) + p2

(p1 %+% w20 + labs(title = "2020/21")) + (p2 %+% w20)

Cumulative snowfall for the 2018/19 winter along with average cumulative snowfall over all years in the record:

p <- ggplot(w18, aes(x = date)) +
    geom_line(aes(y = cumsum(SNOW))) +
    geom_line(aes(y = cumsum(av_snow)), color = "blue")
p + labs(title = "2018/19")

p %+% w20 + labs(title = "2020/21")

A Combined Graph

w18day <- transmute(w18, date, snow = SNOW,
                    base = "2018/19",
                    type = "Daily Snowfall")
nrmday <- transmute(w18, date, snow = av_snow,
                    base = "Normal",
                    type = "Daily Snowfall")
w18cum <- transmute(w18day, date, snow = cumsum(snow),
                    base = "2018/19",
                    type = "Cumulative Snowfall")
nrmcum <- transmute(nrmday, date, snow = cumsum(snow),
                    base = "Normal",
                    type = "Cumulative Snowfall")
ggplot(w18, aes(x = date, y = snow, color = base)) +
    geom_line(data = w18cum) +
    geom_line(data = nrmcum) +
    geom_segment(aes(x = date, xend = date, y = 0, yend = snow),
                 size = 1, data = w18day) +
    geom_line(data = nrmday) +
    facet_wrap(~ type, ncol = 1, scales = "free_y") +
    xlab(element_blank()) + ylab("Snowfall (inches)") +
    ggtitle("2018/19 Snowfall in Iowa City") +
    scale_color_manual(values = c("2018/19" = "navy",
                                  "Normal" = "forestgreen")) +
    guides(color = guide_legend(override.aes = list(size = 1.5))) +
    theme_bw() +
    theme(legend.position = "bottom", legend.title = element_blank(),
          plot.title = element_text(face = "bold", hjust = .5),
          aspect.ratio = 0.4)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Total Snow Fall by Year and Month

It is helpful to define a variable that shifts the year to include a full winter, not the end of one and the start of the next.

What month should be start with? Here are the counts for months with recorded snowfall:

count(filter(ic_data, SNOW > 0), month)
## # A tibble: 9 × 2
##   month     n
##   <int> <int>
## 1     1   559
## 2     2   456
## 3     3   300
## 4     4    73
## 5     5     1
## 6     7     1
## 7    10    20
## 8    11   136
## 9    12   477

The July snowfall seems odd:

filter(ic_data, SNOW > 0, month == 7)
## # A tibble: 1 × 9
##    year month   day date        PRCP  SNOW  SNWD  TMIN  TMAX
##   <int> <int> <dbl> <date>     <dbl> <dbl> <dbl> <dbl> <dbl>
## 1  1896     7    26 1896-07-26 0.110  2.91    NA    68  91.9

This is probably a recording error.

For October, the years with recorded recorded snowfall are shown by

count(filter(ic_data, SNOW > 0, month == 10), year)
## # A tibble: 14 × 2
##     year     n
##    <int> <int>
##  1  1898     1
##  2  1906     1
##  3  1913     2
##  4  1917     3
##  5  1923     1
##  6  1925     1
##  7  1929     1
##  8  1930     1
##  9  1967     1
## 10  1972     1
## 11  1980     2
## 12  1997     2
## 13  2019     2
## 14  2020     1

It seems reasonable to start the winter year in October:

ic_data <- mutate(ic_data, wyear = ifelse(month >= 10, year, year - 1))

The yearly total snowfall values:

tot <- group_by(ic_data, wyear) %>%
    summarize(tsnow = sum(SNOW, na.rm = TRUE))
tot18 <- filter(tot, wyear == 2018)
ggplot(tot, aes(x = wyear, y = tsnow)) +
    geom_line() +
    geom_point(data = tot18, color = "red") +
    geom_text_repel(data = tot18, label = "2018/19", color = "red")

So 2018/19 is not a record, but still pretty high:

tot100 <- mutate(filter(tot, wyear >= 1918), rank = min_rank(desc(tsnow)))
filter(tot100, wyear == 2018)
## # A tibble: 1 × 3
##   wyear tsnow  rank
##   <dbl> <dbl> <int>
## 1  2018  40.8    11

What are the months with the highest average total snowfall?

mtot <- group_by(ic_data, year, wyear, month) %>%
    summarize(msnow = sum(SNOW, na.rm = TRUE)) %>%
    ungroup()
## `summarise()` has grouped output by 'year', 'wyear'. You can override using the
## `.groups` argument.
group_by(mtot, month) %>% summarize(av_msnow = mean(msnow))
## # A tibble: 12 × 2
##    month av_msnow
##    <int>    <dbl>
##  1     1  7.41   
##  2     2  6.54   
##  3     3  4.22   
##  4     4  1.03   
##  5     5  0.00153
##  6     6  0      
##  7     7  0.0228 
##  8     8  0      
##  9     9  0      
## 10    10  0.256  
## 11    11  1.54   
## 12    12  6.27

For 2018/19:

filter(mtot, wyear == 2018)
## # A tibble: 12 × 4
##     year wyear month  msnow
##    <int> <dbl> <int>  <dbl>
##  1  2018  2018    10  0    
##  2  2018  2018    11  8.98 
##  3  2018  2018    12  0    
##  4  2019  2018     1 17.1  
##  5  2019  2018     2 14.2  
##  6  2019  2018     3  0.512
##  7  2019  2018     4  0    
##  8  2019  2018     5  0    
##  9  2019  2018     6  0    
## 10  2019  2018     7  0    
## 11  2019  2018     8  0    
## 12  2019  2018     9  0

For 2020/21:

filter(mtot, wyear == 2020)
## # A tibble: 8 × 4
##    year wyear month msnow
##   <int> <dbl> <int> <dbl>
## 1  2020  2020    10  1.50
## 2  2020  2020    11  0   
## 3  2020  2020    12  9.96
## 4  2021  2020     1 14.0 
## 5  2021  2020     2  9.69
## 6  2021  2020     3  0   
## 7  2021  2020     4  0   
## 8  2021  2020     5  0

A Data Processing Functions

If you want to look at data for another city you would need to go through the same steps starting with a new raw data frame.

Defining a function that encapsulates the steps makes this easy:

processData <- function(data_raw) {
    ## select just the columns we need:
    data <- select(data_raw, year, month, element, starts_with("VALUE"))

    ## reshape from (very) wide to (too) long
    data <- pivot_longer(data, VALUE1 : VALUE31,
                         names_to = "day", values_to = "value")
    ## extract the day as a number:
    data <- mutate(data, day = readr::parse_number(day))

    ## reshape from too long to tidy with one row per day, keeping
    ## only the primary variables:
    corevars <- c("TMAX", "TMIN", "PRCP", "SNOW", "SNWD")
    data <- filter(data, element %in% corevars)
    data <- pivot_wider(data, names_from = "element", values_from = "value")

    ## add a 'date' variable and get rid of bogus days
    data <- mutate(data, date = lubridate::make_date(year, month, day))
    data <- filter(data, ! is.na(date))

    ## make units more standard (American)
    mm2in <- function(x) x / 25.4
    C2F <- function(x) 32 + 1.8 * x
    data <- transmute(data, year, month, day, date,
                      PRCP = mm2in(PRCP / 10),
                      SNOW = mm2in(SNOW),
                      SNWD = mm2in(SNWD),
                      TMIN = C2F(TMIN / 10),
                      TMAX = C2F(TMAX / 10))

    ## add winter year variable 'wyear' starting in October
    data <- mutate(data, wyear = ifelse(month >= 10, year, year - 1))

    data
}

A quick sanity check:

stopifnot(identical(ic_data, processData(ic_data_raw)))

Often it is a good idea to include a sanity check like this in a code chunk marked as include = FALSE.

Other Explorations

pivot_longer(w18, TMIN | TMAX, names_to = "which", values_to = "value") %>%
    ggplot(aes(x = date)) +
    geom_line(aes(y = value, color = which))

wd18 <- filter(ic_data, date >= "2018-10-01")
avs <- group_by(ic_data, month, day) %>%
    summarize(avsnow = mean(SNOW, na.rm = TRUE),
              avsnwd = mean(SNWD, na.rm = TRUE),
              avtmin = mean(TMIN, na.rm = TRUE),
              avtmax = mean(TMAX, na.rm = TRUE)) %>%
    ungroup()
## `summarise()` has grouped output by 'month'. You can override using the
## `.groups` argument.
wd18 <- left_join(wd18, avs)
## Joining with `by = join_by(month, day)`

ggplot(filter(w18, ! is.na(SNWD))) +
    geom_line(aes(x = date, y = SNWD)) +
    geom_line(aes(x = date, y = av_snwd), col = "blue")

plot(SNOW ~ date, data = wd18, type = "h")
lines(wd18$date, wd18$avsnow, col = "red")
lines(smooth.spline(wd18$date, wd18$avsnow), col = "blue")


plot(SNWD ~ date, data = fill(wd18, SNWD), type = "l")
lines(wd18$date, wd18$avsnwd, col = "red")


plot(cumsum(ifelse(is.na(SNOW), 0, SNOW)) ~ date, data = wd18, type = "l")
lines(wd18$date, cumsum(wd18$avsnow), col = "red")

thm <- theme_bw() +
    theme(panel.grid.major.x = element_blank(),
          panel.grid.minor = element_blank())

ggplot(wd18) +
    geom_col(aes(x = date, y = SNOW)) +
    geom_line(aes(x = date, y = avsnow), col = "red") +
    geom_smooth(aes(x = date, y = avsnow), se = FALSE) +
    thm
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 25 rows containing missing values (`position_stack()`).


ggplot(fill(wd18, SNWD)) +
    geom_line(aes(x = date, y = SNWD)) +
    geom_line(aes(x = date, y = avsnwd), col = "red") +
    geom_smooth(aes(x = date, y = avsnwd), se = FALSE) +
    thm
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'


ggplot(wd18) +
    geom_line(aes(x = date, y = cumsum(ifelse(is.na(SNOW), 0, SNOW)))) +
    geom_line(aes(x = date, y = cumsum(wd18$avsnow)), col = "red") +
    thm
## Warning: Use of `wd18$avsnow` is discouraged.
## ℹ Use `avsnow` instead.

ggplot(filter(ic_data, year == 2019), aes(x = date)) +
    geom_line(aes(y = TMAX, color = "TMAX")) +
    geom_line(aes(y = TMIN, color = "TMIN")) +
    thm

d13 <- filter(ic_data, wyear >= 2013, month <= 4 | month >= 10)
d13 <- mutate(d13, wday = as.numeric(date - lubridate::make_date(wyear, 10, 1)),
              SNOW = ifelse(is.na(SNOW), 0, SNOW))
ggplot(d13) +
    geom_col(aes(x = wday, y = SNOW), width = 3, position = "identity") +
    facet_wrap(~ wyear) + thm

ggplot(d13) + geom_line(aes(x = wday, y = SNOW)) + facet_wrap(~ wyear) + thm

d13 <- group_by(d13, wyear) %>%
    mutate(snow_total = cumsum(SNOW)) %>%
    ungroup()
ggplot(d13) +
    geom_line(aes(x = wday, y = snow_total, color = factor(wyear))) + thm

ggplot(d13) +
    geom_line(aes(x = wday, y = snow_total)) +
    facet_wrap(~ wyear) + thm

d13 <- fill(d13, SNWD)
ggplot(d13) + geom_line(aes(x = wday, y = SNWD)) + facet_wrap(~ wyear) + thm

w <- filter(ic_data, wyear >= 2013)
w <- mutate(w, wday = date - lubridate::make_date(wyear, 10, 1))
w <- filter(w, wday <= 180)
w <- mutate(w, wday = lubridate::make_date(2018, 10, 1) + wday)
ggplot(w) +
    geom_line(aes(x = wday, y = TMAX), color = "blue") +
    geom_line(aes(x = wday, y = TMIN), color = "red") +
    facet_wrap(~ wyear)

LS0tCnRpdGxlOiAiU25vd2ZhbGwgaW4gSW93YSBDaXR5IgpvdXRwdXQ6CiAgaHRtbF9kb2N1bWVudDoKICAgIHRvYzogeWVzCiAgICBjb2RlX2Rvd25sb2FkOiB0cnVlCiAgICBjb2RlX2ZvbGRpbmc6ICJzaG93IgotLS0KCkFkYXB0ZWQgZnJvbSBhIFtibG9nCnBvc3RdKGh0dHBzOi8vc2N0eW5lci5naXRodWIuaW8vcmVkb2luZy1ncmFwaHMuaHRtbCkgYnkgW1NhbQpUeW5lcl0oaHR0cHM6Ly9zY3R5bmVyLmdpdGh1Yi5pbykKCmBgYHtyIGdsb2JhbF9vcHRpb25zLCBpbmNsdWRlID0gRkFMU0V9CmtuaXRyOjpvcHRzX2NodW5rJHNldChjb2xsYXBzZSA9IFRSVUUsIGZpZy5hbGlnbiA9ICJjZW50ZXIiKQpgYGAKCmBgYHtyLCBtZXNzYWdlID0gRkFMU0V9CmxpYnJhcnkoZHBseXIpCmxpYnJhcnkodGlkeXIpCmxpYnJhcnkocm5vYWEpCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShnZ3JlcGVsKQpsaWJyYXJ5KHBhdGNod29yaykKdGhlbWVfc2V0KHRoZW1lX21pbmltYWwoKSArIHRoZW1lKHRleHQgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDE2KSkpCmBgYAoKIyMgR2V0dGluZyB0aGUgRGF0YSBmcm9tIE5PQUEKCkNsaW1hdGUgZGF0YSBpcyBhdmFpbGFibGUgdmlhIGEgW3dlYiBzZXJ2aWNlCkFQSV0oaHR0cHM6Ly93d3cubmNkYy5ub2FhLmdvdi9jZG8td2ViL3dlYnNlcnZpY2VzL3YyKSBmcm9tCltOT0FBXShodHRwczovL3d3dy5ub2FhLmdvdi8pLgoKVGhlCltgcm5vYWFgXShodHRwczovL2NyYW4uci1wcm9qZWN0Lm9yZy9wYWNrYWdlPXJub2FhKQpwYWNrYWdlIHByb3ZpZGVzIGFuIGludGVyZmFjZSB0byB0aGlzIEFQSS4gIEEKW3R1dG9yaWFsXShodHRwczovL2RvY3Mucm9wZW5zY2kub3JnL3Jub2FhL2FydGljbGVzL3Jub2FhLmh0bWwpIGlzCmF2YWlsYWJsZSBmcm9tIFtyT3BlblNjaV0oaHR0cHM6Ly9yb3BlbnNjaS5vcmcpLgoKVGhpcyBjb2RlIHdvdWxkIGJlIHVzZWQgdG8gZmluZCB0aGUgbW9zdCBzdWl0YWJsZSBJb3dhIENpdHkgc3RhdGlvbjoKCmBgYHtyLCBldmFsID0gRkFMU0V9Cm9wdGlvbnMoIjxrZXk+IikgIyMgbWF5IG5vdCBiZSBuZWVkZWQKc3RhdCA8LSBnaGNuZF9zdGF0aW9ucygpCmljX3N0YXQgPC0gZmlsdGVyKHN0YXQsIGdyZXBsKCJJT1dBIENJVFkiLCBuYW1lKSAmIGVsZW1lbnQgPT0gIlNOT1ciKQpzZWxlY3QoaWNfc3RhdCwgbmFtZSwgaWQsIGVuZHNfd2l0aCgieWVhciIpKQpgYGAKClRoZSBiZXN0IHN0YXRpb24gc2VlbXMgdG8gYmUKCmBgYHtyfQppY19pZCA8LSAiVVNDMDAxMzQxMDEiCmBgYAoKVGhlIGN1cnJlbnQgZGF0YSBmb3IgdGhpcyBzdGF0aW9uIGNhbiB0aGVuIGJlIG9idGFpbmVkIGZyb20gTk9BQToKCmBgYHtyLCBldmFsID0gRkFMU0V9CmljX2RhdGFfcmF3IDwtIGdoY25kKGljX2lkLCByZWZyZXNoID0gVFJVRSkKIyMgZ29lcyBpbnRvIGEgY2FjaGUgZmlsZSBhcyBjc3YKIyMgL2hvbWUvbHVrZS8uY2FjaGUvcm5vYWEvZ2hjbmQvVVNDMDAxMzQxMDEuZGx5CndyaXRlLmNzdihpY19kYXRhX3Jhdywgcm93Lm5hbWVzID0gRkFMU0UsIGZpbGUgPSAiaWNfbm9hYS5jc3YiKQojIyBnemlwIGFuZCBtb3ZlIHRvIGRhdGEgZGlyZWN0b3J5CmBgYAoKSSBoYXZlIGRvbmUgdGhlc2Ugc3RlcHMgYW5kIHN0b3JlZCB0aGUgZGF0YSBsb2NhbGx5OyBJIHdpbGwgdXBkYXRlCnRoZSBkYXRhIGFzIHRoZSBjb3Vyc2UgcHJvZ3Jlc3Nlcy4KCkRvd25sb2FkIHRoZSBsb2NhbCBkYXRhIGNvcHk6CgpgYGB7cn0KaWYgKCEgZmlsZS5leGlzdHMoImljX25vYWEuY3N2Lmd6IikpCiAgICBkb3dubG9hZC5maWxlKCJodHRwOi8vd3d3LnN0YXQudWlvd2EuZWR1L35sdWtlL2RhdGEvaWNfbm9hYS5jc3YuZ3oiLAogICAgICAgICAgICAgICAgICAiaWNfbm9hYS5jc3YuZ3oiKQojIyBtYWtlIHN1cmUgZGF0YSBpcyB1cCB0byBkYXRlCmljX2RhdGFfcmF3IDwtIGFzX3RpYmJsZShyZWFkLmNzdigiaWNfbm9hYS5jc3YuZ3oiLCBoZWFkID0gVFJVRSkpCmBgYAoKCiMjIFByb2Nlc3NpbmcgdGhlIERhdGEKClRoZSByYXcgZGF0YToKCmBgYHtyfQppY19kYXRhX3JhdwpgYGAKClRoZSBkb2N1bWVudGF0aW9uIGZvciB0aGUgZGF0YSBpcyBhdmFpbGFibGUgYXQKCj4gaHR0cHM6Ly93d3cxLm5jZGMubm9hYS5nb3YvcHViL2RhdGEvZ2hjbi9kYWlseS9yZWFkbWUudHh0CgpUaGUgYXZhaWxhYmxlIF9lbGVtZW50c18sIG9yIG1lYXN1cmVtZW50cywgYXJlOgoKYGBge3J9CnVuaXF1ZShpY19kYXRhX3JhdyRlbGVtZW50KQpgYGAKClRoZSBwcmltYXJ5IG9uZXMgYXJlOgoKYGBgClBSQ1AgPSBQcmVjaXBpdGF0aW9uICh0ZW50aHMgb2YgbW0pClNOT1cgPSBTbm93ZmFsbCAobW0pClNOV0QgPSBTbm93IGRlcHRoIChtbSkKVE1BWCA9IE1heGltdW0gdGVtcGVyYXR1cmUgKHRlbnRocyBvZiBkZWdyZWVzIEMpClRNSU4gPSBNaW5pbXVtIHRlbXBlcmF0dXJlICh0ZW50aHMgb2YgZGVncmVlcyBDKQpgYGAKClRoZSB2YWx1ZXMgZm9yIGRheSAkayQgYXJlIHN0b3JlZCBpbiB0aGUgdmFyaWFibGUgYFZBTFVFYCRrJC4gQSBzZWxlY3Rpb246CgpgYGB7cn0KaWNfZmViMjAxOSA8LSBmaWx0ZXIoaWNfZGF0YV9yYXcsIHllYXIgPT0gMjAxOSwgbW9udGggPT0gMikKc2VsZWN0KGljX2ZlYjIwMTksIGVsZW1lbnQsIFZBTFVFMSwgVkFMVUUyOCwgVkFMVUUzMSkKYGBgCgpXZSBuZWVkIHRvOgoKLSBSZXNoYXBlIHRoZSBkYXRhLgotIERlYWwgd2l0aCB0aGUgYm9ndXMgZGF5cy4KClN0YXJ0IGJ5IHNlbGVjdGluZyB0aGUgY29sdW1ucyB3ZSBuZWVkOgoKYGBge3J9CmljX2RhdGEgPC0gc2VsZWN0KGljX2RhdGFfcmF3LCB5ZWFyLCBtb250aCwgZWxlbWVudCwgc3RhcnRzX3dpdGgoIlZBTFVFIikpCmljX2RhdGEKYGBgCgpSZXNoYXBlIGZyb20gKHZlcnkpIHdpZGUgdG8gKHRvbykgbG9uZzoKYGBge3J9CmljX2RhdGEgPC0gcGl2b3RfbG9uZ2VyKGljX2RhdGEsIFZBTFVFMSA6IFZBTFVFMzEsCiAgICAgICAgICAgICAgICAgICAgICAgIG5hbWVzX3RvID0gImRheSIsIHZhbHVlc190byA9ICJ2YWx1ZSIpCmljX2RhdGEKYGBgCgpFeHRyYWN0IHRoZSBkYXkgYXMgYSBudW1iZXI6CgpgYGB7cn0KaWNfZGF0YSA8LSBtdXRhdGUoaWNfZGF0YSwgZGF5ID0gcmVhZHI6OnBhcnNlX251bWJlcihkYXkpKQppY19kYXRhCmBgYAoKUmVzaGFwZSBmcm9tIHRvbyBsb25nIHRvIHRpZHkgd2l0aCBvbmUgcm93IHBlciBkYXksIGtlZXBpbmcgb25seSB0aGUKcHJpbWFyeSB2YXJpYWJsZXM6CgpgYGB7cn0KY29yZXZhcnMgPC0gYygiVE1BWCIsICJUTUlOIiwgIlBSQ1AiLCAiU05PVyIsICJTTldEIikKaWNfZGF0YSA8LSBmaWx0ZXIoaWNfZGF0YSwgZWxlbWVudCAlaW4lIGNvcmV2YXJzKQppY19kYXRhIDwtIHBpdm90X3dpZGVyKGljX2RhdGEsIG5hbWVzX2Zyb20gPSAiZWxlbWVudCIsIHZhbHVlc19mcm9tID0gInZhbHVlIikKaWNfZGF0YQpgYGAKCkFkZCBhIGBkYXRlYCB2YXJpYWJsZSBmb3IgcGxvdHRpbmcgYW5kIHRvIGhlbHAgZ2V0IHJpZCBvZiBib2d1cyBkYXlzOgoKYGBge3J9CmljX2RhdGEgPC0gbXV0YXRlKGljX2RhdGEsIGRhdGUgPSBsdWJyaWRhdGU6Om1ha2VfZGF0ZSh5ZWFyLCBtb250aCwgZGF5KSkKZmlsdGVyKGljX2RhdGEsIHllYXIgPT0gMjAxOSwgbW9udGggPT0gMiwgZGF5ID49IDI3KQpgYGAKCkdldCByaWQgb2YgdGhlIGJvZ3VzIGRheXM6CmBgYHtyfQppY19kYXRhIDwtIGZpbHRlcihpY19kYXRhLCAhIGlzLm5hKGRhdGUpKQpgYGAKCk1ha2UgdW5pdHMgbW9yZSBzdGFuZGFyZCAoQW1lcmljYW4pOgoKYGBge3J9Cm1tMmluIDwtIGZ1bmN0aW9uKHgpIHggLyAyNS40CkMyRiA8LSBmdW5jdGlvbih4KSAzMiArIDEuOCAqIHgKaWNfZGF0YSA8LSB0cmFuc211dGUoaWNfZGF0YSwgeWVhciwgbW9udGgsIGRheSwgZGF0ZSwKICAgICAgICAgICAgICAgICAgICAgUFJDUCA9IG1tMmluKFBSQ1AgLyAxMCksCiAgICAgICAgICAgICAgICAgICAgIFNOT1cgPSBtbTJpbihTTk9XKSwKICAgICAgICAgICAgICAgICAgICAgU05XRCA9IG1tMmluKFNOV0QpLAogICAgICAgICAgICAgICAgICAgICBUTUlOID0gQzJGKFRNSU4gLyAxMCksCiAgICAgICAgICAgICAgICAgICAgIFRNQVggPSBDMkYoVE1BWCAvIDEwKSkKYGBgCgoKIyMgVGhlIFdpbnRlcnMgb2YgMjAxOC8yMDE5IGFuZCAyMDIwLzIwMjEKCiMjIyBEYWlseSBTbm93ZmFsbAoKRXh0cmFjdCB0aGUgZGF0YSBmb3IgZGF5cyBiZXR3ZWVuIE5vdmVtYmVyIDEgYW5kIE1hcmNoIDMxOgoKYGBge3J9CncxOCA8LSBmaWx0ZXIoaWNfZGF0YSwgZGF0ZSA+PSAiMjAxOC0xMS0wMSIgJiBkYXRlIDw9ICIyMDE5LTAzLTMxIikKdzIwIDwtIGZpbHRlcihpY19kYXRhLCBkYXRlID49ICIyMDIwLTExLTAxIiAmIGRhdGUgPD0gIjIwMjEtMDMtMzEiKQpgYGAKCkEgZmlyc3Qgc2V0IG9mIHBsb3RzOgoKYGBge3IsIGZpZy53aWR0aCA9IDl9CnAxIDwtIGdncGxvdCh3MTgsIGFlcyh4ID0gZGF0ZSkpICsgZ2VvbV9jb2woYWVzKHkgPSBTTk9XKSkKcDIgPC0gZ2dwbG90KHcxOCwgYWVzKHggPSBkYXRlKSkgKyBnZW9tX2xpbmUoYWVzKHkgPSBTTldEKSkKcDEgKyBwMgpgYGAKCkZvciBkYWlseSBzbm93ZmFsbCBpdCBpcyBwcm9iYWJseSByZWFzb25hYmxlIHRvIHJlcGxhY2UgbWlzc2luZyB2YWx1ZXMKYnkgemVyb3MuCgpgYGB7cn0KdzE4IDwtIG11dGF0ZSh3MTgsIFNOT1cgPSBpZmVsc2UoaXMubmEoU05PVyksIDAsIFNOT1cpKQp3MjAgPC0gbXV0YXRlKHcyMCwgU05PVyA9IGlmZWxzZShpcy5uYShTTk9XKSwgMCwgU05PVykpCmBgYAoKRm9yIHNub3cgZGVwdGggaXQgbWF5IG1ha2Ugc2Vuc2UgdG8gX2ZpbGwgZm9yd2FyZF8sIHJlcGxhY2luZyBtaXN1c2luZwp2YWx1ZXMgYnkgdGhlIHByZXZpb3VzIG5vbi1taXNzaW5nIHZhbHVlLiBVc2luZyB0aGUgYGZpbGxgIGZ1bmN0aW9uCmZyb20gYHRpZHlyYCBpcyBvbmUgd2F5IHRvIGRvIHRoaXMuCgpgYGB7cn0KdzE4IDwtIGZpbGwodzE4LCBTTldEKQp3MjAgPC0gZmlsbCh3MjAsIFNOV0QpCmBgYAoKVGhlIHBsb3Qgb2Ygc25vdyBkZXB0aCBsb29rcyBhIGJpdCBiZXR0ZXIgd2l0aCB0aGlzIGNoYW5nZTsgdGhlIGF4aXMKcmFuZ2UgZm9yIHRoZSBwbG90IG9mIGRhaWx5IHNub3dmYWxsIGlzIGFsc28gYWZmZWN0ZWQuIEFsc28sIHRoZQp3YXJuaW5ncyBhcmUgZWxpbWluYXRlZC4KCmBgYHtyLCBmaWcud2lkdGggPSA5fQpwMSA8LSBnZ3Bsb3QodzE4LCBhZXMoeCA9IGRhdGUpKSArIGdlb21fY29sKGFlcyh5ID0gU05PVykpCnAyIDwtIGdncGxvdCh3MTgsIGFlcyh4ID0gZGF0ZSkpICsgZ2VvbV9saW5lKGFlcyh5ID0gU05XRCkpCihwMSArIGxhYnModGl0bGUgPSAiMjAxOC8xOSIpKSArIHAyCihwMSAlKyUgdzIwICsgbGFicyh0aXRsZSA9ICIyMDIwLzIxIikpICsgKHAyICUrJSB3MjApCmBgYAoKT2Z0ZW4gc25vd2ZhbGwgaXMgdmlzdWFsaXplZCBpbiB0ZXJtcyBvZiB0aGUgY3VtdWxhdGl2ZSBzbm93ZmFsbCBmb3IKdGhlIHNlYXNvbiB1cCB0byBlYWNoIGRhdGU6CgpgYGB7cn0KcCA8LSBnZ3Bsb3QodzE4LCBhZXMoeCA9IGRhdGUpKSArIGdlb21fbGluZShhZXMoeSA9IGN1bXN1bShTTk9XKSkpCnAgKyBsYWJzKHRpdGxlID0gIjIwMTgvMTkiKQpwICUrJSB3MjAgKyBsYWJzKHRpdGxlID0gIjIwMjAvMjEiKQpgYGAKCgojIyMgQ29tcGFyaW5nIHRvIExvbmcgVGVybSBBdmVyYWdlcwoKQ29tcHV0ZSB0aGUgYXZlcmFnZSBzbm93ZmFsbCBmb3IgZWFjaCBkYXkgb2YgdGhlIHllYXIgYW5kIG1lcmdlIHRoZQpyZXN1bHRzIGZvciB0aGUgZGF5cyBpbiBgdzE4YCBpbnRvIGB3MThgIHdpdGggYSBgbGVmdF9qb2luYDoKCmBgYHtyLCBtZXNzYWdlID0gRkFMU0V9CmljX2F2ZyA8LSBncm91cF9ieShpY19kYXRhLCBtb250aCwgZGF5KSAlPiUKICAgIHN1bW1hcml6ZShhdl9zbm93ID0gbWVhbihTTk9XLCBuYS5ybSA9IFRSVUUpLAogICAgICAgICAgICAgIGF2X3Nud2QgPSBtZWFuKFNOV0QsIG5hLnJtID0gVFJVRSkpICU+JQogICAgdW5ncm91cCgpCncxOCA8LSBsZWZ0X2pvaW4odzE4LCBpY19hdmcsIGMoIm1vbnRoIiwgImRheSIpKQp3MjAgPC0gbGVmdF9qb2luKHcyMCwgaWNfYXZnLCBjKCJtb250aCIsICJkYXkiKSkKYGBgCgpEYWlseSBzbm93ZmFsbCBhbmQgc25vdyBjb3ZlciBmb3IgdGhlIDIwMTgvMTkgd2ludGVyIGFsb25nIHdpdGgKYXZlcmFnZXMgb3ZlciBhbGwgeWVhcnMgaW4gdGhlIHJlY29yZDoKCmBgYHtyLCBmaWcud2lkdGggPSA5fQpwMSA8LSBnZ3Bsb3QodzE4LCBhZXMoeCA9IGRhdGUpKSArCiAgICBnZW9tX2NvbChhZXMoeSA9IFNOT1cpKSArCiAgICBnZW9tX2xpbmUoYWVzKHkgPSBhdl9zbm93KSwgY29sb3IgPSAiYmx1ZSIpCnAyIDwtIGdncGxvdCh3MTgsIGFlcyh4ID0gZGF0ZSkpICsKICAgIGdlb21fbGluZShhZXMoeSA9IFNOV0QpKSArCiAgICBnZW9tX2xpbmUoYWVzKHkgPSBhdl9zbndkKSwgY29sb3IgPSAiYmx1ZSIpCihwMSArIGxhYnModGl0bGUgPSAiMjAxOC8xOSIpKSArIHAyCihwMSAlKyUgdzIwICsgbGFicyh0aXRsZSA9ICIyMDIwLzIxIikpICsgKHAyICUrJSB3MjApCmBgYAoKQ3VtdWxhdGl2ZSBzbm93ZmFsbCBmb3IgdGhlIDIwMTgvMTkgd2ludGVyIGFsb25nIHdpdGggYXZlcmFnZQpjdW11bGF0aXZlIHNub3dmYWxsIG92ZXIgYWxsIHllYXJzIGluIHRoZSByZWNvcmQ6CgpgYGB7cn0KcCA8LSBnZ3Bsb3QodzE4LCBhZXMoeCA9IGRhdGUpKSArCiAgICBnZW9tX2xpbmUoYWVzKHkgPSBjdW1zdW0oU05PVykpKSArCiAgICBnZW9tX2xpbmUoYWVzKHkgPSBjdW1zdW0oYXZfc25vdykpLCBjb2xvciA9ICJibHVlIikKcCArIGxhYnModGl0bGUgPSAiMjAxOC8xOSIpCnAgJSslIHcyMCArIGxhYnModGl0bGUgPSAiMjAyMC8yMSIpCmBgYAoKCiMjIyBBIENvbWJpbmVkIEdyYXBoCgpgYGB7cn0KdzE4ZGF5IDwtIHRyYW5zbXV0ZSh3MTgsIGRhdGUsIHNub3cgPSBTTk9XLAogICAgICAgICAgICAgICAgICAgIGJhc2UgPSAiMjAxOC8xOSIsCiAgICAgICAgICAgICAgICAgICAgdHlwZSA9ICJEYWlseSBTbm93ZmFsbCIpCm5ybWRheSA8LSB0cmFuc211dGUodzE4LCBkYXRlLCBzbm93ID0gYXZfc25vdywKICAgICAgICAgICAgICAgICAgICBiYXNlID0gIk5vcm1hbCIsCiAgICAgICAgICAgICAgICAgICAgdHlwZSA9ICJEYWlseSBTbm93ZmFsbCIpCncxOGN1bSA8LSB0cmFuc211dGUodzE4ZGF5LCBkYXRlLCBzbm93ID0gY3Vtc3VtKHNub3cpLAogICAgICAgICAgICAgICAgICAgIGJhc2UgPSAiMjAxOC8xOSIsCiAgICAgICAgICAgICAgICAgICAgdHlwZSA9ICJDdW11bGF0aXZlIFNub3dmYWxsIikKbnJtY3VtIDwtIHRyYW5zbXV0ZShucm1kYXksIGRhdGUsIHNub3cgPSBjdW1zdW0oc25vdyksCiAgICAgICAgICAgICAgICAgICAgYmFzZSA9ICJOb3JtYWwiLAogICAgICAgICAgICAgICAgICAgIHR5cGUgPSAiQ3VtdWxhdGl2ZSBTbm93ZmFsbCIpCmBgYAoKYGBge3IsIGZpZy5oZWlnaHQgPSA4fQpnZ3Bsb3QodzE4LCBhZXMoeCA9IGRhdGUsIHkgPSBzbm93LCBjb2xvciA9IGJhc2UpKSArCiAgICBnZW9tX2xpbmUoZGF0YSA9IHcxOGN1bSkgKwogICAgZ2VvbV9saW5lKGRhdGEgPSBucm1jdW0pICsKICAgIGdlb21fc2VnbWVudChhZXMoeCA9IGRhdGUsIHhlbmQgPSBkYXRlLCB5ID0gMCwgeWVuZCA9IHNub3cpLAogICAgICAgICAgICAgICAgIHNpemUgPSAxLCBkYXRhID0gdzE4ZGF5KSArCiAgICBnZW9tX2xpbmUoZGF0YSA9IG5ybWRheSkgKwogICAgZmFjZXRfd3JhcCh+IHR5cGUsIG5jb2wgPSAxLCBzY2FsZXMgPSAiZnJlZV95IikgKwogICAgeGxhYihlbGVtZW50X2JsYW5rKCkpICsgeWxhYigiU25vd2ZhbGwgKGluY2hlcykiKSArCiAgICBnZ3RpdGxlKCIyMDE4LzE5IFNub3dmYWxsIGluIElvd2EgQ2l0eSIpICsKICAgIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSBjKCIyMDE4LzE5IiA9ICJuYXZ5IiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICJOb3JtYWwiID0gImZvcmVzdGdyZWVuIikpICsKICAgIGd1aWRlcyhjb2xvciA9IGd1aWRlX2xlZ2VuZChvdmVycmlkZS5hZXMgPSBsaXN0KHNpemUgPSAxLjUpKSkgKwogICAgdGhlbWVfYncoKSArCiAgICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAiYm90dG9tIiwgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGxvdC50aXRsZSA9IGVsZW1lbnRfdGV4dChmYWNlID0gImJvbGQiLCBoanVzdCA9IC41KSwKICAgICAgICAgIGFzcGVjdC5yYXRpbyA9IDAuNCkKYGBgCgojIyMgVG90YWwgU25vdyBGYWxsIGJ5IFllYXIgYW5kIE1vbnRoCgpJdCBpcyBoZWxwZnVsIHRvIGRlZmluZSBhIHZhcmlhYmxlIHRoYXQgc2hpZnRzIHRoZSB5ZWFyIHRvIGluY2x1ZGUgYQpmdWxsIHdpbnRlciwgbm90IHRoZSBlbmQgb2Ygb25lIGFuZCB0aGUgc3RhcnQgb2YgdGhlIG5leHQuCgpXaGF0IG1vbnRoIHNob3VsZCBiZSBzdGFydCB3aXRoPyBIZXJlIGFyZSB0aGUgY291bnRzIGZvciBtb250aHMgd2l0aApyZWNvcmRlZCBzbm93ZmFsbDoKCmBgYHtyfQpjb3VudChmaWx0ZXIoaWNfZGF0YSwgU05PVyA+IDApLCBtb250aCkKYGBgCgpUaGUgSnVseSBzbm93ZmFsbCBzZWVtcyBvZGQ6CgpgYGB7cn0KZmlsdGVyKGljX2RhdGEsIFNOT1cgPiAwLCBtb250aCA9PSA3KQpgYGAKClRoaXMgaXMgcHJvYmFibHkgYSByZWNvcmRpbmcgZXJyb3IuCgpGb3IgT2N0b2JlciwgdGhlIHllYXJzIHdpdGggIHJlY29yZGVkIHJlY29yZGVkIHNub3dmYWxsIGFyZSBzaG93biBieQoKYGBge3J9CmNvdW50KGZpbHRlcihpY19kYXRhLCBTTk9XID4gMCwgbW9udGggPT0gMTApLCB5ZWFyKQpgYGAKCkl0IHNlZW1zIHJlYXNvbmFibGUgdG8gc3RhcnQgdGhlIF93aW50ZXIgeWVhcl8gaW4gT2N0b2JlcjoKCmBgYHtyfQppY19kYXRhIDwtIG11dGF0ZShpY19kYXRhLCB3eWVhciA9IGlmZWxzZShtb250aCA+PSAxMCwgeWVhciwgeWVhciAtIDEpKQpgYGAKClRoZSB5ZWFybHkgdG90YWwgc25vd2ZhbGwgdmFsdWVzOgoKYGBge3J9CnRvdCA8LSBncm91cF9ieShpY19kYXRhLCB3eWVhcikgJT4lCiAgICBzdW1tYXJpemUodHNub3cgPSBzdW0oU05PVywgbmEucm0gPSBUUlVFKSkKdG90MTggPC0gZmlsdGVyKHRvdCwgd3llYXIgPT0gMjAxOCkKZ2dwbG90KHRvdCwgYWVzKHggPSB3eWVhciwgeSA9IHRzbm93KSkgKwogICAgZ2VvbV9saW5lKCkgKwogICAgZ2VvbV9wb2ludChkYXRhID0gdG90MTgsIGNvbG9yID0gInJlZCIpICsKICAgIGdlb21fdGV4dF9yZXBlbChkYXRhID0gdG90MTgsIGxhYmVsID0gIjIwMTgvMTkiLCBjb2xvciA9ICJyZWQiKQpgYGAKClNvIDIwMTgvMTkgaXMgbm90IGEgcmVjb3JkLCBidXQgc3RpbGwgcHJldHR5IGhpZ2g6CgpgYGB7cn0KdG90MTAwIDwtIG11dGF0ZShmaWx0ZXIodG90LCB3eWVhciA+PSAxOTE4KSwgcmFuayA9IG1pbl9yYW5rKGRlc2ModHNub3cpKSkKZmlsdGVyKHRvdDEwMCwgd3llYXIgPT0gMjAxOCkKYGBgCgpXaGF0IGFyZSB0aGUgbW9udGhzIHdpdGggdGhlIGhpZ2hlc3QgYXZlcmFnZSB0b3RhbCBzbm93ZmFsbD8KCmBgYHtyfQptdG90IDwtIGdyb3VwX2J5KGljX2RhdGEsIHllYXIsIHd5ZWFyLCBtb250aCkgJT4lCiAgICBzdW1tYXJpemUobXNub3cgPSBzdW0oU05PVywgbmEucm0gPSBUUlVFKSkgJT4lCiAgICB1bmdyb3VwKCkKZ3JvdXBfYnkobXRvdCwgbW9udGgpICU+JSBzdW1tYXJpemUoYXZfbXNub3cgPSBtZWFuKG1zbm93KSkKYGBgCgpGb3IgMjAxOC8xOToKCmBgYHtyfQpmaWx0ZXIobXRvdCwgd3llYXIgPT0gMjAxOCkKYGBgCkZvciAyMDIwLzIxOgoKYGBge3J9CmZpbHRlcihtdG90LCB3eWVhciA9PSAyMDIwKQpgYGAKCgojIyBBIERhdGEgUHJvY2Vzc2luZyBGdW5jdGlvbnMKCklmIHlvdSB3YW50IHRvIGxvb2sgYXQgZGF0YSBmb3IgYW5vdGhlciBjaXR5IHlvdSB3b3VsZCBuZWVkIHRvIGdvCnRocm91Z2ggdGhlIHNhbWUgc3RlcHMgc3RhcnRpbmcgd2l0aCBhIG5ldyByYXcgZGF0YSBmcmFtZS4KCkRlZmluaW5nIGEgZnVuY3Rpb24gdGhhdCBlbmNhcHN1bGF0ZXMgdGhlIHN0ZXBzIG1ha2VzIHRoaXMgZWFzeToKCjwhLS0gIyMgbm9saW50IHN0YXJ0OiBvYmplY3RfdXNhZ2UgLS0+CmBgYHtyfQpwcm9jZXNzRGF0YSA8LSBmdW5jdGlvbihkYXRhX3JhdykgewogICAgIyMgc2VsZWN0IGp1c3QgdGhlIGNvbHVtbnMgd2UgbmVlZDoKICAgIGRhdGEgPC0gc2VsZWN0KGRhdGFfcmF3LCB5ZWFyLCBtb250aCwgZWxlbWVudCwgc3RhcnRzX3dpdGgoIlZBTFVFIikpCgogICAgIyMgcmVzaGFwZSBmcm9tICh2ZXJ5KSB3aWRlIHRvICh0b28pIGxvbmcKICAgIGRhdGEgPC0gcGl2b3RfbG9uZ2VyKGRhdGEsIFZBTFVFMSA6IFZBTFVFMzEsCiAgICAgICAgICAgICAgICAgICAgICAgICBuYW1lc190byA9ICJkYXkiLCB2YWx1ZXNfdG8gPSAidmFsdWUiKQogICAgIyMgZXh0cmFjdCB0aGUgZGF5IGFzIGEgbnVtYmVyOgogICAgZGF0YSA8LSBtdXRhdGUoZGF0YSwgZGF5ID0gcmVhZHI6OnBhcnNlX251bWJlcihkYXkpKQoKICAgICMjIHJlc2hhcGUgZnJvbSB0b28gbG9uZyB0byB0aWR5IHdpdGggb25lIHJvdyBwZXIgZGF5LCBrZWVwaW5nCiAgICAjIyBvbmx5IHRoZSBwcmltYXJ5IHZhcmlhYmxlczoKICAgIGNvcmV2YXJzIDwtIGMoIlRNQVgiLCAiVE1JTiIsICJQUkNQIiwgIlNOT1ciLCAiU05XRCIpCiAgICBkYXRhIDwtIGZpbHRlcihkYXRhLCBlbGVtZW50ICVpbiUgY29yZXZhcnMpCiAgICBkYXRhIDwtIHBpdm90X3dpZGVyKGRhdGEsIG5hbWVzX2Zyb20gPSAiZWxlbWVudCIsIHZhbHVlc19mcm9tID0gInZhbHVlIikKCiAgICAjIyBhZGQgYSAnZGF0ZScgdmFyaWFibGUgYW5kIGdldCByaWQgb2YgYm9ndXMgZGF5cwogICAgZGF0YSA8LSBtdXRhdGUoZGF0YSwgZGF0ZSA9IGx1YnJpZGF0ZTo6bWFrZV9kYXRlKHllYXIsIG1vbnRoLCBkYXkpKQogICAgZGF0YSA8LSBmaWx0ZXIoZGF0YSwgISBpcy5uYShkYXRlKSkKCiAgICAjIyBtYWtlIHVuaXRzIG1vcmUgc3RhbmRhcmQgKEFtZXJpY2FuKQogICAgbW0yaW4gPC0gZnVuY3Rpb24oeCkgeCAvIDI1LjQKICAgIEMyRiA8LSBmdW5jdGlvbih4KSAzMiArIDEuOCAqIHgKICAgIGRhdGEgPC0gdHJhbnNtdXRlKGRhdGEsIHllYXIsIG1vbnRoLCBkYXksIGRhdGUsCiAgICAgICAgICAgICAgICAgICAgICBQUkNQID0gbW0yaW4oUFJDUCAvIDEwKSwKICAgICAgICAgICAgICAgICAgICAgIFNOT1cgPSBtbTJpbihTTk9XKSwKICAgICAgICAgICAgICAgICAgICAgIFNOV0QgPSBtbTJpbihTTldEKSwKICAgICAgICAgICAgICAgICAgICAgIFRNSU4gPSBDMkYoVE1JTiAvIDEwKSwKICAgICAgICAgICAgICAgICAgICAgIFRNQVggPSBDMkYoVE1BWCAvIDEwKSkKCiAgICAjIyBhZGQgd2ludGVyIHllYXIgdmFyaWFibGUgJ3d5ZWFyJyBzdGFydGluZyBpbiBPY3RvYmVyCiAgICBkYXRhIDwtIG11dGF0ZShkYXRhLCB3eWVhciA9IGlmZWxzZShtb250aCA+PSAxMCwgeWVhciwgeWVhciAtIDEpKQoKICAgIGRhdGEKfQpgYGAKPCEtLSAjIyBub2xpbnQgZW5kIC0tPgoKQSBxdWljayBzYW5pdHkgY2hlY2s6CgpgYGB7cn0Kc3RvcGlmbm90KGlkZW50aWNhbChpY19kYXRhLCBwcm9jZXNzRGF0YShpY19kYXRhX3JhdykpKQpgYGAKCk9mdGVuIGl0IGlzIGEgZ29vZCBpZGVhIHRvIGluY2x1ZGUgYSBzYW5pdHkgY2hlY2sgbGlrZSB0aGlzIGluIGEgY29kZQpjaHVuayBtYXJrZWQgYXMgYGluY2x1ZGUgPSBGQUxTRWAuCgoKIyMgT3RoZXIgRXhwbG9yYXRpb25zCgpgYGB7cn0KcGl2b3RfbG9uZ2VyKHcxOCwgVE1JTiB8IFRNQVgsIG5hbWVzX3RvID0gIndoaWNoIiwgdmFsdWVzX3RvID0gInZhbHVlIikgJT4lCiAgICBnZ3Bsb3QoYWVzKHggPSBkYXRlKSkgKwogICAgZ2VvbV9saW5lKGFlcyh5ID0gdmFsdWUsIGNvbG9yID0gd2hpY2gpKQpgYGAKCgpgYGB7cn0Kd2QxOCA8LSBmaWx0ZXIoaWNfZGF0YSwgZGF0ZSA+PSAiMjAxOC0xMC0wMSIpCmF2cyA8LSBncm91cF9ieShpY19kYXRhLCBtb250aCwgZGF5KSAlPiUKICAgIHN1bW1hcml6ZShhdnNub3cgPSBtZWFuKFNOT1csIG5hLnJtID0gVFJVRSksCiAgICAgICAgICAgICAgYXZzbndkID0gbWVhbihTTldELCBuYS5ybSA9IFRSVUUpLAogICAgICAgICAgICAgIGF2dG1pbiA9IG1lYW4oVE1JTiwgbmEucm0gPSBUUlVFKSwKICAgICAgICAgICAgICBhdnRtYXggPSBtZWFuKFRNQVgsIG5hLnJtID0gVFJVRSkpICU+JQogICAgdW5ncm91cCgpCndkMTggPC0gbGVmdF9qb2luKHdkMTgsIGF2cykKCmdncGxvdChmaWx0ZXIodzE4LCAhIGlzLm5hKFNOV0QpKSkgKwogICAgZ2VvbV9saW5lKGFlcyh4ID0gZGF0ZSwgeSA9IFNOV0QpKSArCiAgICBnZW9tX2xpbmUoYWVzKHggPSBkYXRlLCB5ID0gYXZfc253ZCksIGNvbCA9ICJibHVlIikKYGBgCgoKYGBge3J9CnBsb3QoU05PVyB+IGRhdGUsIGRhdGEgPSB3ZDE4LCB0eXBlID0gImgiKQpsaW5lcyh3ZDE4JGRhdGUsIHdkMTgkYXZzbm93LCBjb2wgPSAicmVkIikKbGluZXMoc21vb3RoLnNwbGluZSh3ZDE4JGRhdGUsIHdkMTgkYXZzbm93KSwgY29sID0gImJsdWUiKQoKcGxvdChTTldEIH4gZGF0ZSwgZGF0YSA9IGZpbGwod2QxOCwgU05XRCksIHR5cGUgPSAibCIpCmxpbmVzKHdkMTgkZGF0ZSwgd2QxOCRhdnNud2QsIGNvbCA9ICJyZWQiKQoKcGxvdChjdW1zdW0oaWZlbHNlKGlzLm5hKFNOT1cpLCAwLCBTTk9XKSkgfiBkYXRlLCBkYXRhID0gd2QxOCwgdHlwZSA9ICJsIikKbGluZXMod2QxOCRkYXRlLCBjdW1zdW0od2QxOCRhdnNub3cpLCBjb2wgPSAicmVkIikKYGBgCgpgYGB7cn0KdGhtIDwtIHRoZW1lX2J3KCkgKwogICAgdGhlbWUocGFuZWwuZ3JpZC5tYWpvci54ID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgICAgcGFuZWwuZ3JpZC5taW5vciA9IGVsZW1lbnRfYmxhbmsoKSkKCmdncGxvdCh3ZDE4KSArCiAgICBnZW9tX2NvbChhZXMoeCA9IGRhdGUsIHkgPSBTTk9XKSkgKwogICAgZ2VvbV9saW5lKGFlcyh4ID0gZGF0ZSwgeSA9IGF2c25vdyksIGNvbCA9ICJyZWQiKSArCiAgICBnZW9tX3Ntb290aChhZXMoeCA9IGRhdGUsIHkgPSBhdnNub3cpLCBzZSA9IEZBTFNFKSArCiAgICB0aG0KCmdncGxvdChmaWxsKHdkMTgsIFNOV0QpKSArCiAgICBnZW9tX2xpbmUoYWVzKHggPSBkYXRlLCB5ID0gU05XRCkpICsKICAgIGdlb21fbGluZShhZXMoeCA9IGRhdGUsIHkgPSBhdnNud2QpLCBjb2wgPSAicmVkIikgKwogICAgZ2VvbV9zbW9vdGgoYWVzKHggPSBkYXRlLCB5ID0gYXZzbndkKSwgc2UgPSBGQUxTRSkgKwogICAgdGhtCgpnZ3Bsb3Qod2QxOCkgKwogICAgZ2VvbV9saW5lKGFlcyh4ID0gZGF0ZSwgeSA9IGN1bXN1bShpZmVsc2UoaXMubmEoU05PVyksIDAsIFNOT1cpKSkpICsKICAgIGdlb21fbGluZShhZXMoeCA9IGRhdGUsIHkgPSBjdW1zdW0od2QxOCRhdnNub3cpKSwgY29sID0gInJlZCIpICsKICAgIHRobQpgYGAKCmBgYHtyfQpnZ3Bsb3QoZmlsdGVyKGljX2RhdGEsIHllYXIgPT0gMjAxOSksIGFlcyh4ID0gZGF0ZSkpICsKICAgIGdlb21fbGluZShhZXMoeSA9IFRNQVgsIGNvbG9yID0gIlRNQVgiKSkgKwogICAgZ2VvbV9saW5lKGFlcyh5ID0gVE1JTiwgY29sb3IgPSAiVE1JTiIpKSArCiAgICB0aG0KYGBgCgpgYGB7cn0KZDEzIDwtIGZpbHRlcihpY19kYXRhLCB3eWVhciA+PSAyMDEzLCBtb250aCA8PSA0IHwgbW9udGggPj0gMTApCmQxMyA8LSBtdXRhdGUoZDEzLCB3ZGF5ID0gYXMubnVtZXJpYyhkYXRlIC0gbHVicmlkYXRlOjptYWtlX2RhdGUod3llYXIsIDEwLCAxKSksCiAgICAgICAgICAgICAgU05PVyA9IGlmZWxzZShpcy5uYShTTk9XKSwgMCwgU05PVykpCmdncGxvdChkMTMpICsKICAgIGdlb21fY29sKGFlcyh4ID0gd2RheSwgeSA9IFNOT1cpLCB3aWR0aCA9IDMsIHBvc2l0aW9uID0gImlkZW50aXR5IikgKwogICAgZmFjZXRfd3JhcCh+IHd5ZWFyKSArIHRobQpnZ3Bsb3QoZDEzKSArIGdlb21fbGluZShhZXMoeCA9IHdkYXksIHkgPSBTTk9XKSkgKyBmYWNldF93cmFwKH4gd3llYXIpICsgdGhtCmBgYAoKYGBge3J9CmQxMyA8LSBncm91cF9ieShkMTMsIHd5ZWFyKSAlPiUKICAgIG11dGF0ZShzbm93X3RvdGFsID0gY3Vtc3VtKFNOT1cpKSAlPiUKICAgIHVuZ3JvdXAoKQpnZ3Bsb3QoZDEzKSArCiAgICBnZW9tX2xpbmUoYWVzKHggPSB3ZGF5LCB5ID0gc25vd190b3RhbCwgY29sb3IgPSBmYWN0b3Iod3llYXIpKSkgKyB0aG0KZ2dwbG90KGQxMykgKwogICAgZ2VvbV9saW5lKGFlcyh4ID0gd2RheSwgeSA9IHNub3dfdG90YWwpKSArCiAgICBmYWNldF93cmFwKH4gd3llYXIpICsgdGhtCmBgYAoKYGBge3J9CmQxMyA8LSBmaWxsKGQxMywgU05XRCkKZ2dwbG90KGQxMykgKyBnZW9tX2xpbmUoYWVzKHggPSB3ZGF5LCB5ID0gU05XRCkpICsgZmFjZXRfd3JhcCh+IHd5ZWFyKSArIHRobQpgYGAKCmBgYHtyfQp3IDwtIGZpbHRlcihpY19kYXRhLCB3eWVhciA+PSAyMDEzKQp3IDwtIG11dGF0ZSh3LCB3ZGF5ID0gZGF0ZSAtIGx1YnJpZGF0ZTo6bWFrZV9kYXRlKHd5ZWFyLCAxMCwgMSkpCncgPC0gZmlsdGVyKHcsIHdkYXkgPD0gMTgwKQp3IDwtIG11dGF0ZSh3LCB3ZGF5ID0gbHVicmlkYXRlOjptYWtlX2RhdGUoMjAxOCwgMTAsIDEpICsgd2RheSkKZ2dwbG90KHcpICsKICAgIGdlb21fbGluZShhZXMoeCA9IHdkYXksIHkgPSBUTUFYKSwgY29sb3IgPSAiYmx1ZSIpICsKICAgIGdlb21fbGluZShhZXMoeCA9IHdkYXksIHkgPSBUTUlOKSwgY29sb3IgPSAicmVkIikgKwogICAgZmFjZXRfd3JhcCh+IHd5ZWFyKQoKYGBgCgotIEJldHRlciB4IGF4aXMgbGFiZWwgZm9yIHdpbnRlciBkYXlzCi0gc2V0IHRoZW1lIGdsb2JhbGx5Ci0gcmUtY3JlYXRlIG9yaWdpbmFsIHBsb3QgaW4gc29tZSB2YXJpYW50cwotIGJldHRlciBmYWNldCBsYWJlbHMK