3 min read

Snowfall maps code

This code builds maps of snowfall since 2018 and places which currently have snow on the ground. You can see the updated version of these maps here.

# If you haven't installed these packages yet, do so now
library(tidyverse)
library(sf)
library(tmap)

# Read station locations
stations <- read_table("https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/ghcnd-stations.txt",
                       col_names = c("STNID","latitude","longitude","elevation","state","name","gsnflag",
                                     "hcnflag","wmoid","method")) %>%
  filter(state != "")
stations <- as_data_frame(state.abb) %>%
  filter(value != "AK" & value != "HI") %>% # remove Alaska and Hawaii
  inner_join(stations, by = c("value"="state")) # just keep the contiguous US stations

# You'll need to upload a shapefile of US states. This file is from the Census
state <- st_read("~/SHPfiles/cb_2017_us_state_500k/cb_2017_us_state_500k.shp")
state <- state[as.numeric(as.character(state$STATEFP))<57 &
                 state$STATEFP!="15" & state$STATEFP!="02",]

##  snowfall and snowdepth are in mm
# Get current year data
orig <- read_csv(paste0("https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/by_year/",
                        str_sub(Sys.Date(),1,4),
                        ".csv.gz"),
                 col_names = c("STNID","date","type","amount","flag1","flag2","flag3","flag4")) %>%
  inner_join(stations[,c(2)])

# Remove flags
d <- orig[,c(1:4)]
d$date <- as.Date(paste(str_sub(d$date,1,4),
                        str_sub(d$date,5,6),
                        str_sub(d$date,7,8), sep = "-"))

# Get a list of everywhere that has reported in the last 7 days
has.reported <- d %>%
  filter(date >= Sys.Date()-7) %>%
  group_by(STNID) %>% summarise()

# Fix obvious error
# Sometimes there are data entries which are obviously wrong. I fix one here
d$amount[d$STNID=="US1MIET0005" & d$date==as.Date("2018-10-06")] <- 0

has.snowed <- d %>%
  filter(type == "SNOW") %>%
  mutate(amount = amount*0.03937008) %>%
  inner_join(has.reported) %>%
  filter(date > as.Date("2018-07-31"),
         amount > 0) %>%
  group_by(STNID) %>%
  mutate(cumulative = cumsum(amount)) %>%
  filter(date == last(date)) %>%
  ungroup() %>%
  inner_join(stations) %>%
  st_as_sf(coords = c("longitude","latitude"), crs = 4326)

snow.on.ground <- d %>%
  filter(type == "SNWD",
         date >= Sys.Date()-7) %>%
  mutate(amount = amount*0.03937008) %>%
  filter(amount > 0) %>%
  group_by(STNID) %>%
  filter(date == max(date)) %>%
  ungroup() %>%
  inner_join(stations) %>%
  st_as_sf(coords = c("longitude","latitude"), crs = 4326)

has.reported.shp <- has.reported %>%
  inner_join(stations) %>%
  st_as_sf(coords = c("longitude","latitude"), crs = 4326)

tm.has.snowed <- tm_shape(state) +
  tm_polygons(alpha = 0) +
  tm_shape(has.reported.shp) +
  tm_dots(col = "grey") +
  tm_shape(has.snowed) +
  tm_bubbles(col = "blue", size = "cumulative", border.col = "blue",
             scale = 0.75, legend.size.is.portrait = T,
             title.size = "Cumulative inches") +
  tm_layout(title = "Places where it has snowed since August 1, 2018",
            title.position = c(.6,.95), frame = F, fontfamily = "serif",
            legend.position = c("left","bottom"), title.size = 0.8) +
  tm_credits(paste("Updated on", Sys.Date()), position = c(.75,.25)) +
  tm_credits("Universe: NOAA weather stations\nwhich have posted updates\nin the last week",
             position = c(.77,.1)) +
  tm_credits("Abnormal locations may be NOAA data entry errors.",
             position = c(0.07,0.01))
print(tm.has.snowed)

tm.snow.on.ground <- tm_shape(state) +
  tm_polygons(alpha = 0) +
  tm_shape(has.reported.shp) +
  tm_dots(col = "grey") +
  tm_shape(snow.on.ground) +
  tm_bubbles(col = "blue", border.col = "blue",
             size = "amount", scale = 0.75, legend.size.is.portrait = T,
             title.size = "Inches", size.lim = c(0.1,30)) +
  tm_layout(title = "Places with snow on the ground",
            title.position = c(.7,.95), frame = F, fontfamily = "serif",
            legend.position = c("left","bottom"), title.size = 0.8) +
  tm_credits(paste("Updated on", Sys.Date()), position = c(.75,.25)) +
  tm_credits("Universe: NOAA weather stations\nwhich have posted updates\nin the last week",
             position = c(.77,.1)) +
  tm_credits("Note: extreme outliers are not sized proportionally. Abnormal locations may be NOAA data entry errors.",
             position = c(0.07,0.01))
print(tm.snow.on.ground)