Create a gist now

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment