Continental Counties vs Health Service Areas (HSA) [USA]

.chart_usa_counties_vs_hsa.r

library(tigris)
library(ggplot2)
library(dplyr)
library(readxl)
library(sf)

sf <- 2
width <- 1920 * sf
height <- 1080 * sf

options(vsc.dev.args = list(width = width, height = height, res = 72 * sf))

# Load the US counties shapefile for 2020
us_counties <- counties(cb = TRUE, year = 2020, class = "sf")

# Define FIPS codes for the continental US (48 contiguous states)
continental_fips <- c(
  "01", "04", "05", "06", "08", "09", "10", "11", "12", "13", "16", "17", "18",
  "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31",
  "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42", "44", "45",
  "46", "47", "48", "49", "50", "51", "53", "54", "55", "56"
)
# continental_fips <- c("06")

# Filter to only include counties in the continental US
continental_counties <- us_counties |> filter(STATEFP %in% continental_fips)

# USA Continental Counties 2020
chart <- ggplot(data = continental_counties) +
  geom_sf() +
  theme_minimal() +
  labs(
    title = "Continental US Counties Map (2020)",
    caption = "Data source: tigris package"
  ) +
  theme(
    axis.text = element_blank(), # Remove axis labels
    axis.ticks = element_blank()
  ) # Remove axis ticks

ggplot2::ggsave(
  filename = "chart1.png", plot = chart, width = width, height = height,
  units = "px", dpi = 72 * sf, device = grDevices::png, type = c("cairo")
)

# USA Continental HSA's 2020
# Download the Excel file with HSA data
download.file(
  paste0(
    "https://seer.cancer.gov/seerstat/variables/countyattribs/",
    "Health.Service.Areas.xls"
  ),
  "/tmp/data.xlsx"
)
hsa_data <- read_xls("/tmp/data.xlsx", sheet = "HSA (NCI Modified)") |>
  select(1, 4) |>
  setNames(c("hsa", "fips")) |>
  mutate(
    fips = sprintf("%05s", fips), # Zero-pad FIPS codes to 5 digits
    hsa = as.character(hsa) # Ensure HSA is a character vector
  )

continental_hsa <- continental_counties |> left_join(hsa_data, by = c("GEOID" = "fips"))

continental_hsa_merged <- continental_hsa |>
  group_by(hsa) |>
  summarize(geometry = st_union(geometry)) |>
  ungroup()

chart <- ggplot() +
  geom_sf(data = continental_counties, fill = "grey80", color = "white") + # Base county map
  geom_sf(data = continental_hsa_merged, fill = NA, color = "blue", size = 0.5) + # HSA boundaries
  theme_minimal() +
  labs(
    title = "Continental US Counties with Health Service Areas (2020)",
    caption = "Data source: tigris package"
  ) +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank()
  )

ggplot2::ggsave(
  filename = "chart2.png", plot = chart, width = width, height = height,
  units = "px", dpi = 72 * sf, device = grDevices::png, type = c("cairo")
)

continental_hsa_merged <- continental_hsa |>
  filter(STATEFP == "06") |>
  group_by(hsa) |>
  summarize(geometry = st_union(geometry)) |>
  ungroup()

# Combined map: Colored counties with blue HSA boundaries
chart <- ggplot() +
  # Color the counties based on the HSA region
  geom_sf(data = continental_counties |> filter(STATEFP == "06"), fill = "grey80", color = "white") + # Base county map
  geom_sf(data = continental_hsa_merged, fill = NA, color = "blue", size = 0.5) + # HSA boundaries
  # scale_fill_viridis_d(option = "plasma") +
  theme_minimal() +
  labs(
    title = "California Counties with Health Service Areas (2020)",
    fill = "HSA Region",
    caption = "Data source: tigris package; seer.cancer.gov"
  ) +
  theme(
    axis.text = element_blank(), # Remove axis labels
    axis.ticks = element_blank() # Remove axis ticks
  )


ggplot2::ggsave(
  filename = "chart3.png", plot = chart, width = width, height = height,
  units = "px", dpi = 72 * sf, device = grDevices::png, type = c("cairo")
)

https://www.mortality.watch/charts/list.html#continental-counties-vs-health-service-areas-hsa-usa