Step 1 - Set up the packages

Package(s): dplyr, tigris, acs, stringr, leaflet, maptools

knitr::opts_chunk$set(collapse=TRUE)

require(dplyr)
require(maptools) # required for rgdal to work correctly
require(tigris)
require(acs)
require(stringr) # to pad FIPS codes
require(leaflet)

Step 2 - Get the spatial data

Package(s): tigris

Loads a SpatialPolygonsDataFrame from the US Census website. Although I believe that the actual FIPS codes are 5 digits (the first two digits are the state, followed by three for the county), tigris seems to only work by removing the state digits.

## Grab the spatial data for all counties
tracts <- tracts(state = 'TX', cb=TRUE)
## Warning in readOGR(dsn = cache_dir, layer = shape, encoding = "UTF-8",
## verbose = FALSE, : Z-dimension discarded
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
# Note: You can use county names in the tigris 
# package but not in the acs.fetch function. 
# For simplicity, I use the FIPS ID for both

Step 3 - Get the tabular data

Package(s): acs

Pull the data that we’re investigating from the US Census website. This requires an API key, and you’ll need to know the specific table that contains the variable you’re investigating. In this example, we’ll explore income level

If needed, get an API key from the Census website in order to use the acs package

# ## Create a geographic set to grab tabular data (acs)
# geo.tab <- geo.make(
#   state=c("TX"),   # List of states
#   county="Travis", # Here's our list of counties
#   tract="*")       # All tracts

## Fetch the data from the Census website
fetched <- acs.fetch(
  geography = geo.make(state = "TX", county="*", tract = "*"),
  endyear = 2012, span = 5,# Package only goes to 2013, so end=2012
  table.number = "B19001", # Table showing 'Income'
  col.names = "pretty")    # Gives the full column definitions


## The resulting "fetched" object is not a data.frame; it's a list
names(attributes(fetched))    # see what's available
##  [1] "endyear"        "span"           "acs.units"      "currency.year" 
##  [5] "modified"       "geography"      "acs.colnames"   "estimate"      
##  [9] "standard.error" "class"
attr(fetched, "acs.colnames") # see column names
##  [1] "Household Income: Total:"              
##  [2] "Household Income: Less than $10,000"   
##  [3] "Household Income: $10,000 to $14,999"  
##  [4] "Household Income: $15,000 to $19,999"  
##  [5] "Household Income: $20,000 to $24,999"  
##  [6] "Household Income: $25,000 to $29,999"  
##  [7] "Household Income: $30,000 to $34,999"  
##  [8] "Household Income: $35,000 to $39,999"  
##  [9] "Household Income: $40,000 to $44,999"  
## [10] "Household Income: $45,000 to $49,999"  
## [11] "Household Income: $50,000 to $59,999"  
## [12] "Household Income: $60,000 to $74,999"  
## [13] "Household Income: $75,000 to $99,999"  
## [14] "Household Income: $100,000 to $124,999"
## [15] "Household Income: $125,000 to $149,999"
## [16] "Household Income: $150,000 to $199,999"
## [17] "Household Income: $200,000 or more"

## Convert to a data.frame for merging
acs_df <- data.frame(
  paste0(
    str_pad(fetched@geography$state,  2, "left", pad="0"), 
    str_pad(fetched@geography$county, 3, "left", pad="0"), 
    str_pad(fetched@geography$tract,  6, "left", pad="0")),
  fetched@estimate[,c("Household Income: Total:", "Household Income: $200,000 or more")], 
  stringsAsFactors = FALSE)

acs_df           <- select(acs_df, 1:3) %>% tbl_df()
rownames(acs_df) <- 1:nrow(acs_df)
names(acs_df)    <- c("GEOID", "total", "over_200")
acs_df$percent   <- 100*(acs_df$over_200/acs_df$total)

Step 4 - Do the merge

Package(s): tigris

df_merged <- geo_join(tracts, acs_df, "GEOID", "GEOID")

# there are some tracts with no land that we should exclude
df_merged <- df_merged[df_merged$ALAND>0,]

Step 5 - Make the map

Package(s): leaflet

popup <- paste0("GEOID: ", df_merged$GEOID, "<br>", "Percent of Households above $200k: ", round(df_merged$percent,2))
pal <- colorNumeric(
  palette = "YlGnBu",
  domain = df_merged$percent
)
leaflet() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(data = df_merged, 
              fillColor = ~pal(percent), 
              color = "#b2aeae", # you need to use hex colors
              fillOpacity = 0.7, 
              weight = 1, 
              smoothFactor = 0.2,
              popup = popup) %>%
  addLegend(pal = pal, 
            values = df_merged$percent, 
            position = "bottomright", 
            title = "Percent of Households<br>above $200k",
            labFormat = labelFormat(suffix = "%"))