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 = "yourkeyhere", install = T)

One time to install your key for use in tidycensus

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

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 = "POVERTY", x = label))%>%
  select(name, label)
## # A tibble: 38 × 2
##    name       label                                                             
##    <chr>      <chr>                                                             
##  1 DP03_0119  Estimate!!PERCENTAGE OF FAMILIES AND PEOPLE WHOSE INCOME IN THE P…
##  2 DP03_0119P Percent!!PERCENTAGE OF FAMILIES AND PEOPLE WHOSE INCOME IN THE PA…
##  3 DP03_0120  Estimate!!PERCENTAGE OF FAMILIES AND PEOPLE WHOSE INCOME IN THE P…
##  4 DP03_0120P Percent!!PERCENTAGE OF FAMILIES AND PEOPLE WHOSE INCOME IN THE PA…
##  5 DP03_0121  Estimate!!PERCENTAGE OF FAMILIES AND PEOPLE WHOSE INCOME IN THE P…
##  6 DP03_0121P Percent!!PERCENTAGE OF FAMILIES AND PEOPLE WHOSE INCOME IN THE PA…
##  7 DP03_0122  Estimate!!PERCENTAGE OF FAMILIES AND PEOPLE WHOSE INCOME IN THE P…
##  8 DP03_0122P Percent!!PERCENTAGE OF FAMILIES AND PEOPLE WHOSE INCOME IN THE PA…
##  9 DP03_0123  Estimate!!PERCENTAGE OF FAMILIES AND PEOPLE WHOSE INCOME IN THE P…
## 10 DP03_0123P Percent!!PERCENTAGE OF FAMILIES AND PEOPLE WHOSE INCOME IN THE PA…
## # … with 28 more rows
v15_Profile%>%
  filter(grepl(pattern = "Built 2000 to 2009", x = label))%>%
  select(name, label)
## # A tibble: 2 × 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 2019 for Bexar County, TX Census Tracts

Here is a query where we extract several variables from the 2019 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("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)

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)

# rename variables and filter missing cases
sa_acs2<-sa_acs %>%
  mutate(mhi = DP03_0062E) %>%
  st_transform(crs = 2919) %>%
  na.omit()
st_geometry(sa_acs2) <- NULL

foreign::write.dbf(as.data.frame(sa_acs2), file="/Users/briansurratt/Library/CloudStorage/OneDrive-UniversityofTexasatSanAntonio/DEM 7093 GIS/Homework 2/acs_lab1.dbf")

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="/Users/briansurratt/Library/CloudStorage/OneDrive-UniversityofTexasatSanAntonio/DEM 7093 GIS/Homework 2/sa_tr_dp03.gpkg",
             layer="sa_tr_dp03",
             driver="GPKG") ## Save as geopackage format - QGIS likes this
mydat<- st_read("/Users/briansurratt/Library/CloudStorage/OneDrive-UniversityofTexasatSanAntonio/DEM 7093 GIS/Homework 2/sa_tr_dp03.gpkg")
## Reading layer `sa_tr_dp03' from data source 
##   `/Users/briansurratt/Library/CloudStorage/OneDrive-UniversityofTexasatSanAntonio/DEM 7093 GIS/Homework 2/sa_tr_dp03.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 362 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 2029921 ymin: 13589620 xmax: 2249502 ymax: 13821860
## Projected CRS: NAD83(HARN) / Texas South Central (ftUS)
names(mydat)
## [1] "GEOID"      "NAME"       "DP03_0062E" "DP03_0062M" "county"    
## [6] "mhi"        "geom"

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(xplorerr)
source("https://raw.githubusercontent.com/coreysparks/Rcode/master/mutate_map_funs.R")

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

