Get a Census developer API Key

Obtain one at the Census Developer website

Save your API key to your working directory

use tidycensus::census_api_key(key = "d42eebdfb8a5be15b37eb8cef2a3abc37a71f12b", install = T)

One time to install your key for use in tidycensus

library(tidycensus)
library(tidyverse)
library(sf)
library(ggplot2)
library(classInt)

Examine data profile tables

The load_variables() function will load all available variables in the ACS for a given year

v15_Profile <- load_variables(year = 2019 ,
                              dataset =  "acs5/profile",
                              cache = TRUE)

Calling View(v15_Profile) will let you interactively browse and filter the ACS variables, this is one way to search for what you’re looking for.

Search for variables by keywords in the label

v15_Profile%>%
  filter(grepl(pattern = "inc", x = label))%>%
  select(name, label)
## # A tibble: 32 x 2
##    name         label                                                           
##    <chr>        <chr>                                                           
##  1 DP02_0062    Estimate!!EDUCATIONAL ATTAINMENT!!Population 25 years and over!~
##  2 DP02_0062P   Percent!!EDUCATIONAL ATTAINMENT!!Population 25 years and over!!~
##  3 DP02PR_0062  Estimate!!EDUCATIONAL ATTAINMENT!!Population 25 years and over!~
##  4 DP02PR_0062P Percent!!EDUCATIONAL ATTAINMENT!!Population 25 years and over!!~
##  5 DP03_0049    Estimate!!CLASS OF WORKER!!Civilian employed population 16 year~
##  6 DP03_0049P   Percent!!CLASS OF WORKER!!Civilian employed population 16 years~
##  7 DP03_0062    Estimate!!INCOME AND BENEFITS (IN 2019 INFLATION-ADJUSTED DOLLA~
##  8 DP03_0062P   Percent!!INCOME AND BENEFITS (IN 2019 INFLATION-ADJUSTED DOLLAR~
##  9 DP03_0063    Estimate!!INCOME AND BENEFITS (IN 2019 INFLATION-ADJUSTED DOLLA~
## 10 DP03_0063P   Percent!!INCOME AND BENEFITS (IN 2019 INFLATION-ADJUSTED DOLLAR~
## # ... with 22 more rows
## grep searches for patterns of characters like in a word eg. "pov" to search entries with labels including the word "poverty"
v15_Profile%>%
  filter(grepl(pattern = "Built 2000 to 2009", x = label))%>%
  select(name, label)
## # A tibble: 2 x 2
##   name       label                                                              
##   <chr>      <chr>                                                              
## 1 DP04_0019  Estimate!!YEAR STRUCTURE BUILT!!Total housing units!!Built 2000 to~
## 2 DP04_0019P Percent!!YEAR STRUCTURE BUILT!!Total housing units!!Built 2000 to ~

Extract from ACS summary file data

The tidycensus package has a function get_acs() that will download data from the Census API for you automatically assuming you’ve installed your key from above

Here I get data profile variables from 2017 for Bexar County, TX Census Tracts

Here is a query where we extract several variables from the 2017 ACS for Bexar County, Texas. We can also get the spatial data by requesting geometry=TRUE.

Using output="wide" will put each variable in a column of the data set, with each row being a census tract.

sa_acs<-get_acs(geography = "tract",
                state="TX",
                county = "Bexar",
                year = 2019,
                variables=c("DP05_0001E", "DP03_0062E") ,
                geometry = T,
                output = "wide")
## Getting data from the 2015-2019 5-year ACS
## Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
## Using the ACS Data Profile
head(sa_acs)

## state FIP codes are first two numbers in the GEOID, county FIP codes are the next three digits, census tract FIP codes are the last 6 digits

Here, I create some other variables that we may need later

# create a county FIPS code - 5 digit
sa_acs$county<-substr(sa_acs$GEOID, 1, 5)
## subtr = "parts of that string"

# rename variables and filter missing cases
sa_acs2<-sa_acs%>%
  mutate(totpop= DP05_0001E,
         medincome=DP03_0062E) %>%
  st_transform(crs = 2919)%>%
  na.omit()

