Add site gages to moving state boundaries

gistfile1.txt

library(USAboundaries)  # contains data (shape files)
library(dplyr)          # for wrangling
library(lubridate)      # for wrangling dates
library(rgdal)          # for wrangling shape files / map data
library(maptools)       # for wrangling shape files / map data
library(ggplot2)        # for plotting
library(gganimate)      # for gif animation
library(viridis)        # for lovely color palettes
library(dataRetrieval)

# site_sum_all <- data.frame()
# for(i in stateCd$STUSAB[1:51]){
#   sites <- readNWISdata(service = "site",
#                         seriesCatalogOutput=TRUE,
#                         parameterCd="00060",
#                         stateCd = i)
#   sites_sum <- filter(sites, parm_cd == "00060",
#                       data_type_cd == "dv") %>%
#     select(site_no, station_nm, dec_lat_va, dec_long_va, begin_date, end_date, count_nu) %>%
#     data.frame() %>%
#     mutate(begin_date = as.Date(begin_date),
#            end_date = as.Date(end_date))
#   
#   site_sum_all <- bind_rows(site_sum_all, sites_sum)
#   
# }
# 
# site_sum_all$begin_year <- as.numeric(format(site_sum_all$begin_date, "%Y"))
# site_sum_all$end_year <- as.numeric(format(site_sum_all$end_date, "%Y"))
# 
# site_data <- data.frame()
# for(i in 1:nrow(site_sum_all)){
#   years <- data.frame(years = seq(site_sum_all$begin_year[i],site_sum_all$end_year[i]))
#   site_data_sub <- bind_cols(site_sum_all[rep(i,nrow(years)),], years)
#   site_data <- bind_rows(site_data, site_data_sub)
# }
#   
# saveRDS(site_data, file="D:/LADData/RCode/Archive/pics/site_data.rds")

site_data <- readRDS("D:/LADData/RCode/Archive/pics/site_data.rds")

get_sites <- function(date, site_data=site_data){

  site_data_sub <- filter(site_data, years == as.numeric(format(as.Date(date),"%Y"))) %>%
    filter(!is.na(dec_lat_va)) %>%
    filter(!is.na(dec_long_va))

  coordinates(site_data_sub) <- ~ dec_long_va + dec_lat_va
  proj4string(site_data_sub) <- CRS("+proj=longlat +ellps=GRS80 +no_defs")
  
  site_data_sub <- spTransform(site_data_sub, CRS("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"))

  # extract, then rotate & shift hawaii
  hawaii <-  site_data_sub[site_data_sub$dec_long_va < -3000000,]
  if(length(hawaii) > 0){
    hawaii <- elide(hawaii, rotate=-35)
    hawaii <- elide(hawaii, shift=c(5400000, -1400000))
    proj4string(hawaii) <- proj4string(site_data_sub)    
  }
  
  alaska <-  site_data_sub[site_data_sub$dec_lat_va > 1000000,]
  if(length(alaska) > 0){

    # if(as.numeric(format(as.Date(date),"%Y")) <= 2000){
    #   sp <- us_boundaries(as.Date(date), type = "state")
    # } else {
    #   sp <- us_boundaries(as.Date("2000-01-01"), type = "state")
    # }
    # sp_t <- spTransform(sp, CRS("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"))
    # rownames(sp_t@data) <- as.character(as.numeric(rownames(sp_t@data)) - 1)
    # 
    # alaskaSF <-  sp_t[grep("Alaska", sp_t$full_name, fixed = TRUE),]
    # 
    # alaskaSF <- elide(alaskaSF, rotate=-50)
    
    alaska <- elide(alaska, rotate=-50)
    if(nrow(alaska) > 1){
      alaska <- elide(alaska, scale=max(apply(bbox(alaska), 1, diff)) / 2.3)
    }
    alaska <- elide(alaska, shift=c(-2100000, -2500000))
    proj4string(alaska) <- proj4string(site_data_sub)
  }
  
  site_data_sub_noAKHI <- site_data_sub[site_data_sub$dec_lat_va < 1000000 &
                                          site_data_sub$dec_long_va > -3000000,]
  if(length(site_data_sub_noAKHI) != length(site_data_sub)){
    if(length(alaska) > 0){
      site_data_sub_t <- rbind(site_data_sub_noAKHI, alaska)
    }
    if(length(hawaii) > 0){
      site_data_sub_t <- rbind(site_data_sub_noAKHI, hawaii)
    }

  } else {
    site_data_sub_t <- site_data_sub_noAKHI
  }
  
  sites <- as.data.frame(coordinates(site_data_sub_t)) %>%
    mutate(year = as.numeric(format(as.Date(date),"%Y")))  
  
}

