DEM 5093/7093 Lab 1 using R - Basic Map making
Jacob M. Souch.
2023-02-01
library(tidycensus)
library(tidyverse)
library(sf)
library(ggplot2)
library(classInt)
library(dplyr)setwd("C:/Users/jacob/OneDrive - University of Texas at San Antonio/Courses/GIS for Population Science")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 = "INCOME", x = label))%>%
select(name, label)## # A tibble: 174 × 2
## name label
## <chr> <chr>
## 1 DP03_0051 Estimate!!INCOME AND BENEFITS (IN 2019 INFLATION-ADJUSTED DOLLARS…
## 2 DP03_0051P Percent!!INCOME AND BENEFITS (IN 2019 INFLATION-ADJUSTED DOLLARS)…
## 3 DP03_0052 Estimate!!INCOME AND BENEFITS (IN 2019 INFLATION-ADJUSTED DOLLARS…
## 4 DP03_0052P Percent!!INCOME AND BENEFITS (IN 2019 INFLATION-ADJUSTED DOLLARS)…
## 5 DP03_0053 Estimate!!INCOME AND BENEFITS (IN 2019 INFLATION-ADJUSTED DOLLARS…
## 6 DP03_0053P Percent!!INCOME AND BENEFITS (IN 2019 INFLATION-ADJUSTED DOLLARS)…
## 7 DP03_0054 Estimate!!INCOME AND BENEFITS (IN 2019 INFLATION-ADJUSTED DOLLARS…
## 8 DP03_0054P Percent!!INCOME AND BENEFITS (IN 2019 INFLATION-ADJUSTED DOLLARS)…
## 9 DP03_0055 Estimate!!INCOME AND BENEFITS (IN 2019 INFLATION-ADJUSTED DOLLARS…
## 10 DP03_0055P Percent!!INCOME AND BENEFITS (IN 2019 INFLATION-ADJUSTED DOLLARS)…
## # … with 164 more rowsv15_Profile%>%
filter(grepl(pattern = "household", x = label))%>%
select(name, label)## # A tibble: 194 × 2
## name label
## <chr> <chr>
## 1 DP02_0001 Estimate!!HOUSEHOLDS BY TYPE!!Total households
## 2 DP02_0001P Percent!!HOUSEHOLDS BY TYPE!!Total households
## 3 DP02_0002 Estimate!!HOUSEHOLDS BY TYPE!!Total households!!Married-couple fa…
## 4 DP02_0002P Percent!!HOUSEHOLDS BY TYPE!!Total households!!Married-couple fam…
## 5 DP02_0003 Estimate!!HOUSEHOLDS BY TYPE!!Total households!!Married-couple fa…
## 6 DP02_0003P Percent!!HOUSEHOLDS BY TYPE!!Total households!!Married-couple fam…
## 7 DP02_0004 Estimate!!HOUSEHOLDS BY TYPE!!Total households!!Cohabiting couple…
## 8 DP02_0004P Percent!!HOUSEHOLDS BY TYPE!!Total households!!Cohabiting couple …
## 9 DP02_0005 Estimate!!HOUSEHOLDS BY TYPE!!Total households!!Cohabiting couple…
## 10 DP02_0005P Percent!!HOUSEHOLDS BY TYPE!!Total households!!Cohabiting couple …
## # … with 184 more rowsExtract 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("DP03_0062E","DP05_0001E") ,
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 Profilehead(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(mhhinc = DP03_0062E,
totpop = DP05_0001E) %>%
st_transform(crs = 2919)%>%
na.omit()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.
#eval=F means eval=FALSE which means this code is not running, it is here for you to view though
#st_geometry(sa_acs2)<-NULL
foreign::write.dbf(as.data.frame(sa_acs2), file="C:/Users/jacob/OneDrive - University of Texas at San Antonio/Courses/GIS for Population Science/shapefileTX48029")#change the directory for your computer
sf::st_write(sa_acs2,
dsn="C:/Users/jacob/OneDrive - University of Texas at San Antonio/Courses/GIS for Population Science/sa_tr_dp03.gpkg",
layer="sa_tr_dp03",
driver="GPKG", append=F) ## Save as geopackage format - QGIS likes thismydat<- st_read("C:/Users/jacob/OneDrive - University of Texas at San Antonio/Courses/GIS for Population Science/sa_tr_dp03.gpkg")
names(mydat)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")
hhinc_map <-sa_acs2 %>%
mutate_map_brks(mhhinc, n=6, style="quantile")%>%
mutate_map_brks(mhhinc, n=6, style="jenks")
p1<-ggplot(hhinc_map, aes(fill = mhhinc_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(hhinc_map)+
scalebar(hhinc_map, location="bottomleft",
dist=5, transform = T,
dist_unit = "km",
model="WGS84",
st.size =2 )
p1
p2<-ggplot(hhinc_map, aes(fill = mhhinc_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(hhinc_map)+
scalebar(hhinc_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="C:/Users/xee291/OneDrive - University of Texas at San Antonio/Documents/UTSA 2022-2023/GIS Class/Lab 1/Lab 1_V2/lab1map1.png",
dpi = "print")## Saving 7 x 5 in image## Warning in dir.create(dir, recursive = TRUE): cannot create dir
## 'C:\Users\xee291', reason 'Permission denied'## Warning in normalizePath(path.expand(path), winslash, mustWork):
## path[1]="C:/Users/xee291/OneDrive - University of Texas at San
## Antonio/Documents/UTSA 2022-2023/GIS Class/Lab 1/Lab 1_V2": The system cannot
## find the path specified## Warning in grDevices::dev.off(): agg could not write to the given fileSlightly 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.
Note: View this website for more style
and other tmap tricks: https://bookdown.org/nicohahn/making_maps_with_r5/docs/tmap.html
library(tmap)
library(tmaptools)
tm_shape(sa_acs2)+
tm_polygons("mhhinc",
title="Median Household Income, Bexar County, 2019",
palette="Blues",
style="pretty", n=5 )+
tm_format("World",
legend.outside=T)+ tm_layout(main.title ="Bexar County, Census Tract-Level Median Household Income, 2019 - Pretty Breaks",main.title.size = 0.8, title.position = c('center', 'top'))+
tm_scale_bar()+
tm_compass()
Interactive map with mapview
(Not necessary for Homework 2, but here to view)
library(mapview)
library(RColorBrewer)
hhinc_map$mhhinc_jbrks<-relevel(hhinc_map$mhhinc_jbrks, ref = "0-10463")
pal <- colorRampPalette(brewer.pal(7, "Blues")) #set colors
mapview(hhinc_map,
zcol="mhhinc_qbrks",
legend=T,
map.types="OpenStreetMap",
layer.name="Median Household Income")