Author: Elhakim Ibrahim

Instructor: Corey Sparks, PhD

February 10, 2020

Objective of the exercise

The focus of this exercise is to create basic GIS map to show spatial distribution of median household income in Bexar County, TX census tracts. To facilitate this exercise, the 2017 5-year American Community Survey census-tracts level data for Bexar County were accessed from US Census Bureau (https://data.census.gov/).

Set working directory and load required packages

setwd("C:/PhD/YR1/02SP20/DEM7093/HWK/01/OUTPUT")

suppressPackageStartupMessages({
    library(sf)
    library(ggsn)
    library(dplyr)
    library(mapview)
    library(ggplot2)
    library(classInt)
    library(patchwork)
    library(tidyverse)
    library(tidycensus)
    library(RColorBrewer)
})

Data preparation

Load data file

# demographic profile tables
v15_profile <- load_variables(2015, "acs5/profile", cache = T)

# View(v15_profile)

# search for variables by keywords in the label
v15_profile[grep(x = v15_profile$label, "Median"), c("name", "label")]
v15_profile[grep(x = v15_profile$label, "Population"), c("name", "label")]

Extract 2013-2017 5-year ACS summary file data on tottal population (DP05_0001E) and median income (DP03_0062E)

Option geometry=T extracts spatial elements of the data file, while option output="wide" arranges each variable in a column of the data set with each row being a census tract.

sa_acsmedinc2 <- get_acs(geography = "tract", state = "TX", county = c("Bexar"), 
    year = 2017, variables = c("DP05_0001E", "DP03_0062E"), geometry = T, output = "wide")
Getting data from the 2013-2017 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

Create 5-digit county FIPS code, rename variables and filter missing cases

sa_acsmedinc2$county <- substr(sa_acsmedinc2$GEOID, 1, 5)

sa_acsmedinc <- sa_acsmedinc2 %>% mutate(totpop = DP05_0001E, medinc = DP03_0062E/1000) %>% 
    # median income is divided by 1000 to avoid ranges of scientific numbers in
# map legend
na.omit()

Write data out to shapefile, change the directory

sf::st_write(sa_acsmedinc, dsn = "C:/PhD/YR1/02SP20/DEM7093/HWK/01/DATA/lab01data01r", 
    layer = "sa_tract_dp03_R", driver = "ESRI Shapefile", delete_layer = T, 
    update = T)
Deleting layer `sa_tract_dp03_R' using driver `ESRI Shapefile'
Writing layer `sa_tract_dp03_R' to data source `C:/PhD/YR1/02SP20/DEM7093/HWK/01/DATA/lab01data01r' using driver `ESRI Shapefile'
features:       362
fields:         9
geometry type:  Multi Polygon

Basic mapping of spatial data

Creating non-interactive maps

Following steps generate a Quantile break for median household income by census tracts and compare it to a Jenks break

# Population in poverty:
medinc_map <- sa_acsmedinc %>% mutate(qmedinc = cut(medinc, breaks = quantile(medinc, 
    na.rm = T, p = seq(0, 1, length.out = 6)), include.lowest = T), jmedinc = cut(medinc, 
    breaks = data.frame(classIntervals(var = sa_acsmedinc$medinc, n = 5, style = "pretty")[2])[, 
        1], include.lowest = T))

# Quantile break map parameters:
map1 <- ggplot(medinc_map, aes(fill = qmedinc, color = qmedinc)) + geom_sf() + 
    ggtitle("Spatial distribution of median household income", subtitle = "Bexar County, Texas, 2017 ACS [Quantile Breaks]") + 
    scale_fill_brewer(palette = "YlOrRd") + scale_color_brewer(palette = "YlOrRd") + 
    theme(axis.text.x = element_blank(), axis.text.y = element_blank()) + north(medinc_map) + 
    scalebar(medinc_map, location = "bottomleft", dist = 5, transform = T, dist_unit = "km", 
        model = "WGS84", st.size = 2)

# Jenks break map parameters:
map2 <- ggplot(medinc_map, aes(fill = jmedinc, color = jmedinc)) + geom_sf() + 
    ggtitle("Spatial distribution of median household income", subtitle = "Bexar County, Texas, 2017 ACS [Jenks Breaks]") + 
    scale_fill_brewer(palette = "YlOrRd") + scale_color_brewer(palette = "YlOrRd") + 
    theme(axis.text.x = element_blank(), axis.text.y = element_blank()) + north(medinc_map) + 
    scalebar(medinc_map, location = "bottomleft", dist = 5, transform = T, dist_unit = "km", 
        model = "WGS84", st.size = 2)

Meanwhile, if I wish to choose different colour depending on the nature of my data, I could first apply display.brewer.all() function that renders all colour palette options in RColorBrewer package. The colours are split into three groups; namely, sequential, qualitative and diverging.

display.brewer.all()

ggsave(filename = "C:/PhD/YR1/02SP20/DEM7093/HWK/01/OUTPUT/brewercolourpalette.png")
Saving 8 x 12 in image

Then, I output separate Quantile break and Jenks break maps and save.

map1

ggsave(filename = "C:/PhD/YR1/02SP20/DEM7093/HWK/01/OUTPUT/hmk01map01.png")
Saving 7 x 5 in image
map2

ggsave(filename = "C:/PhD/YR1/02SP20/DEM7093/HWK/01/OUTPUT/hmk01map02.png")
Saving 7 x 5 in image

Next, I generate stacked Quantile break and Jenks break maps and save

map1/map2

ggsave(filename = "C:/PhD/YR1/02SP20/DEM7093/HWK/01/OUTPUT/hmk01map03.png")
Saving 14 x 10 in image

Creating interactive map with mapview

I generate my interactive map based on Quantile break function. Also, I set color pallette for my map using brewer.pal function together with colorRampPalatte command which offers more flexibility for extending scale of colours. Assuming, for instance, that the variable to map contains 30 levels, parsing colorRampPalette(brewer.pal(9,“Blues”))(30) line of code will generate 30 colours, one for each level, based on the 9 from the “Blues” colour palette.

mapview(medinc_map["qmedinc"], col.regions = colorRampPalette(brewer.pal(6, 
    "YlOrRd")), legend = T, map.types = "OpenStreetMap", layer.name = "Median Income")
ggsave(filename = "C:/PhD/YR1/02SP20/DEM7093/HWK/01/OUTPUT/hwk01map04.png")
Saving 7 x 5 in image
############################################################################################################################