Plotting Census Data by PUMA

A simple way to download Census data within R as well as shapefiles can be found here.

This guided example will pull demographic data to create a map of the distribution of Arabs in the state of Michigan by Public Use Microdata Area (PUMA).

To get started, you will need to obtain an API key from the US Census. To do this, go here. (This doesn’t take long).

First, install the following packages:

# install.packages(c("tigris","acs","stringr"))
library(tigris)
library(acs)
library(stringr)
library(leaflet)
library(htmlwidgets)

Second, we get the spatial data (tigris):

# grab the spatial data (tigris)
pumas <- pumas(state='MI',cb=TRUE)

Different spatial data is available using different commands like: blocks, counties, county_subdivisions, places, school_districts, states, tracts, zctas.

Third, we pull tabular data (acs). This is where you will need your API key.

api.key.install(key="YOUR_API_KEY")
# create a geographic set to grab tabular data (acs)
geo <- geo.make(state=c('MI'),
                puma='*')

The ’*’ tells R to pull all available PUMAs. If we were interested in a subset of PUMAs, we can create a vector of the PUMA numbers of interest.

endyr=2015
sp=5
tabnum = "B04006"
ancestry <- acs.fetch(endyear=endyr,span=sp,geography=geo,
                      table.number=tabnum,col.names="pretty")

# use of col.names="pretty" abve gives the full column definitions
# if you want Census variable IDs use col.names="auto".

# The resulting 'ancestry' object is not a data.frame -
# it's a list to see what's available
names(attributes(ancestry))
attr(ancestry,"acs.colnames")
attr(ancestry,"geography")
attr(ancestry,"estimate")

# convert to a data.frame for merging
ancestry.df <- data.frame(ancestry@geography$NAME,
                          paste0(str_pad(ancestry@geography$state,2,"left",pad="0"),
                                 str_pad(ancestry@geography$publicusemicrodataarea,5,"left",pad="0")),
                          ancestry@estimate[,c("B04006. People Reporting Ancestry: Total:",
                                               "B04006. People Reporting Ancestry: Arab:",
                                               "B04006. People Reporting Ancestry: Arab: Egyptian",
                                               "B04006. People Reporting Ancestry: Arab: Iraqi",
                                               "B04006. People Reporting Ancestry: Arab: Jordanian",
                                               "B04006. People Reporting Ancestry: Arab: Lebanese",
                                               "B04006. People Reporting Ancestry: Arab: Moroccan",
                                               "B04006. People Reporting Ancestry: Arab: Palestinian",
                                               "B04006. People Reporting Ancestry: Arab: Syrian",
                                               "B04006. People Reporting Ancestry: Arab: Arab",
                                               "B04006. People Reporting Ancestry: Arab: Other Arab",
                                               "B04006. People Reporting Ancestry: Assyrian/Chaldean/Syriac")],
                                 stringsAsFactors=FALSE)

# do some cleaning of the data frame and create variables of interest
names(ancestry.df) <- c("name","GEOID","total","arabtotal","egyptian","iraqi","jordanian",
                        "lebanese","moroccan","palestinian","syrian","arab","otharab",
                        "chaldean")
ancestry.df$man.arabtotal<-ancestry.df$egyptian+ancestry.df$iraqi+ancestry.df$jordanian+ancestry.df$lebanese+
                           ancestry.df$moroccan+ancestry.df$palestinian+ancestry.df$syrian+ancestry.df$arab+ancestry.df$otharab
ancestry.df$allarab<-ancestry.df$egyptian+ancestry.df$iraqi+ancestry.df$jordanian+ancestry.df$lebanese+
                     ancestry.df$moroccan+ancestry.df$palestinian+ancestry.df$syrian+ancestry.df$arab+ancestry.df$otharab+
                     ancestry.df$chaldean
ancestry.df$pct.allarab <- 100*(ancestry.df$allarab/ancestry.df$total)

The final step before mapping the data: merge the data with the geography file.

# do the merge
ancestry.merged <- geo_join(pumas,ancestry.df,"GEOID10","GEOID")

Here is where we make the map. First we’ll specify the pop-up, then well create the map object, and then we’ll finally print the map.

# make the map
popup <- paste0("GEOID: ", ancestry.merged$NAME10, "<br>", "Percent of Population Arab:", round(ancestry.merged$pct.allarab,2))
pal <- colorNumeric(
  palette = "YlGnBu",
  domain = ancestry.merged$pct.allarab
)
map3<-leaflet() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(data = ancestry.merged, 
              fillColor = ~pal(pct.allarab), 
              color = "#b2aeae", # you need to use hex colors
              fillOpacity = 0.7, 
              weight = 1, 
              smoothFactor = 0.2,
              popup = popup) %>%
  addLegend(pal = pal, 
            values = ancestry.merged$pct.allarab, 
            position = "bottomright", 
            title = "Percent of Population<br>Arab",
            labFormat = labelFormat(suffix = "%")) 
map3

# don't forget to save
saveWidget(map3,file="map3.html",selfcontained = F)
# this will save the map in your directory and will allow
# you to open it in a web browser to peruse

And there you go!