america_is_on_the_move.R

library(rgdal)    # read and store spatial data
library(dplyr)    # manipulate data frames for data wrangling
library(leaflet)  # connect to LeafletJS library from R
library(magrittr) # pipe workflows
library(openxlsx) # read xlsx files
library(readxl)   # read xls  files

# prepare the data ------------------------------------------------------------
# read MSA layer 
capture.output({
  base <- readOGR("data/cb_2016_us_csa_5m/cb_2016_us_csa_5m.shp",
                  layer = "cb_2016_us_csa_5m", GDAL1_integer64_policy = TRUE)
}) -> captured_output

# # read migration calculations within MSAs
# load('data/all_county_county_summarized_on_msa_level.rda')
# 
# # calculate the net (from source to target) for each MSA as a source
# msa_net_source_to_target <-
#   all_county_county_summarized_on_msa_level %>%
#   group_by(source) %>%
#   summarize(net_source_to_target = sum(net_source_to_target, na.rm = TRUE))

msa_csa <- 
  read_excel('data/list1.xls', skip = 2) %>%
  select(`CBSA Title`, `CSA Title`) %>%
  unique() %>%
  filter(!is.na(`CSA Title`)) %>%
  rename(msa = `CBSA Title`,
         csa = `CSA Title`)

msa_net_source_to_target <- 
  read.xlsx('data/metro-to-metro-ins-outs-nets-gross.xlsx', startRow = 2) %>%
  magrittr::extract(-1, ) %>%
  select(Metro.Statistical.Area.of.Geography.A,
         Metro.Statistical.Area.of.Geography.B,
         Net.Migration.from.Geography.B.to.Geography.A2) %>%
  rename(target               = Metro.Statistical.Area.of.Geography.A,
         source               = Metro.Statistical.Area.of.Geography.B,
         net_source_to_target = Net.Migration.from.Geography.B.to.Geography.A2) %>%
  mutate(target = gsub(' Metro Area', '', target),
         source = gsub(' Metro Area', '', source)) %>%
  filter(!is.na(net_source_to_target)) %>%
  select(source, target, net_source_to_target) %>%
  left_join(msa_csa, by = c('source' = 'msa')) %>% 
  select(-source) %>%
  rename(source = csa) %>%
  group_by(source) %>%
  summarize(net_source_to_target = sum(as.numeric(as.character(net_source_to_target)), na.rm = TRUE)) %>%
  filter(!is.na(source))



# add net calcualtions to the MSA layer
# remember about non ascii strings in names that need to be removed
base$msa_net_source_to_target <- 
  data.frame(source = iconv(base$NAME, from = 'UTF-8', to = 'ASCII', sub = "")) %>% 
  left_join(
    msa_net_source_to_target %>%
      mutate(source = iconv(source, from = 'UTF-8', to = 'ASCII', sub = "") %>% gsub("?", "", ., fixed = TRUE))) %>%  
  magrittr::use_series(net_source_to_target)



base <- base[
  !is.na(base$msa_net_source_to_target),
]
#dim(base)


base$msa_net_source_to_target <- -base$msa_net_source_to_target

# prepare palette for colouring the map
bins <- quantile(base$msa_net_source_to_target, probs = seq(0, 1, 0.25))
pal <- colorBin("YlOrRd",
                domain = base$msa_net_source_to_target,
                bins = bins)


# prepare the leaflet map -------------------------------------------------
# create a basemmap 

base_decreased <-
  base[ # you can't use filters, arranges and top_n on spatial objects....
    base$NAME %>% # so I refer to identificators of rows that corresponds to the top 15 net source to target
      is_in(msa_net_source_to_target %>%
              top_n(15, net_source_to_target) %>%
              use_series(source)),
    ]
base_increased <-
  base[ # you can't use filters, arranges and top_n on spatial objects....
    base$NAME %>% # so I refer to identificators of rows that corresponds to the top 15 net source to target
      is_in(msa_net_source_to_target %>% filter(!is.na(source)) %>%
              top_n(15, -net_source_to_target) %>%
              use_series(source)),
    ]

