# Automating census data pulls with r

I had a pretty cool little quantitative micro-targeting of policy issues model I wanted to show you guys today. But I had to get this Census Bureau API-R relationship smoothed out for a work thing and I’m kinda thinking it will have more universal appeal. So I’m posting it first.

Here is the gist: The Census Bureau collects lots of really cool data: demographic, socio-economic, social structure type stuff:

• educational attainment by state, congressional district, city, census tract, ext
• percent of population speaking a language other than English at home…by state, census tract, school district, whatever
• number of homeowners versus renters by various geographies

I know everybody is all beside themselves about using Twitter/Facebook/whatever for some micro-analysis. But we shouldn’t forget that for a lot of applications (political issue tracking, environmental justice, economic development) we can get a lot of milage out of the stuff the Census Bureau already collects.

The rub is that the Census Bureau’s data organization is kinda poor and their web GUIs for data retrieval are really confusing and hard to use.

Most of the stuff in this post will rely on API calls rather than static data files…but there are a few .csv files you may want and the full annotated R code is also available here:

# Pulling Census Data in R

Have you ever wanted to know how many people in your city speak Taglog at home? Or how many black people live in Vermont? Or what percent of the people in Jackson, Mississippi have a college degree? Then the American Fact Finder web GUI is for you.

But if you’ve ever wondered how the age structure of Franklin County, TN has changed over the last 7 years (perhaps as job opportunities in rural areas continue to decline and young people increasingly migrate to urban areas) AND you’d like to compare that change to changes in other rural geographies in the U.S. AND you’d like your analysis to be reproducible AND possibly applicable to other similar research questions…then you owe it to yourself to get friendly with the Census Bureau’s API.

Here’s a great resources that I borrowed from pretty heavily

## Prerequisites

There are two crucial prerequisites for this section:

1. An API key from the Census Bureau…get one here

2. The R package RJSONIO.

A quick word on item 2: there is an R package called acs or something like that which is designed to work with data from the Census Bureau’s American Communities Survey. Try it out if you want. I monkeyed with it but found it was overkill. The fromJSON() function in the RJSONIO library is pretty much all I need…or want for that matter

## An Example

The sample API call on the Census Bureau’s Developers page to get population by state from the 1 year estimates of the 2015 American Communities Survey is:

api.census.gov/data/2015/acs1?get=NAME,B01001_001E&for=state:*&key=…

In the code below I pass this API call to the Census Bureau, use the fromJSON() function to collect the JSON results in an R list, and clean up the list

library(RJSONIO)
key <- mykey
resURL <- paste('http://api.census.gov/data/2015/acs1?get=NAME,B01001_001E&for=state:*&key=',
+                 key,sep="")
ljson<-fromJSON(resURL)
[[1]]
[1] "NAME"        "B01001_001E" "state"

[[2]]
[1] "Alabama" "4858979" "01"

[[3]]

[[4]]
[1] "Arizona" "6828065" "04"

[[5]]
[1] "Arkansas" "2978204"  "05"

[[6]]
[1] "California" "39144818"   "06"



I can drop this into a dataframe pretty easily with:

ljson<-ljson[2:length(ljson)]
pop <- data.frame(statename=unlist(lapply(ljson,function(x){x[1]})),
totalpop=unlist(lapply(ljson,function(x){x[2]})),
state=unlist(lapply(ljson,function(x){x[3]})))
statename totalpop state
1    Alabama  4858979    01
3    Arizona  6828065    04
4   Arkansas  2978204    05
5 California 39144818    06


## Another Example

I like examples. Especially when you are coding, I think it’s hyper-useful to have a stock of examples. It ups the probability that, for whatever specific thing you want to do, you’ll just have to alter a couple lines of existing code.

Here’s what’s going on in this example:

• I used the API call to pull total Hispanic population by congressional district for the U.S.
• used a separate API call to pull total population by congressional district for the U.S.
• combined those 2 to create the % Hispanic population in each congressional district
• merged that with results of the 2016 Congressional Elections
#set up the API call for Hispanic Population by Congressional District
resURL <- paste("http://api.census.gov/data/2015/acs1?get=NAME,B01001I_001E&for=congressional+district:*&key=",
key,sep="")
latinos <- fromJSON(resURL)
latinos <- latinos[2:length(latinos)]

#Set up the API call for total population
resURL <- paste("http://api.census.gov/data/2015/acs1?get=NAME,B01001_001E&for=congressional+district:*&key=",
key,sep="")
pop <- fromJSON(resURL)
pop <- pop[2:length(pop)]
pop.L <- data.frame(name=unlist(lapply(pop,function(x){unlist(x)[1]})),
pop=unlist(lapply(pop,function(x){unlist(x)[2]})),
state=unlist(lapply(pop,function(x){unlist(x)[3]})),
cd=unlist(lapply(pop,function(x){unlist(x)[4]})))

latinos <- data.frame(name=unlist(lapply(latinos,function(x){unlist(x)[1]})),
latino.pop=unlist(lapply(latinos,function(x){unlist(x)[2]})))

pop.L <- pop.L %>% inner_join(latinos,by=c('name')) %>%
mutate(pct_latino=as.numeric(as.character(latino.pop))/
as.numeric(as.character(pop)))

#get the 2016 Congressional Election results
#get the numeric codes for each state
names(state.codes) <- c('statename','state.num','state')

inner_join(state.codes,by=c('state'))

#fix pop.L congressional districts
pop.L$cd <- as.numeric(as.character(pop.L$cd))
pop.L$state <- as.numeric(as.character(pop.L$state))

#order it by latino population and see if there are any notable 'R's

state cd party                                                     name    pop latino.pop pct_latino
1      6 21     R   Congressional District 21 (114th Congress), California 716371     535377 0.74734600
2     48 23     R        Congressional District 23 (114th Congress), Texas 747732     521424 0.69734076
3     12 26     R      Congressional District 26 (114th Congress), Florida 776959     541186 0.69654383
4      6 10     R   Congressional District 10 (114th Congress), California 739784     316229 0.42746126
5      6 25     R   Congressional District 25 (114th Congress), California 720316     268191 0.37232409
6      6 24     D   Congressional District 24 (114th Congress), California 736757     264994 0.35967626
7     32  4     D        Congressional District 4 (114th Congress), Nevada 722406     201492 0.27891795
8      6 49     R   Congressional District 49 (114th Congress), California 735828     189640 0.25772327
9      8  3     R      Congressional District 3 (114th Congress), Colorado 737812     181423 0.24589326
10     4  1     D       Congressional District 1 (114th Congress), Arizona 759663     173716 0.22867508
11    17 10     D     Congressional District 10 (114th Congress), Illinois 714251     160421 0.22460032
12    12  7     D       Congressional District 7 (114th Congress), Florida 738367     157720 0.21360651
13     8  6     R      Congressional District 6 (114th Congress), Colorado 796156     163880 0.20583906
14    32  3     D        Congressional District 3 (114th Congress), Nevada 758676     125431 0.16532881
15     6  7     D    Congressional District 7 (114th Congress), California 739069     117592 0.15910828
16     6 52     D   Congressional District 52 (114th Congress), California 755498     118641 0.15703682
17    12 18     R      Congressional District 18 (114th Congress), Florida 748028     111507 0.14906795
18    36  1     R      Congressional District 1 (114th Congress), New York 726589     107873 0.14846495
19    51 10     R     Congressional District 10 (114th Congress), Virginia 807670     105618 0.13076875
20    34  5     D    Congressional District 5 (114th Congress), New Jersey 743424      94463 0.12706477
21    20  3     R        Congressional District 3 (114th Congress), Kansas 755396      88598 0.11728683
22    31  2     R      Congressional District 2 (114th Congress), Nebraska 652870      71614 0.10969106
23    36  3     D      Congressional District 3 (114th Congress), New York 721764      74127 0.10270255
24    12 13     D      Congressional District 13 (114th Congress), Florida 716429      72211 0.10079296
25    36 19     R     Congressional District 19 (114th Congress), New York 705044      51710 0.07334294
26    19  3     R          Congressional District 3 (114th Congress), Iowa 812464      54484 0.06706020
27    27  2     R     Congressional District 2 (114th Congress), Minnesota 695029      40193 0.05782924
28    55  8     R     Congressional District 8 (114th Congress), Wisconsin 727491      36678 0.05041712
29    42  8     R  Congressional District 8 (114th Congress), Pennsylvania 707083      35083 0.04961652
30    26  8     R      Congressional District 8 (114th Congress), Michigan 728781      35125 0.04819692
31    36 24     R     Congressional District 24 (114th Congress), New York 708959      31211 0.04402370
32    27  3     R     Congressional District 3 (114th Congress), Minnesota 700121      30328 0.04331823
33    26  7     R      Congressional District 7 (114th Congress), Michigan 697627      29680 0.04254422
34    19  1     R          Congressional District 1 (114th Congress), Iowa 771888      28650 0.03711678
35    36 23     R     Congressional District 23 (114th Congress), New York 708372      26115 0.03686622
36    36 22     R     Congressional District 22 (114th Congress), New York 711391      26099 0.03668728
37    33  1     D Congressional District 1 (114th Congress), New Hampshire 671640      24624 0.03666250
38    18  9     R       Congressional District 9 (114th Congress), Indiana 741982      24680 0.03326226
39    36 21     R     Congressional District 21 (114th Congress), New York 710842      23524 0.03309315
40    17 12     R     Congressional District 12 (114th Congress), Illinois 699369      23119 0.03305694
41    26  1     R      Congressional District 1 (114th Congress), Michigan 700136      13902 0.01985614
42    27  8     D     Congressional District 8 (114th Congress), Minnesota 663670      11133 0.01677490
43    23  2     R         Congressional District 2 (114th Congress), Maine 657319       8938 0.01359766
>


