<p>Back in 2014 I blogged about <a target="_blank" rel="nofollow" href="https://rud.is/b/2014/11/26/visualizing-historical-most-likely-first-snowfall-dates-for-u-s-regions/">first snowfall dates</a> for a given U.S. state. It’s April 1, 2017 and we’re slated to get 12-18″ of snow up here in Main and <a target="_blank" rel="nofollow" href="http://twitter.com/mrshrbrmstr">@mrshrbrmstr</a> asked how often this — snow in May — has occurred near us.</p>
<p>As with all of these “R⁴ posts, expository is minimal and the focus is generally to demonstrate one small concept. </p>
<p>What I’ve done here (first) is make a full <code>tidyverse</code> update to the <a target="_blank" rel="nofollow" href="https://github.com/hrbrmstr/snowfirst"><code>snowfirst</code></a> code posted in the aforementioned blog post. You’ll need to clone that repo if you’re trying to work verbatim from the code below (otherwise just change file path code):</p>
<pre><code>library(rprojroot)
library(stringi)
library(hrbrthemes)
library(tidyverse)

pre &lt;- find_rstudio_root_file()

# Get and read in Maine precip ------------------------------------------------------

URL &lt;- &quot;http://cdiac.ornl.gov/ftp/ushcn_daily/state17_ME.txt.gz&quot;
fil &lt;- file.path(pre, &quot;data&quot;, basename(URL))
if (!file.exists(fil)) download.file(URL, fil)

read_fwf(file = fil,
         col_positions = fwf_widths(c(6, 4, 2, 4, rep(c(5, 1, 1, 1), 31)),
                                    col_names = c(&quot;coop_id&quot;, &quot;year&quot;, &quot;month&quot;, &quot;element&quot;,
                                                  flatten_chr(map(1:31, ~paste(&quot;r_&quot;, c(&quot;v&quot;, &quot;fm&quot;, &quot;fq&quot;, &quot;fs&quot;),
                                                                               .x, sep=&quot;&quot;))))),
         col_types = paste0(&quot;ciic&quot;, paste0(rep(&quot;iccc&quot;, 31), collapse=&quot;&quot;), collapse=&quot;&quot;),
         na = c(&quot;&quot;, &quot;NA&quot;, &quot;-&quot;, &quot;-9999&quot;)) %>%
  gather(day, value, starts_with(&quot;r_v&quot;)) %>%
  select(-starts_with(&quot;r_&quot;)) %>%
  mutate(day = as.numeric(stri_replace_first_fixed(day, &quot;r_v&quot;, &quot;&quot;))) %>%
  mutate(date = sprintf(&quot;%s-%02d-%02d&quot;, year, month, day)) -> daily_wx

# Read in stations ------------------------------------------------------------------

URL &lt;- &quot;http://cdiac.ornl.gov/ftp/ushcn_daily/ushcn-stations.txt&quot;
fil &lt;- file.path(pre, &quot;data&quot;, basename(URL))
if (!file.exists(fil)) download.file(URL, fil)

read_fwf(file = file.path(pre, &quot;data&quot;, &quot;ushcn-stations.txt&quot;),
         col_positions = fwf_widths(c(6, 9, 10, 7, 3, 31, 7, 7, 7, 3),
                                    col_names = c(&quot;coop_id&quot;, &quot;latitude&quot;, &quot;longitude&quot;,
                                                  &quot;elevation&quot;, &quot;state&quot;, &quot;name&quot;,
                                                  &quot;component_1&quot;, &quot;component_2&quot;,
                                                  &quot;component_3&quot;, &quot;utc_offset&quot;)),
         col_types = &quot;cdddcccccc&quot;) -> stations

closestStation &lt;- function(stations, lat, lon, restrict_to = NULL) {
  if (!is.null(restrict_to)) stations &lt;- filter(stations, state == restrict_to)
  index &lt;- which.min(sqrt((stations$latitude-lat)^2 +
                            (stations$longitude-lon)^2))
  stations[index,]
}

# compute total snow amounts per month ----------------------------------------------