## ------------------------------------------------------------------------ ##
get_boundaries <- function(date, type = "state", site_data=site_data){
  # Get shapefile for US states on specified date
  if(as.numeric(format(as.Date(date),"%Y")) <= 2000){
    sp <- us_boundaries(as.Date(date), type = type)
  } else {
    sp <- us_boundaries(as.Date("2000-01-01"), type = type)
  }
  
  # Extract metadata. Includes info about the territory type for
  # each region (i.e., state, unorganized territory, etc.)
  metadata <- sp %>%
    as.data.frame() %>%
    mutate(id = as.character(id_num -1))# %>%
  # select(id, name, terr_type)
  
  # ----------------------------------------------
  # Convert coordinate system to Albers equal area, 
  # and then reposition Alaska and Hawaii.
  # 
  # thanks to Bob Rudis (hrbrmstr):
  # https://github.com/hrbrmstr/rd3albers
  sp_t <- spTransform(sp, CRS("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"))
  rownames(sp_t@data) <- as.character(as.numeric(rownames(sp_t@data)) - 1)
  
  # extract, then rotate, shrink & move alaska (and reset projection)
  alaska <-  
    if(type == "state"){
      sp_t[grep("Alaska", sp_t$full_name, fixed = TRUE),]
    } else if (type == "county"){
      sp_t[grep("Alaska", sp_t$state_terr, fixed = TRUE),]
    }
  
  alaska <- elide(alaska, rotate=-50)
  alaska <- elide(alaska, scale=max(apply(bbox(alaska), 1, diff)) / 2.3)
  alaska <- elide(alaska, shift=c(-2100000, -2500000))
  proj4string(alaska) <- proj4string(sp_t)
  
  # extract, then rotate & shift hawaii
  hawaii <-  
    if(type == "state"){
      sp_t[grep("Hawaii", sp_t$full_name, fixed = TRUE),]
    } else if (type == "county"){
      sp_t[grep("Hawaii", sp_t$state_terr, fixed = TRUE),]
    }
  
  hawaii <- elide(hawaii, rotate=-35)
  hawaii <- elide(hawaii, shift=c(5400000, -1400000))
  proj4string(hawaii) <- proj4string(sp_t)
  
  # remove old states and put new ones back in
  sp_t <-
    if(type == "state"){
      sp_t[!grepl("Alaska|Hawaii", sp_t$full_name),]
    } else if(type == "county") {
      sp_t[!grepl("Alaska|Hawaii", sp_t$state_terr),]
    }
  
  sp_t <- rbind(sp_t, alaska, hawaii, makeUniqueIDs = TRUE)
  # ----------------------------------------------  
  
  # get borders for each region
  borders <- sp_t %>%
    fortify() %>%
    as.data.frame() %>%
    mutate(year = year(date)) %>%
    left_join(metadata, by = "id")
}

create_timeline <- function(boundary_df){
  # determine how many year's worth of data will be plotted
  years <- unique(boundary_df$year)
  
  # create sequence that spans the width of the map (from min to max longitude),
  # with evenly spaced intervals for years
  timeline_range <- with(boundary_df, 
                         seq(from = min(long),
                             to = max(long),
                             length.out = length(years)))
  
  # determine latitude coordinate for positioning timeline.
  # (max latitude * 1.2 places timeline just above map)
  timeline_y <- max(boundary_df$lat) * 1.2
  
  # collect timeline data into a df
  timeline <- data.frame(
    x = timeline_range,
    y = timeline_y,
    year = years
  )
  
  return(timeline)
}

# specify dates for which historic map data is to be gathered
state_dates <- seq(as.Date("1864-01-01"), as.Date("2016-01-01"), by = "year")

# for faster rendering, try plotting a reduced date range
# state_dates <- seq(as.Date("1784-01-01"), as.Date("1789-01-01"), by = "year")

# get state-level data for every date, bind all data into one df,
# and set factor levels for territory type
states <- lapply(state_dates, 
                 function(x) get_boundaries(x, type = "state")) %>%
  do.call(rbind, .) %>%
  mutate(
    terr_type = factor(terr_type,
                       levels = c("State", 
                                  "District of Columbia", 
                                  "Territory",
                                  "Unorganized Territory", 
                                  "Other")))

sites <- lapply(state_dates, 
                 function(x) get_sites(x,site_data)) %>%
  do.call(rbind, .)

# create timeline for state data
state_timeline <- create_timeline(states)


# plot the map
p_state <- ggplot() +
  # Add territory polygons
  geom_polygon(aes(x = long, y = lat, group = group, fill = terr_type, frame = year),
               data = states,
               alpha = 0.5) +
  geom_point(aes(x = dec_long_va, y=dec_lat_va, frame = year), 
             data = sites, colour = "red", size = 0.5) +
  # Add timeline
  geom_segment(aes(x = min(x), xend = max(x),
                   y = y, yend = y,
                   frame = year),
               data = state_timeline,
               colour = "grey50", size = 1) +
  geom_point(aes(x = x, y = y, frame = year),
             data = state_timeline,
             colour = "grey50", size = 3) +
  
  # styling  
  coord_equal() +
  scale_fill_viridis(discrete = TRUE) +
  labs(title = "Map of USGS Streamgages\nJanuary 1,",
       fill = "") +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        legend.position = "bottom")

gg_animate(p_state, interval = 0.1)