ggplot2-map-example.r

# load packages ----
library(sf)
library(dplyr)
library(httr)
library(ggplot2)
library(ggthemes)
library(ggview)
library(mapview)

# load spatial data ----
# load ky county polygons from KY Open Data arcgis api
# from here: https://opengisdata.ky.gov/datasets/kygeonet::kentucky-county-polygons/about
url <- parse_url("https://services3.arcgis.com/ghsX9CKghMvyYjBU/arcgis/rest/services/Ky_County_Polygons_WM/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson")
request <- build_url(url)

ky_county_poly <- st_read(request) |>
  st_transform(6473) |>
  janitor::clean_names()

ky_county_poly_select <- ky_county_poly |>
  select(fips_id, name2, pop10)

# load county level data ----
# from here: https://www.countyhealthrankings.org/health-data/kentucky/data-and-resources
# Define the URL of the Excel file
url <- "https://www.countyhealthrankings.org/sites/default/files/media/document/2025%20County%20Health%20Rankings%20Kentucky%20Data%20-%20v3.xlsx"  # Replace with your Office View URL

# Specify the destination file path
destfile <- "data/ky_county_health_rankings_2025.xlsx"

# Download the file
GET(url, write_disk(destfile, overwrite = TRUE))

# Read the Excel file
ky_county_health_rankings_2025 <- readxl::read_excel(destfile, sheet = "Additional Measure Data", skip = 1) |>
  janitor::clean_names()

ky_county_health_rankings_2025_select <- ky_county_health_rankings_2025 |>
  select(fips, life_expectancy) |>
  mutate(
    fips = as.integer(fips),
    life_expectancy = round(life_expectancy, 2),
    # create intervals using quantile breaks
    life_expectancy_bin = cut(
      life_expectancy,
      breaks = quantile(
        life_expectancy,
        probs = seq(0, 1, length.out = 6),
        na.rm = TRUE
      ),
      include.lowest = TRUE
    ),
    # create intervals based on pretty breaks
    life_expectancy_pretty = cut(
      life_expectancy,
      breaks = classIntervals(life_expectancy, n = 5, style = "pretty")$brks,
      digits = 2
    )
  )

# combine county data to spatial data ----
ky_county_poly_join <- left_join(
  ky_county_poly_select,
  ky_county_health_rankings_2025_select,
  by = join_by(fips_id == fips)
)

# create map ----
# create data for labels
myLabels <- ky_county_poly_join |>
  filter(
    life_expectancy %in% head(sort(ky_county_poly_join$life_expectancy), 20) |
      pop10 %in% tail(sort(ky_county_poly_join$pop10), 10)
  )

myPlot <- ggplot(ky_county_poly_join) +
  geom_sf(aes(fill = life_expectancy), color = "grey50", linewidth = 0.75) +
  ggrepel::geom_label_repel(
    data = myLabels,
    aes(label = name2, geometry = geometry),
    size = 2,
    segment.color = "black",
    stat = "sf_coordinates",
    min.segment.length = 0
  ) +
  scale_fill_viridis_c() +
  theme_void() +
  theme(
    plot.margin = margin(0, 1, 0, 1, "cm"),
    legend.position = "bottom"
  ) +
  labs(
    title = "Life Expectancy by County (Kentucky)",
    fill = NULL
  ) +
  canvas(
    width = 8,
    height = 4,
    units = "in"
  )

myPlot

save_ggplot(myPlot, file = "export/myPlot.png")

myPlotBinned <- ggplot(ky_county_poly_join) +
  geom_sf(aes(fill = life_expectancy_pretty), color = "grey50", linewidth = 0.75) +
  ggrepel::geom_label_repel(
    data = myLabels,
    aes(label = name2, geometry = geometry),
    size = 2,
    segment.color = "black",
    stat = "sf_coordinates",
    min.segment.length = 0
  ) +
  scale_fill_brewer(
    palette = "GnBu",
    direction = -1,
    guide = guide_legend(
      keyheight = unit(3, units = "mm"),
      keywidth = unit(12, units = "mm"),
      label.position = "bottom",
      nrow = 1
    )
  ) +
  theme_void() +
  theme(
    plot.margin = margin(0, 1, 0, 1, "cm"),
    legend.position = "bottom"
  ) +
  labs(
    title = "Life Expectancy by County (Kentucky)",
    fill = NULL
  ) +
  canvas(
    width = 8,
    height = 4,
    units = "in"
  )

myPlotBinned

save_ggplot(myPlotBinned, file = "export/myPlotBinned.png")