## drops any census tracts that are missing estimates

Write data out to shapefile

You may need to create or modify some data in R and then use it in the desktop GIS (QGIS), we can write any data from R into a variety of data formats using the sf::st_write() function.

#change the directory for your computer
sf::st_write(sa_acs2,
             dsn="../r_labs/sa_tr_dp03",
             layer="sa_tr_dp03",
             driver="GPKG") ## Save as geopackage format - QGIS likes this
## geopackage is a more compact version of a shapefile which includes 5 or so files

Some basic mapping of variables

Here I generate a quantile break for % black in census tracts and compare it to a Jenks break. Note in ggplot, the Jenks break is harder to do

library(patchwork)
library(ggsn)
#library(XploreR)
source("https://raw.githubusercontent.com/coreysparks/Rcode/master/mutate_map_funs.R")

income_map<-sa_acs2 %>%
  mutate_map_brks(medincome, n=6, style="quantile")%>%
           mutate_map_brks(medincome, n=6, style="jenks")



p1<-ggplot(income_map, aes(fill = medincome_qbrks)) + 
  geom_sf(color="grey50") + 
  ggtitle("Median Household Income", 
          subtitle = "Bexar County Texas, 2017 - Quantile Breaks")+
    scale_fill_brewer(palette = "Reds") + 
  scale_color_brewer(palette = "Reds")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank())+
  north(income_map)+
  scalebar(income_map, location="bottomleft",
           dist=5, transform = T,
           dist_unit = "km",
           model="WGS84",
           st.size =2 )
p1

p2<-ggplot(income_map, aes(fill = medincome_jbrks)) + 
  geom_sf(color="grey50") + 
  ggtitle("Median Household Income", 
          subtitle = "Bexar County Texas, 2017 - Jenks Breaks")+
  scale_fill_brewer(palette = "Reds") + 
  scale_color_brewer(palette = "Reds")+
    theme(axis.text.x = element_blank(),
          axis.text.y = element_blank())+
  north(income_map)+
  scalebar(income_map, 
           location="bottomleft",
           dist=5,
           transform = T,
           dist_unit = "km",
           model="WGS84",
           st.size =2)
p2

p1+ p2

You can save the image from above to your computer by using ggsave()

ggsave(filename="../images/lab1map1.png",
       dpi = "print")

Slightly easier way using tmap

The tmap package is an attractive alternative to using ggplot() when making maps, and makes basic cartographic principles easier.

Note style="fisher" is equivalent to style="jenks" and scales better to larger data.

library(tmap)
library(tmaptools)

tm_shape(sa_acs2)+
  tm_polygons("medincome",
              title="Median Household Income",
              palette="Blues",
              style="quantile", n=5 )+
  tm_format("World",
            title="San Antonio Income Estimates - Quantile Breaks",
            legend.outside=T)+
  tm_scale_bar()+
  tm_compass()

tm_shape(sa_acs2)+
  tm_polygons("medincome",
              title="Median Household Income",
              palette="Blues",
              style="fisher",
              n=5,
              legend.hist=T )+
  tm_format("World",
            title="San Antonio Income Estimates - Jenks Breaks",
            legend.outside=T)+
  tm_scale_bar()+
  tm_compass()

Interactive map with mapview

library(mapview)
library(RColorBrewer)
income_map$medincome_jbrks<-relevel(income_map$medincome_jbrks,ref = "0-10463" )

##ppov_map$ppov_jbrks was from the ggplot section

pal <- colorRampPalette(brewer.pal(7, "Blues")) #set colors

# to change the pivot table that populates when you click on a geo element on the map:
# sub <- ppov_map %>% 
  # select(Poverty_Rate = ppov_jbrks, ppov, Name = NAME)

mapview(income_map,
        zcol="medincome_jbrks",
        legend=T,
        map.types="OpenStreetMap",
        layer.name="Household Median Income")
## context is key to human behavior! context matters to the non-expert

## function map.shot will create an image of a map, saves the map as an html or an image