Maps using sf

maps_using_sf.R

# Lec 27: Friday 2018/11/9
library(tidyverse)
library(USAboundaries)
library(leaflet)
library(maps)
library(fivethirtyeight)
library(sf)

# Forget Lec26! sf package for maps is better!
USA_sf <- maps::map("state", plot = FALSE, fill = TRUE) %>%
  sf::st_as_sf() %>%
  sf::st_transform(4326)
pres_2016_trail_sf <- pres_2016_trail %>%
  sf::st_as_sf(coords = c("lng", "lat"), crs = 4326)

ggplot() +
  geom_sf(data = USA_sf) +
  geom_sf(data = pres_2016_trail_sf, aes(col = candidate), size = 2) +
  facet_wrap(~candidate)

# Download .zip file containing shapefiles, unzip the contents of the .zip file
# to a new directory/folder called "biketrails_folder".
url <- "http://download.massgis.digital.mass.gov/shapefiles/state/biketrails_arc.zip"
download.file(url, destfile = "biketrails_arc.zip")
unzip("biketrails_arc.zip", exdir = "biketrails_folder")

# After doing the above, in the RStudo Files panel navigate to the 
# "biketrails_folder" directory/folder. There should be 9 files with prefix
# "biketrails_arc". These are your shapefiles

# Read in shapefiles into R using 
bike_trails <- "biketrails_folder" %>%
  sf::read_sf()

# Plot this data! If you get an error that says "Error in grid.Call..."
# Try plotting a few more times; it for some reason works for me
ggplot() +
  geom_sf(data = bike_trails, col = "red")

# The above map lacks context however. Where are these bikepaths?
# Let's obtain a map of Massachusetts counties from the USAboundaries package
MA_counties <- 
  USAboundaries::us_counties(resolution = "high", states = c("massachusetts")) 

# Plot new layer!
ggplot() +
  geom_sf(data = bike_trails, col = "red") + 
  geom_sf(data = MA_counties) 

# Where did our bike paths go? The order you plot in matters! Switch the two
# geom_sf()
ggplot() +
  geom_sf(data = MA_counties) +
  geom_sf(data = bike_trails, col = "red")

# Better! But what are these counties? Let's compute each counties' centroid
MA_counties <- MA_counties %>%
  # Add centroids to each region using purrr package. I have no idea how 
  # this works! 
  mutate(
    lon = purrr::map_dbl(geometry, ~sf::st_centroid(.x)[[1]]),
    lat = purrr::map_dbl(geometry, ~sf::st_centroid(.x)[[2]])
    )

# Plot!
ggplot() +
  geom_sf(data = MA_counties, size = 0.5) +
  geom_sf(data = bike_trails, col = "red") +
  geom_point(data = MA_counties, aes(x=lon, y=lat)) 

# Can we do better? Yes! Use geom_text() instead with the label aesthetic!
ggplot() +
  geom_sf(data = MA_counties, size = 0.5) +
  geom_sf(data = bike_trails, col = "red") +
  geom_text(data = MA_counties, aes(x=lon, y=lat, label = name)) 



# Lec 29 Code: Wednesday 2018/11/14
library(tidyverse)
library(tidycensus)
library(sf)
library(stringr)
# Get MA county sf object, that includes population estimate column in data
mass_pop_orig <- 
  tidycensus::get_acs(
    geography = "county", 
    variables = "B01003_001", 
    state = "MA",
    geometry = TRUE
  ) %>% 
  # Add centroids to each region using purrr package (I have no idea how this 
  # works!)
  mutate(
    lon = purrr::map_dbl(geometry, ~sf::st_centroid(.x)[[1]]),
    lat = purrr::map_dbl(geometry, ~sf::st_centroid(.x)[[2]])
  )  
mass_pop_orig

# Unfortunately the original NAME variable above is very long and redundant.
# Create a new data frame with the just the county name using the separate()
# function from the tidyr package. 
mass_pop <- mass_pop_orig %>% 
  tidyr::separate(NAME, c("County", "State"), sep = ", ") %>% 
  tidyr::separate(County, c("County", "Fluff"), sep = " ") 

# Plot gradient choropleth map based on numerical variable estimate
ggplot() +
  geom_sf(data = mass_pop, aes(fill = estimate), size = 0.5) +
  geom_text(data = mass_pop, aes(x=lon, y=lat, label = County)) +
  labs(fill = "Population\nEstimate")

# What if instead we want a bin based choropleth map. Use the cut_number()
# function to convert the numerical variable estimate to a categorical 
# variable with 3 levels!
mass_pop <- mass_pop %>%
  mutate(estimate_discrete = cut_number(estimate, 3))
mass_pop

# Plot!
ggplot() +
  geom_sf(data = mass_pop, aes(fill = estimate_discrete), size = 0.5) +
  geom_text(data = mass_pop, aes(x=lon, y=lat, label = County)) +
  labs(fill = "Population\nEstimate")

# Question: How do we change the color palette for a sequential/hierarchical
# categorical variable again?