Notes 10-25

Use the development version of ggplot2 to get access to geom_sf().

# devtools::install_github("tidyverse/ggplot2")

library(tidycensus)

## Warning: package 'tidycensus' was built under R version 3.4.2

library(tidyverse)

## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr

## Warning: package 'dplyr' was built under R version 3.4.2

## Conflicts with tidy packages ----------------------------------------------

## filter(): dplyr, stats
## lag():    dplyr, stats

library(viridis)

## Loading required package: viridisLite

library(stringr)
census_api_key("7491cf955005ef7221778ebe1228115c53a12df6")

## To install your API key for use in future sessions, run this function with `install = TRUE`.

Get a list of the acs variables.

v15 <- load_variables(2015, "acs5", cache = TRUE)

Let’s look at median household income in the last 12 months for the State of Washington using ACS 2015 5 year estimates.

options(tigris_use_cache = TRUE)

medhh <- get_acs(state = "WA", geography = "county", 
                  variables = "B19013_001", geometry = TRUE)

head(medhh)

## Simple feature collection with 6 features and 5 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -124.7631 ymin: 45.83596 xmax: -117.8194 ymax: 48.55072
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs
## # A tibble: 6 x 6
##   GEOID                        NAME   variable estimate   moe
##   <chr>                       <chr>      <chr>    <dbl> <dbl>
## 1 53005   Benton County, Washington B19013_001    60251  1370
## 2 53007   Chelan County, Washington B19013_001    51837  1886
## 3 53009  Clallam County, Washington B19013_001    47253  1606
## 4 53015  Cowlitz County, Washington B19013_001    47452  1861
## 5 53037 Kittitas County, Washington B19013_001    46458  3093
## 6 53043  Lincoln County, Washington B19013_001    46069  1791
## # ... with 1 more variables: geometry <S3: sfc_MULTIPOLYGON>

medhh %>%
  mutate(NAME = str_replace(NAME,'County, Washington','')) %>% 
  ggplot(aes(x = estimate, y = reorder(NAME, estimate))) + 
  geom_point()

medhh %>% 
  ggplot(aes(fill = estimate, color = estimate)) + 
  geom_sf() + 
  coord_sf(crs = 26911) + 
  scale_fill_viridis(option = "magma") + 
  scale_color_viridis(option = "magma")

Focus on Thurston County.

thursinc <- get_acs(state = "WA", county = "Thurston", geography = "tract", variables = "B19013_001", geometry = TRUE)

head(thursinc)

## Simple feature collection with 6 features and 5 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -123.0239 ymin: 46.89342 xmax: -122.7017 ymax: 47.168
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs
## # A tibble: 6 x 6
##         GEOID                                             NAME   variable
##         <chr>                                            <chr>      <chr>
## 1 53067010400    Census Tract 104, Thurston County, Washington B19013_001
## 2 53067010510 Census Tract 105.10, Thurston County, Washington B19013_001
## 3 53067011100    Census Tract 111, Thurston County, Washington B19013_001
## 4 53067011822 Census Tract 118.22, Thurston County, Washington B19013_001
## 5 53067012330 Census Tract 123.30, Thurston County, Washington B19013_001
## 6 53067990100   Census Tract 9901, Thurston County, Washington B19013_001
## # ... with 3 more variables: estimate <dbl>, moe <dbl>, geometry <S3:
## #   sfc_MULTIPOLYGON>

thursinc %>%
  mutate(NAME = str_replace(NAME,'Thurston County, Washington','')) %>% 
  mutate(NAME = str_replace(NAME, "Census Tract","")) %>% 
  ggplot(aes(x = estimate, y = reorder(NAME, estimate))) + 
  geom_point()

## Warning: Removed 1 rows containing missing values (geom_point).

thursinc %>% 
  ggplot(aes(fill = estimate, color = estimate)) + 
  geom_sf() + 
  coord_sf(crs = 26911) + 
  scale_fill_viridis(option = "magma") + 
  scale_color_viridis(option = "magma")

Let’s redo this using Silge’s code.

library(leaflet)
library(stringr)
library(sf)

## Linking to GEOS 3.6.1, GDAL 2.1.3, proj.4 4.9.3

pal <- colorQuantile(palette = "viridis", domain = thursinc$estimate, n = 10)

thursinc %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "Esri.WorldStreetMap") %>%
    addPolygons(popup = ~ str_extract(NAME, "^([^,]*)"),
                stroke = FALSE,
                smoothFactor = 0,
                fillOpacity = 0.4,
                color = ~ pal(estimate)) %>%
    addLegend("bottomright", 
              pal = pal, 
              values = ~ estimate,
              title = "Income percentiles",
              opacity = 1)