I didn’t really do this for any reason other than it sounded fun…but it does seem pretty interesting that the top congressional districts ranked by latino population are represented by Republicans…Although it deserves mention that I didn’t do any voting age population filters or anything like that here. It’s pretty off-the-cuff. But interesting nonetheless.

## The Shitty Stuff

After toggling my way through the Census Bureau’s data GUI a few times to try and find basic shit like population, educational attainment, etc., I loved working with the API. However, the drawback to automating data retrieval with the API is that you have to know what series you are looking for. The Census Bureau does not make this easy.

Here is an example of where the API is still cool…but maybe not SUPER cool:

Suppose I want population by age and sex for each state in the U.S. and I start with the American Fact Finder. I can use either the ‘guided search’ or ‘advanced search’ to work my way through several filtering screens and arrive at an html table that will have a row/column (depending on how I decide to format the web query output) for each state and a row/column for each of like 100 variables (total population, total male population, males ages 5 and under, males ages 5-17, males ages 18-19,….total females, females ages 5 and under, etc., etc.). I can choose to download this table into a .csv format.

Alternatively, I can try to save myself several points and clicks and get the same data using the API. Here’s the rub: when I used the web tool I discovered that this particular table was Table B01001. THERE IS NO API CALL FOR TABLE B01001. There is an API call for:

• B01001_002E - total male population
• B01001_007E - total males ages 18 and 19
• B01001_026E - total female population
• B01001_040E - females age 50-54

Basically, using the API requires you to not only know the precise table you want but also the precise data series within that table.

Here is how I got 2012, 2014, and 2015 population by age and sex for each congressional district disregarding ages under 18 and over 55:

series.males <- c('001E','002E','007E','008E','009E','010E','011E','012E','013E','014E','015E','016E')
series.females <- c('026E','031E','032E','033E','034E','035E','036E','037E','038E','039E','040E')
series <- c(series.males,series.females)

series <- paste('B01001_',series,sep="")
series.names<- c('total pop','total male','m18_19','m20','m21','m22_24','m25_29','m30_34',
'm35_39','m40_44','m45_49','m50_54','total female','f18_19','f20','f21','f22_24',
'f25_29','f30_34',
'f35_39','f40_44','f45_49','f50_54')

