Add VT2-examples: all 6 chapters of 'Visualize This' (110 files) — data scraping, ETL, trends, categorical, multivariate, spatial, distribution analysis

VT2-examples/vt2-ch7-codedata-dir/map-distributions.R

@@ -0,0 +1,145 @@
+library(sf)
+library(cartogram)
+
+
+#
+# Load data
+#
+
+# Boundaries
+statefp <- "data/cb_states/cb_2022_us_state_20m.shp"
+statebnds <- st_read(statefp)
+
+# Conterminous subset
+inconterm <- !(statebnds$STUSPS %in% c("AK", "HI", "PR"))
+conterm <- statebnds[inconterm,]
+conterm_p4s <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lon_0=97.2w"
+conterm_albers <- st_transform(conterm, conterm_p4s)
+
+# Load population data
+fileloc <- "data/DECENNIALPL2020.P1/DECENNIALPL2020.P1.csv"
+pop <- read.csv(fileloc, stringsAsFactors = FALSE)
+
+
+#
+# Join datasets
+#
+
+# Look at data structure
+pop[1:10, 1:5]
+conterm_albers[1:10, 1:6]
+
+# Match
+conterm_albers$asianpop <-
+    pop$P1_006N[match(conterm_albers$AFFGEOID, pop$GEO_ID)]
+conterm_albers$totalpop <-
+    pop$P1_001N[match(conterm_albers$AFFGEOID, pop$GEO_ID)]
+
+# Calculate rate
+conterm_albers$asianrate <- 
+    conterm_albers$asianpop / conterm_albers$totalpop
+
+# Check data
+conterm_albers[1:10, 
+    c("NAME", "totalpop", "asianpop", "asianrate")]
+
+
+
+#
+# Choropleth
+#
+
+# Summary statistics
+summary(conterm_albers$asianrate)
+
+# Define color
+pickCol <- function(x) {
+    propbreaks <- c(0, .02, .04, .06, .08, 1)
+    shades <- c('#d9e0e1', '#acc0c4', '#7fa1a8', 
+                '#51828c', '#196572')
+    i <- max(which(propbreaks <= x))
+    return(shades[i])
+}
+conterm_albers$col <- 
+    sapply(conterm_albers$asianrate, pickCol)
+
+# Draw map
+par(mar=c(1,1,1,1))
+plot(st_geometry(conterm_albers), 
+     col = conterm_albers$col,
+     border = "white")
+legend("bottomleft", 
+       legend = c("0-.01", ".02-.03", ".04-.05", 
+                  ".06-.07", ".08+"), 
+       fill = c('#d9e0e1', '#acc0c4', '#7fa1a8', 
+                '#51828c', '#196572'),
+       title = "Asian Population")
+
+
+
+#
+# Cartogram
+#
+
+# Dorling calculation
+conterm_dorling <- cartogram_dorling(conterm_albers, 
+                                     "asianpop")
+plot(st_geometry(conterm_albers),
+     border = "#cccccc", 
+     lwd = .6)
+plot(st_geometry(conterm_dorling), 
+     col = conterm_albers$col,
+     add = TRUE)
+
+# Non-contiguous
+conterm_ncont <- cartogram_ncont(conterm_albers, 
+                                 "asianpop")
+plot(st_geometry(conterm_albers),
+     border = "#cccccc", 
+     lwd = .6)
+plot(st_geometry(conterm_ncont), 
+     col = conterm_albers$col,
+     add = TRUE)
+
+# Contiguous cartogram
+conterm_cont <- cartogram_cont(conterm_albers,
+                               "asianpop")
+plot(st_geometry(conterm_cont), 
+     col = conterm_albers$col)
+
+
+#
+# Dot density
+#
+
+# Load block boundaries
+bgfp <- "data/cb_2022_32_bg_500k/cb_2022_32_bg_500k.shp"
+nvblocks <- st_read(bgfp)
+
+# Block level population
+popfp <- "data/DEC2020.P8_bg/DECENNIALDHC2020.P8-Data.csv"
+nvpop <- read.csv(popfp)
+
+# Clark county blocks
+clarkblocks <- nvblocks[nvblocks$COUNTYFP == "003",]
+
+# County level ids in nvpop
+clark_geoid <- "1500000US32003"
+nvpop$countyid <- substr(nvpop$GEO_ID, 1, 
+                         nchar(clark_geoid))
+clarkpop <- nvpop[nvpop$countyid == clark_geoid,]
+
+# Merge population by county
+imatch <- match(clarkblocks$AFFGEOID, clarkpop$GEO_ID)
+clarkblocks$asianpop <- clarkpop$P8_006N[imatch]
+
+# Dots
+clarkdots <- st_sample(clarkblocks, 
+    as.integer(clarkblocks$asianpop/10))
+
+# Map and dot
+plot(st_geometry(clarkblocks),
+     border = "#cccccc", lwd = .5)
+plot(clarkdots, pch=19, cex=0.1, 
+     col="#19657250", 
+     add = TRUE)