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)