pop.fn <- function(i,yr){
resURL <- paste('http://api.census.gov/data/',yr,
'/acs1?get=NAME,',
series[i],'&for=congressional+district:*&key=',key,sep="")
ljson <- fromJSON(resURL)
ljson <- ljson[2:length(ljson)]
tmp <- data.frame(unlist(lapply(ljson,function(x)x[1])),
unlist(lapply(ljson,function(x)x[2])),
unlist(lapply(ljson,function(x)x[3])),
unlist(lapply(ljson,function(x)x[4])),
series.names[i])
names(tmp) <- c('name','variable','state','congressional district','series_name')
return(tmp)
}

pop.df2015 <- data.frame(rbindlist(lapply(c(1:length(series)),pop.fn,yr=2015)))
pop.df2015$source <- 'ACS 2015 1 yr' pop.df2014 <- data.frame(rbindlist(lapply(c(1:length(series)),pop.fn,yr=2014))) pop.df2014$source <- 'ACS 2014 1 yr'
pop.df2012 <- data.frame(rbindlist(lapply(c(1:length(series)),pop.fn,yr=2012)))
pop.df2012$source <- 'ACS 2012 1 yr'  I don’t have much wisdom to offer on how to efficiently search for the precise series you are looking for (e.g. B01001_032E). I will tell you what I have found to be wildly inefficient: opening the variables list html table and doing a command+f search for some keyword. There are over 62,000 data series in the summary tables and another ~60,000 in the subject tables and my chrome browser hangs up for minutes trying to search the html table. I haven’t totally flushed this out yet but I was able to use R’s XML package to read the variable list html table into R. ################################################################## #Try to parse out the variables list so I can do some manner # of informed search for searies I'm looking for: library(XML) theurl <- "http://api.census.gov/data/2015/acs1/variables.html" tables <- readHTMLTable(theurl) #################################################################  From here there is an object called Concept. This field has keywords that can help you narrow your search for the right series. At this point I’m not entirely sure how to access this object…but I do know that Concept has keywords like ‘Educational Attainment’ which can, at least in theory, help you locate the stuff you want. ## A Mapping Example Now, I want to combine some of the stuff I showed off in my last post. Specifically, I’m going to use the API to pull educational attainment for a chunk of western states at the geographic unit of congressional districts and crank out a map. ################################################################ ################################################################ ################################################################ # A quick mapping example library(rgdal) library(ggplot2) library(ggmap) library(scales) library(maptools) library(rgeos) library(dplyr) #use the function I wrote above to pull total population and population with a # bachelors degree by congressional district from Census Bureau edu.fn <- function(yr){ resURL <- paste('http://api.census.gov/data/',yr, '/acs1/subject?get=NAME,S1501_C01_006E&for=congressional+district:*&', 'key=f5a32f694a14b28acf7301f4972eaab8551eafda',sep="") ljson <- fromJSON(resURL) name <- unlist(lapply(ljson[2:length(ljson)],function(x)x[1])) pop25 <- unlist(lapply(ljson[2:length(ljson)],function(x)x[2])) state <- unlist(lapply(ljson[2:length(ljson)],function(x)x[3])) congressional_district <- unlist(lapply(ljson[2:length(ljson)],function(x)x[4])) df <- data.frame(name=name,pop25=value,state=state,congressional_district=congressional_district, source=paste('ACS_1yr_',yr,sep="")) #add in the number of people 25 and over with a bachelor's degree resURL <- paste('http://api.census.gov/data/',yr, '/acs1/subject?get=NAME,S1501_C01_012E&for=congressional+district:*&', 'key=f5a32f694a14b28acf7301f4972eaab8551eafda',sep="") ljson <- fromJSON(resURL) pop25_bachelors <- unlist(lapply(ljson[2:length(ljson)],function(x)x[2])) df$pop25_bachelors <- pop25_bachelors
return(df)
}

#get educational attainment from the 2015 ACS
edu2015 <- edu.fn(yr=2015)