p1<-ggplot(mhi_map, aes(fill = mhi_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(mhi_map)+
  scalebar(mhi_map, location="bottomleft",
           dist=5, transform = T,
           dist_unit = "km",
           model="WGS84",
           st.size =2 )

p1

p2<-ggplot(mhi_map, aes(fill = mhi_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(mhi_map) +
  scalebar(mhi_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="/Users/briansurratt/Library/CloudStorage/OneDrive-UniversityofTexasatSanAntonio/DEM 7093 GIS/Homework 2/homework2map1.png",
       dpi = "print")
## Saving 7 x 5 in image

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("mhi",
              title="Median Household Income",
              palette="Reds",
              style="quantile", n=5 )+
  tm_format("World",
            title="San Antonio Median Household Income - Quantile Breaks",
            legend.outside=T)+
  tm_scale_bar()+
  tm_compass()

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

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

Interactive map with mapview

head(mhi_map)
## Simple feature collection with 6 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 2106567 ymin: 13681230 xmax: 2144180 ymax: 13744320
## Projected CRS: NAD83(HARN) / Texas South Central (ftUS)
##         GEOID                                      NAME DP03_0062E DP03_0062M
## 1 48029190601 Census Tract 1906.01, Bexar County, Texas      46818      14436
## 2 48029181100    Census Tract 1811, Bexar County, Texas      64930       6521
## 3 48029140700    Census Tract 1407, Bexar County, Texas      36222       6656
## 4 48029180602 Census Tract 1806.02, Bexar County, Texas      38389       7038
## 5 48029140800    Census Tract 1408, Bexar County, Texas      43299       3656
## 6 48029170200    Census Tract 1702, Bexar County, Texas      25755       7141
##                         geometry county   mhi         mhi_qbrks   mhi_jbrks
## 1 MULTIPOLYGON (((2119318 137...  48029 46818    42899.33-52230 38175-54289
## 2 MULTIPOLYGON (((2107703 137...  48029 64930    52230-66891.67 54289-73902
## 3 MULTIPOLYGON (((2136161 136...  48029 36222 35210.33-42899.33 10463-38175
## 4 MULTIPOLYGON (((2106567 137...  48029 38389 35210.33-42899.33 38175-54289
## 5 MULTIPOLYGON (((2134799 136...  48029 43299    42899.33-52230 38175-54289
## 6 MULTIPOLYGON (((2118664 136...  48029 25755    10463-35210.33 10463-38175
summary(mhi_map)
##     GEOID               NAME             DP03_0062E       DP03_0062M   
##  Length:362         Length:362         Min.   : 10463   Min.   : 2192  
##  Class :character   Class :character   1st Qu.: 38946   1st Qu.: 6365  
##  Mode  :character   Mode  :character   Median : 52230   Median : 8758  
##                                        Mean   : 60357   Mean   :10103  
##                                        3rd Qu.: 76107   3rd Qu.:12261  
##                                        Max.   :215083   Max.   :41980  
##                                                                        
##           geometry      county               mhi                     mhi_qbrks 
##  MULTIPOLYGON :362   Length:362         Min.   : 10463   0-10463          : 1  
##  epsg:2919    :  0   Class :character   1st Qu.: 38946   10463-35210.33   :60  
##  +proj=lcc ...:  0   Mode  :character   Median : 52230   35210.33-42899.33:60  
##                                         Mean   : 60357   42899.33-52230   :60  
##                                         3rd Qu.: 76107   52230-66891.67   :60  
##                                         Max.   :215083   66891.67-84893.33:60  
##                                                          84893.33-215083  :61  
##          mhi_jbrks  
##  0-10463      :  1  
##  10463-38175  : 83  
##  38175-54289  :106  
##  54289-73902  : 77  
##  73902-98311  : 58  
##  98311-131278 : 28  
##  131278-215083:  9
library(mapview)
library(RColorBrewer)
mhi_map$mhi_jbrks<-relevel(mhi_map$mhi_jbrks, ref = "0-10463" )
pal <- colorRampPalette(brewer.pal(7, "Blues")) #set colors
mapview(mhi_map,
        zcol="mhi_jbrks",
        legend=T,
        map.types="OpenStreetMap",
        layer.name="Median Household Income")