Code for my map for day 12 of the 2021 30 Day Map Challenge (November 12 2021: "Population")

20211112.R

library(sf)
library(tidycensus)
library(ggplot2)
library(dplyr)

acs_pop <- get_acs(
  geography = "block group",
  variables = "B01003_001",
  state = "MA",
  year = 2019,
  geometry = TRUE
)

acs_pop <- acs_pop |> 
  mutate(area = units::drop_units(units::set_units(st_area(geometry), "miles^2")),
         area_adjust = estimate / area,
         area_adjust = area_adjust + 1)


ggplot(acs_pop) + 
  geom_sf(aes(fill = area_adjust), color = "black", size = 0.035) + 
  scale_fill_distiller(palette = "Purples", direction = 1, trans = "log", labels = scales::number_format(big.mark = ","), breaks = c(1, 10, 100, 1000, 10000), name = expression("Population density (# / mile"^2 ~ ")")) + 
  theme_void() + 
  theme(legend.position = "bottom",
        legend.key.width = unit(0.6, "in"),
        plot.background = element_rect(fill = "white", color = NA))

ggsave("20211112.png", width = 16, height = 10)