Example choropleth using {bcmaps} & data from the B.C. Data Catalogue sourced with {bcdata}

choropleth-example.txt

library(bcdata) #get data from the BC Data Catalogue
library(bcmaps) #get maps
library(tidyr) #clean popn data file
library(dplyr) #summarize data
library(stringr) #clean and align RD names
library(ggplot2) #make plot
library(sf) #clip RD with BC boundaries 

#get and tidy some population change data from the BC Data Catalogue
rd_popn <- bcdc_get_data(
  "bc-population-estimates",
  sheet = "Mun-RD",
  skip = 2) |> 
  janitor::clean_names() |> 
  select(-sgc) |> 
  drop_na("name") |> 
  filter(area_type %in% c("RD", "R")) |> 
  mutate(name = tolower(name),
         name = str_replace(name, " ", "-"))

# get the RD map and match RD names to popn data file (hardest part!)
rdm <- combine_nr_rd() |>
  mutate(
    name = str_remove(ADMIN_AREA_NAME, "Regional District of "),
    name = str_remove(name, " Regional District")
  ) |>
  mutate(name = tolower(name),
         name = str_replace(name, " ", "-")) |>
  mutate(
    name = case_when(
      name == "comox-valley" ~ "comox",
      name == "stikine-region (unincorporated)" ~ "stikine",
      name == "northern-rockies regional municipality" ~ "northern-rockies",
      .default = name
    )
  ) |> 
  st_intersection(bc_bound_hres()) |> #clip RDs with bc boundary file
  rmapshaper::ms_simplify(keep_shapes = TRUE) |> 
  st_make_valid()

# join data and calculate % change
df <- rdm |>
  left_join(rd_popn, keep = TRUE) |> 
  mutate(change = x2019-x2011,
         change2 = change/x2019,
         percent_change = round(100 * change2, digits = 2))

# make a choropleth map
ggplot() +
  geom_sf(data = df, aes(fill = percent_change)) +
  labs(fill = "Population Change\n2011-2019 (%)") +
  theme_minimal() +
  coord_sf(datum = NA)