Congressional district maps

cong-dist-figs.R

# LA Times figures

library(rgdal)
library(RColorBrewer)
library(raster)
library(viridis)

# sources for spatial datasets
#   congressional districts:
#     https://www.census.gov/geo/maps-data/data/tiger-line.html
#   Global Lakes and Wetlands
#     https://www.worldwildlife.org/publications/ \
#     global-lakes-and-wetlands-database-large-lake-polygons-level-1
#   NC sub-counties
#     https://www2.census.gov/geo/tiger/TIGER2016/COUSUB/

# congressional districts
cd115 <- readOGR("tl_2016_us_cd115/tl_2016_us_cd115.shp",
                 layer = "tl_2016_us_cd115")
cd114 <- readOGR("tl_2015_us_cd114/tl_2015_us_cd114.shp",
                 layer = "tl_2015_us_cd114")
cd112 <- readOGR("tl_2012_us_cd112/tl_2012_us_cd112.shp",
                 layer = "tl_2012_us_cd112")

# Global Lakes Shapefile, to confine Ohio 9 to land only
lakes <- readOGR("glwd_1.shp", layer = "glwd_1")
erie <- lakes[which(lakes$LAKE_NAME == "Lake Erie"), ]
erieproj <- "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"
erie@proj4string <- CRS(erieproj)
erie <- spTransform(erie, cd115@proj4string)

# State and district identifiers
stfp <- c(48, 17, 42, 26, 39, 37)
stnm <- c("TX", "IL", "PA", "MI", "OH", "NC")
cds <- c("35", "04", "07", "14", "09")

# Figure 1
squiggles <- lapply(1:5, function(i) {
  cd115[cd115$STATEFP == stfp[i] & cd115$CD115FP == cds[i], ]
})
# clip out Lake Erie from Ohio 9
squiggles[[5]] <- rgeos::gDifference(squiggles[[5]], erie)

stcol <- brewer.pal(8, "Set2")[7]
pdf("squiggly_dists.pdf", height = 1.5, width = 7,
    bg = "transparent")
par(mfrow = c(1, 5), mar = c(0, 0, 0, 0))
for(i in 1:5) {
  plot(squiggles[[i]], col = stcol, border = stcol)
  mtext(paste0(stnm[i], "-", cds[i]), side = 1, line = -1)
}
dev.off()

# Figure 2
nc_cd115 <- cd115[cd115$STATEFP == stfp[6], ]  # NC dists only
dem_dists <- c("01", "04", "12")
nc_cd114 <- cd114[cd114$STATEFP == stfp[6], ]
nc_cd112 <- cd112[cd112$STATEFP == stfp[6], ]

dcol <- brewer.pal(9, "Blues")[9]
rcol <- brewer.pal(9, "Reds")[8]
pdf("nc_dists.pdf", height = 4, width = 7)
par(mfrow = c(1, 1), mar = c(0, 0, 0, 0), bg = "transparent")
plot(nc_cd115[!nc_cd115$CD115FP %in% dem_dists, ], col = rcol, border = "grey",
     lwd = 0.5)
plot(nc_cd115[nc_cd115$CD115FP %in% dem_dists, ], col = dcol, border = "grey",
     add = TRUE, lwd = 0.5)
# mtext(paste0(stnm[i], "-", cds[i]), side = 1, line = -1)
dev.off()

# Figure 3 - NC population density versus congressional districts
ctyshp <- "CouSub_2010Census_DP1/CouSub_2010Census_DP1.shp"
cntypop <- readOGR(ctyshp, layer = "CouSub_2010Census_DP1")  # slow (sf maybe
# faster to read in large polygons file)
cntypop@data <- cntypop@data[, c(1:10, length(names(cntypop@data)))]

# http://www-geography.jsu.edu/NC/NCprojection.html
alb <- paste0("+proj=aea +lat_1=34.73 +lat_2=35.672 +lat_0=35.214",
              " +lon_0=-80.391 +x_0=0",
              " +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")

# Project population and district shapes to Albers
nc_cty <- readOGR("tl_2016_37_cousub/tl_2016_37_cousub.shp",
	              layer = "tl_2016_37_cousub")
ncpop <- cntypop[which(substr(cntypop$GEOID10, 1, 2) == "37"), ]
ncpopalb <- spTransform(ncpop, CRS(alb))
nc_ctyalb <- spTransform(nc_cty, CRS(alb))
nc_cdalb <- lapply(list(nc_cd112, nc_cd114, nc_cd115), function(x) {
  spTransform(x, CRS(alb))
})

# Calculate population density
ncpopalb$area <- rgeos::gArea(ncpopalb, byid = TRUE) / 10000 / 100
ncpopalb$popdens <- round(ncpopalb$DP0010001 / ncpopalb$area, 2)

# Rasterize pop density
ncrast <- raster(extent(nc_ctyalb), res = c(1000, 1000))  # 1 km resolution
ncrast[] <- 1:ncell(ncrast)
crs(ncrast) <- crs(ncpopalb)
ncpopr <- rasterize(ncpopalb, ncrast, field = "popdens", fun = "max")
ncpopr[ncpopr > 1000] <- 1000 # set density to 1000 (one pixel seems anomalous)

pdf("nc_pop_sub.pdf", height = 9, width = 7)
par(mfrow = c(3, 1), mar = c(0, 0, 0, 1), bg = "transparent")
cdyr <- paste(c("112th", "114th", "115th"), "Congress")
for(i in 1:3) {
  plot(ncpopr, col = inferno(20), axes = FALSE, box = FALSE,
       interpolate = TRUE)
  plot(nc_cdalb[[i]], border = "cyan", add = TRUE, lwd = 0.1)
  mtext(text = cdyr[i], side = 3, line = -2)
}
dev.off()