(near_me &lt;- closestStation(stations, 43.2672, -70.8617, restrict_to=&quot;ME&quot;))</code></pre>
<p>Now that we have the data, the short lesson here is just exposing the fact that you can get blank facets for free with ggplot2. I’m pointing this out as many folks seem to not like reading R documentation or miss things in said documentation (in fact, I had to be instructed today by <a target="_blank" rel="nofollow" href="http://twitter.com/thomasp85">@thomasp85</a> about a <code>ggplot2</code> theme element setting that I didn’t know about and <em>should have</em> since I do try to keep up).</p>
<pre><code>filter(daily_wx, coop_id == near_me$coop_id, element==&quot;SNOW&quot;, value>0) %>%
  count(year, month, wt=value) %>%
  ungroup() %>%
  mutate(
    n = n / 10, # readings are in 10ths of inches
    date = as.Date(sprintf(&quot;%s-%02d-01&quot;, year, month)),
    month_name = lubridate::month(date, TRUE, FALSE)
  ) %>%
  ggplot(aes(x=date, y=n)) +
  geom_segment(aes(xend=date, yend=0), size=0.75, color=&quot;#9ecae1&quot;) +
  scale_y_continuous(limits=c(0, 65)) +
  facet_wrap(~month_name, ncol=3, drop=FALSE, scales=&quot;free&quot;) +
  labs(x=NULL, y=&quot;inches&quot;, title=&quot;Total snowfall in a given month by year&quot;,
       subtitle=&quot;Data for Station id 176905 — Portland (Maine) Jetport&quot;) +
  theme_ipsum_rc(grid=&quot;Y&quot;, axis_text_size=8)</code></pre>
<p><img alt="" class="wp-image-5566" src="https://i2.wp.com/rud.is/b/wp-content/uploads/2017/04/Cursor_and___projects_snowfirst_-_master_-_RStudio.png?resize=780%2C714&ssl=1" data-recalc-dims="1"></p>
<p>Without <code>ggplot2</code> helping us out we would have had to do some work to have those no-value facets to show up. I also like how there are no x-axis labels since there’s no data. <code>ggplot2::facet_wrap()</code> has many, very granular options for customizing the appearance of facets:</p>
<pre><code>facet_wrap(facets, nrow = NULL, ncol = NULL, scales = &quot;fixed&quot;,
           shrink = TRUE, labeller = &quot;label_value&quot;, as.table = TRUE,
           switch = NULL, drop = TRUE, dir = &quot;h&quot;, strip.position = &quot;top&quot;)</code></pre>
<p>If you haven’t played with them, you can use this example to try them out.</p>
<h3>Fin</h3>
<p>Even though that visualization gets the message across, I kinda like this view a bit better:</p>
<pre><code>filter(daily_wx, coop_id == near_me$coop_id, element==&quot;SNOW&quot;, value>0) %>%
  count(year, month, wt=value) %>%
  ungroup() %>%
  mutate(n = n / 10) %>%
  complete(year, month=1:12) %>%
  mutate(
    date = as.Date(sprintf(&quot;%s-%02d-01&quot;, year, month)),
    month_name = factor(lubridate::month(date, TRUE, FALSE), levels=rev(month.name))
  ) %>% 
  ggplot(aes(year, month_name)) +
  geom_tile(aes(fill=n), color=&quot;#b2b2b2&quot;, size=0.15) +
  scale_x_continuous(expand=c(0,0.15), position=&quot;top&quot;) +
  viridis::scale_fill_viridis(name = &quot;Total inches&quot;, na.value=&quot;white&quot;) +
  labs(x=NULL, y=NULL, title=&quot;Total snowfall in a given month by year&quot;,
       subtitle=&quot;Data for Station id 176905 — Portland (Maine) Jetport&quot;) +
  theme_ipsum_rc(grid=&quot;&quot;, axis_text_size = 10) +
  guides(fill=guide_colourbar(label.position = &quot;top&quot;, direction = &quot;horizontal&quot;, title.vjust = 0)) +
  theme(legend.title = element_text(size=10)) +
  theme(legend.key.height = unit(0.5, &quot;lines&quot;)) +
  theme(legend.position = c(0.9, 1.25))</code></pre>
<p><img alt="" class="wp-image-5568" src="https://i1.wp.com/rud.is/b/wp-content/uploads/2017/04/Cursor_and___projects_snowfirst_-_master_-_RStudio-1.png?resize=780%2C213&ssl=1" data-recalc-dims="1"></p>
<p>The precision is lacking in the heatmap view, but you get a quick impression of when it has/hasn’t snowed. Plus you get to use <code>viridis</code> <img alt="💙" class="wp-smiley" src="https://s.w.org/images/core/emoji/2/72x72/1f499.png"></p>
<p>All the updated code in in the <a target="_blank" rel="nofollow" href="https://github.com/hrbrmstr/snowfirst"><code>snowfirst</code> repo</a>.</p>
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/kT-4aj0zaYU">