labels <- sprintf(
  "<strong>%s</strong><br/> Net migration change: %s",
  base$NAME, format(base$msa_net_source_to_target, big.mark = ",",scientific = FALSE)
) %>% lapply(htmltools::HTML)

labels_decreased <- sprintf(
  "<strong>%s</strong><br/> Net migration change: %s",
  base_decreased$NAME, format(base_decreased$msa_net_source_to_target, big.mark = ",",scientific = FALSE)
) %>% lapply(htmltools::HTML)


labels_increased <- sprintf(
  "<strong>%s</strong><br/> Net migration change: %s",
  base_increased$NAME, format(base_increased$msa_net_source_to_target, big.mark = ",",scientific = FALSE)
) %>% lapply(htmltools::HTML)


qpal <- colorQuantile("RdYlBu", base$msa_net_source_to_target, n = 5)

map <-
leaflet(base) %>%
  setView(lng = -98, lat = 40 , zoom = 4) %>%
  addTiles() %>%
  addPolygons(
              group       = 'All Metro Areas',
              color       = ~qpal(msa_net_source_to_target), 
              weight      = 1,
              smoothFactor= 0.5,
              opacity     = 0.2,
              fillOpacity = 0.5,
              # fillColor   = ~pal(msa_net_source_to_target),
              highlightOptions = highlightOptions(
                color        = "white",
                weight       = 2, 
                bringToFront = TRUE),
              label          = labels,
              labelOptions = labelOptions(
                style        = list("font-weight" = "normal", padding = "3px 8px"),
                textsize     = "15px",
                direction    = "auto")) %>%
  addPolygons(data = base_decreased,
              group = 'Top 15 Origin Metro Areas',
              color       = "#444444",
              weight      = 1,
              smoothFactor= 0.5,
              opacity     = 0.2,
              fillOpacity = 0.5,
              fillColor   = ~qpal(msa_net_source_to_target),
              highlightOptions = highlightOptions(
                color        = "white",
                weight       = 2,
                bringToFront = TRUE),
              label          = labels_decreased,
              labelOptions = labelOptions(
                style        = list("font-weight" = "normal", padding = "3px 8px"),
                textsize     = "15px",
                direction    = "auto")) %>%
  addPolygons(data = base_increased,
              group = 'Top 15 Destination Metro Areas',
              color       = "#444444",
              weight      = 1,
              smoothFactor= 0.5,
              opacity     = 0.2,
              fillOpacity = 0.5,
              fillColor   = ~qpal(msa_net_source_to_target),
              highlightOptions = highlightOptions(
                color        = "white",
                weight       = 2,
                bringToFront = TRUE),
              label          = labels_increased,
              labelOptions = labelOptions(
                style        = list("font-weight" = "normal", padding = "3px 8px"),
                textsize     = "15px",
                direction    = "auto")) %>%
  # addLegend(position = "bottomleft",
  #           pal     = qpal,
  #           values  = ~msa_net_source_to_target,
  #           title   = "Population Change",
  #           # labFormat = labelFormat(prefix = "$"),
  #           opacity = 0.6 ) %>%
  
  addLegend(position = "bottomleft",
            #pal     = qpal,
            values   = ~msa_net_source_to_target,
            title    = "Net migration",
            colors   = c("#df7f81",  "#f5cda3", "#f8f6d4", "#cfe3ea", "#90b5d1"),
            labels   = c("Loss", "","","", "Gain"),
            opacity  = 1.0 ) %>%
  addLayersControl(baseGroups = 
                     c("All Metro Areas",
                       "Top 15 Origin Metro Areas",
                       "Top 15 Destination Metro Areas"),
                  #overlayGroups = ,
                  options = layersControlOptions(collapsed = FALSE)) 
map

library(htmlwidgets)
saveWidget(map, file = "map.html")