#read the shapefile with cartographic boundaries for congressional districts
#let's see if we can map some shit by congressional district
layer = "cb_2015_us_cd114_500k")
#pull out state, CD, and AFFGEOID so we can merge with the education stuff
df.tmp <- congress@data[,c('STATEFP','CD114FP','AFFGEOID')]
df.tmp$state_cd <- paste(df.tmp$STATEFP,df.tmp$CD114FP,sep="") congress <- fortify(congress, region="AFFGEOID") #keep edu from some western states (CA, AZ, NV, UT) edu.tmp <- edu2015 %>% mutate(state_cd=paste(state,congressional_district,sep="")) %>% inner_join(df.tmp,by=c('state_cd')) %>% filter(state %in% c('06','04','49','32')) %>% rename(id=AFFGEOID) %>% mutate(pct_college=as.numeric(as.character(pop25_bachelors))/ as.numeric(as.character(pop25))) plotData <- inner_join(congress,edu.tmp) #trim the congress data frame ids <- unique(edu.tmp$id)
congress <- congress %>% filter(congress$id %in% ids) ggplot() + geom_polygon(data = plotData, aes(x = long, y = lat, group = group, fill = 1-pct_college)) + geom_polygon(data = congress, aes(x = long, y = lat, group = group), fill = NA, color = "black", size = 0.25) + coord_map() + scale_fill_distiller(palette = "YlOrRd", labels = percent, breaks = pretty_breaks(n = 15)) + guides(fill = guide_legend(reverse = FALSE)) + theme_bw()  Ok, so the map is a little campy and maybe not super informative without some city-type context. It’s also probably worth noting that here I plotted 1 - percent of age 25 and over population wiht a Bachelor’s Degree. I wanted to display the education ‘hotspots’ in the darkest shades but didn’t want to spend all day monkeying with the palette brewer…so the darker reds are the congressional districts where relatively more of the population holds Bachelor’s Degrees. I’m pretty sure there is some really cool quantitive policy I could do with this ability to quickly pull data Census data at the Congressional District level with the API and map it in R. Hopefully, I’ll be able to illustrate some of the analysis I have in mind in the next few blog posts. I’ll try one more with a static google map basemap just to see what it looks like: ## One last map…I promise This one comes directly from Kevin Johnson’s blog post on clipping static Google maps to Census Bureau boundaries to create basemaps. Here I’m reverting back to the Census Tract boundaries and the percent of population with health insurance for California that I used in my last post. I’ll work on amending this to work with the Congressional District Boundaries later today. Note that I’m using a pretty tight zoom around Fresno California. That’s because the larger the extent of the map, the less you will see of the features on the underlying basemap….and if the map is so big that you can’t see the word Fresno it doesn’t really make a ton of sense to add a context-rich basemap at all. ################################################################ ################################################################ ################################################################ #try it with a base map #load the raster library library(raster) #go back to California Census Tract data because the map tract <- readOGR(dsn = "data/gz_2010_06_140_00_500k", layer = "gz_2010_06_140_00_500k") #load back in the insurance data we pulled ins <- read.csv("data/CA_insured_tract.csv") ins$id <- as.character(ins$GEO.id) ins$percent <- ins$Insured18_64/ins$Pop18_64

#get a basemap
map <- get_map("Fresno", zoom = 11, maptype = "roadmap")
p <- ggmap(map)
p

#establish a bounding box around our static google map
box <- as(extent(as.numeric(attr(map, 'bb'))[c(2,4,1,3)] +
c(.001,-.001,.001,-.001)), "SpatialPolygons")
proj4string(box) <- CRS(summary(tract)[[4]])
tractSub <- gIntersection(tract, box, byid = TRUE,
id = as.character(tract\$GEO_ID))
tractSub <- fortify(tractSub, region = "GEO_ID")
plotData <- left_join(tractSub, ins, by = "id")

ggmap(map) +
geom_polygon(data = plotData, aes(x = long, y = lat, group = group,
fill = percent), colour = NA, alpha = 0.5) +
scale_fill_distiller(palette = "YlOrRd", breaks = pretty_breaks(n = 10),
labels = percent) +
labs(fill = "") +
theme_nothing(legend = TRUE) +
guides(fill = guide_legend(reverse = TRUE, override.aes =
list(alpha = 1)))


Written on February 17, 2017