1) Set up the packages

knitr::opts_chunk$set(collapse=TRUE)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(maptools) # required for rgdal to work correctly
## Warning: package 'maptools' was built under R version 3.3.2
## Loading required package: sp
## Checking rgeos availability: TRUE
library(tigris)
## Warning: package 'tigris' was built under R version 3.3.2
## 
## Attaching package: 'tigris'
## The following object is masked from 'package:graphics':
## 
##     plot
library(acs)
## Warning: package 'acs' was built under R version 3.3.2
## Loading required package: stringr
## Loading required package: plyr
## -------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## -------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## Loading required package: XML
## 
## Attaching package: 'acs'
## The following object is masked from 'package:dplyr':
## 
##     combine
## The following object is masked from 'package:base':
## 
##     apply
library(stringr) # to pad FIPS codes
library(leaflet)
## Warning: package 'leaflet' was built under R version 3.3.2

2) Get the spatial data (tigris)

## Grab the spatial data for all counties
tracts <- tracts(state = 'NY', 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

3) Get the tabular data (acs)

In order to do this you will need an API key from the US Census. Go to this site to request one (takes a minute or two) and then use the api.key.install function in the acs package to use the key.


# ## 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 = "NY", 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"
##  [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"
##  [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)
## Warning: Setting row names on a tibble is deprecated.
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,]
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
)

Step 5 - Make the map

Package(s): leaflet
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 = "%"))