How to create a tigris map with zip codes shaded based on a value?

Question

I'm trying to use the tigris package in R to make a map of Douglas County, Colorado. The goal is to show Douglas County divided by the zip codes inside of it, and to have each zip code shaded based on a value. The variables I have right now in my dataset are: Zip code, town, and lead level (continuous variable).

library(tigris)
library(sf)

douglas_zips <- zctas(cb=TRUE, starts_with= c("80108","80109", "80104", "80116", "80126", "80129", "80130", "80118", "80124", "80131", "80134", "80138", "80125", "80135"))
plot(douglas_zips)

When I do this, I get this strange plot: zip code map

Any thoughts or ideas of where to go from here?

Answer 1

The information from Spacedman is correct, but Matt is asking how to execute a data join between the douglas_zips feature and a separate table containing a list of lead levels.

The 'join' is an essential function in GIS and relational databases broadly, so it's a critical skill.

We want to take our douglas_zips feature and match all the records from the lead_levels table by the zipcode, adding the columns from lead_levels to douglas_zips. We can use merge() to accomplish that.

https://www.rdocumentation.org/packages/base/versions/3.6.2/topics/merge

NOTE: specify either 'cb=TRUE, year=2020' or 'cb=FALSE' (the default) for ZCTAs from tigris.

library(tigris)
library(sf)
library(plyr)

douglas_zips <- zctas(cb=TRUE, year=2020, starts_with=c("80108","80109", "80104", "80116", "80126", "80129", "80130", "80118", "80124", "80131", "80134", "80138", "80125", "80135"))

## made up data.frame to simulate what you described as your dataset
zipcode = c("80108","80109", "80104", "80116", "80126", "80129", "80130", "80118", "80124", "80131", "80134", "80138", "80125", "80135")
town = c("town1","town2","town3","town4","town5","town6","town7","town8","town9","town10","town11","town12","town13","town14")
lead_level = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14)
dataset = data.frame(zipcode, town, lead_level)

douglas_lead_levels <- merge(douglas_zips, dataset, by.x="ZCTA5CE20",by.y="zipcode")

The result will add two columns to douglas_zips: "town" and "lead_level".

From there, we can plot the lead_level variable as Spacedman explained:

plot(douglas_lead_levels["lead_level"])

Answer 2

That's showing a map of those zip areas for every column in the data frame. Each row of the data frame is a zip area, and each row has a few columns of info with it. By default, plotting such an object does a map for each column (up to a smallish limit). If you print the object out you'll see its data-frame nature:

> douglas_zips
Simple feature collection with 14 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -105.3987 ymin: 39.12947 xmax: -104.5877 ymax: 39.56851
Geodetic CRS:  NAD83
First 10 features:
      ZCTA5CE20     AFFGEOID20 GEOID20 NAME20 LSAD20   ALAND20 AWATER20
5633      80124 860Z200US80124   80124  80124     Z5  24303768     5755
5760      80135 860Z200US80135   80135  80135     Z5 964291604  4925227
8362      80129 860Z200US80129   80129  80129     Z5  18234617     9646
9575      80138 860Z200US80138   80138  80138     Z5 160361611   162369
12467     80118 860Z200US80118   80118  80118     Z5 339356454   295396
14329     80131 860Z200US80131   80131  80131     Z5   1436637        0
16056     80125 860Z200US80125   80125  80125     Z5 101565571  3920161
24583     80104 860Z200US80104   80104  80104     Z5 183043141        0
24641     80109 860Z200US80109   80109  80109     Z5  67311248    15254
24704     80134 860Z200US80134   80134  80134     Z5 148527309   320476

To plot just one column, do plot(douglas_zips["ZCTA5CE20"]):