Homework 4

sf Maps; Animation Plots

Author

Byeong-Hak Choe

Published

May 16, 2025

Modified

May 16, 2025

Direction

  • Please submit your R script for Homework 4 to Brightspace with the name below:

    • danl-310-hw4-LASTNAME-FIRSTNAME.R
      ( e.g., danl-310-hw4-choe-byeonghak.R )

  • The due is May 5, 2025, 2:00 P.M.

  • Please send Byeong-Hak an email (bchoe@geneseo.edu) if you have any questions.




Part 1. ggmap visualization

Question 1

Replicate the following map using Classwork 10 example.

  • Below Google font families are used:

    • Title and Label: Pacifico
    • Subtitle and Legend: Annie Use Your Telescope
    • Caption: Alegreya Sans

  • scale_*_colorblind() is used to make a plot colorblind-friendly.

  • The following option is used to create the figure.

ggsave("PATHNAME_OF_GGPLOT_OUTPUT.png",
       plot = p, height = 12, width = 8, units = "in")
Click to Check the Answer! 
# install.packages(c("osmdata", "osrm"))      # ← run once, then comment-out
library(sf)          # simple-features objects & spatial ops
library(ggmap)       # basemap tiles + ggplot2 integration
library(tidyverse)   
library(ggthemes)    # extra ggplot2 themes (incl. theme_map)
library(osmdata)     # bounding boxes
library(tidygeocoder) # tidy-style geocoding (Nominatim by default)
library(ggrepel)     
library(ggimage)     # drop-in PNG icons
library(showtext)    # Google-font rendering inside R graphics

#–– fonts ––-------------------------------------------------------------------

font_add_google("Roboto", "roboto")
font_add_google("Lato", "lato")
font_add_google("Poppins", "poppins")
font_add_google("Nunito", "nunito")
font_add_google("Annie Use Your Telescope", "annie")
font_add_google("Pacifico", "pacifico")
font_add_google("Alegreya Sans", "aleg")
showtext_auto()                                # enable showtext globally

#–– Stadia Maps key ––----------------------------------------------------------
stadia_api <- "YOUR_STADIA_API_KEY"     # replace w/ your key
register_stadiamaps(stadia_api, write = TRUE)            # store in env file

#–– Geocode landmarks ––--------------------------------------------------------
locations <- 
  tibble(address = c("SUNY Geneseo", "Mount Morris Dam, NY")) |>
  geocode(address) |>               # adds long / lat columns (CRS 4326)
  mutate(address = str_replace_all(address, ", NY", ""))

# osrm & ggmap both attach a function named “mutate()”; avoid a clash by
# not attaching osrm – instead, call it with osrm:: prefix.
# (If you *must* attach osrm, load it *before* ggmap.)

# Pull lon/lat into a plain matrix expected by osrmRoute()
coords_mat <- as.matrix(locations[, c("long","lat")])   # 2×2 matrix

#–– Bike route ––------------------------------------------------------------
route_bike <- osrm::osrmRoute(
  src      = coords_mat[1, ],   # SUNY Geneseo  (lon, lat)
  dst      = coords_mat[2, ],   # Mount Morris Dam, NY (lon, lat)
  overview = "full",            # full-resolution polyline
  osrm.profile = "bike"         # routing profile: foot / bike / car
) |>
  mutate(mode = "Bike")         # tag for legend

#–– Foot route ––------------------------------------------------------------
route_foot <- osrm::osrmRoute(
  src      = coords_mat[1, ],   # SUNY Geneseo  (lon, lat)
  dst      = coords_mat[2, ],   # Mount Morris Dam, NY (lon, lat)
  overview = "full",            # full-resolution polyline
  osrm.profile = "foot"         # routing profile: foot / bike / car
) |>
  mutate(mode = "Foot")         # tag for legend

#–– Driving route ––------------------------------------------------------------
route_car <- osrm::osrmRoute(
  src      = coords_mat[1, ],
  dst      = coords_mat[2, ],
  overview = "full",
  osrm.profile = "car"
) |>
  mutate(mode = "Car")

# combine & prettify
routes <- rbind(route_bike, route_foot, route_car) |>
  mutate(duration = round(duration, 2),                  # minutes
         mode = str_c(mode, ": ", duration, " min."))    # legend label

#–– Build plotting extent ––----------------------------------------------------
bb_routes <- st_bbox(routes)        # xmin, ymin, xmax, ymax (degrees)
pad_x <- (bb_routes$xmax - bb_routes$xmin) * 0.4   # 40 % padding
pad_y <- (bb_routes$ymax - bb_routes$ymin) * 0.4

xlim <- c(bb_routes$xmin - pad_x, bb_routes$xmax + pad_x)
ylim <- c(bb_routes$ymin - pad_y, bb_routes$ymax + pad_y)


#–– Fetch basemap tiles ––------------------------------------------------------
bb_from_limits <- matrix(
  c(xlim["xmin"], xlim["xmax"],   # x-row (min, max)
    ylim["ymin"], ylim["ymax"]),  # y-row (min, max)
  nrow = 2, byrow = TRUE,
  dimnames = list(c("x", "y"), c("min", "max"))
)

# bb_from_limits

man_map <- get_stadiamap(
  bb_from_limits,
  zoom    = 13,                 # ↑ zoom = ↑ detail (and ↑ tiles & time)
  maptype = "stamen_terrain"    # other options: stamen_terrain_labels, stamen_toner, stamen_watercolor, …
)

#–– Plot ––---------------------------------------------------------------------
icon_url <- "https://bcdanl.github.io/lec_figs/marker-icon-red.png"  # pin PNG

p <- ggmap(man_map) +
  geom_sf(
    data = routes,
    aes(color = mode),
    inherit.aes = FALSE,   # don’t map long/lat twice
    linewidth = rel(2),
    alpha = 0.75
  ) +
  geom_image(
    data = locations,
    aes(x = long, y = lat, image = icon_url),
    size = .025
  ) +
  geom_label_repel(
    data = locations,
    aes(long, lat, label = address),
    size = rel(12),
    box.padding = 1.75,
    family = "pacifico"
  ) +
  coord_sf(                      # keep CRS 4326; clamp to padded bbox
    crs   = st_crs(4326),
    xlim  = xlim,
    ylim  = ylim,
    expand = FALSE
  ) +
  scale_color_colorblind() +
  guides(
    color = guide_legend(
      title.position = "top",
      label.position = "bottom",
      keywidth = rel(10.5)
    )
  ) +
  labs(
    color    = "",
    title    = "Bike vs. Car vs. Foot",
    subtitle = "How would you like to go from SUNY Geneseo to Mount Morris Dam?",
    caption  = "Source: OpenStreetMap (OSM) and Open Source Routing Machine (OSRM)"
  ) +
  theme_map() +
  theme(
    legend.position  = "top",
    plot.title       = element_text(face = "bold", 
                                    size = rel(9), 
                                    family = 'pacifico', 
                                    hjust = 0),
    plot.subtitle    = element_text(face = "bold.italic", 
                                    size = rel(7), 
                                    family = 'annie', 
                                    lineheight = .25,
                                    hjust = 0),
    plot.caption    = element_text(face = "italic", 
                                   size = rel(4), 
                                   family = 'aleg'),
    legend.text      = element_text(face = "bold", 
                                    family = 'annie', 
                                    size = rel(5.5)),
    legend.background = element_rect(fill = NA, color = NA),
    legend.box.margin = margin(-15,0,-7.5,0),
  )

# set a plotting device size manually for forcing the right canvas size
dev.new(width = 8, height = 12)
p

ggsave("/Users/bchoe/Documents/websites/bcdanl.github.io/lec_figs/map_geneseo_mtmorrisdam3.png",
       plot = p, height = 12, width = 8, units = "in")



Part 2. Simple Feature (sf) and Aninimation Plots

Consider the county-level climate opinion data from Yale Program on Climate Change Communication:

climate_opinion <- read_csv(
  'https://bcdanl.github.io/data/climate-change-public-opinion.csv')


Question 2

  • Replicate the following map about variable human:
    • human: Estimated percentage who think that global warming is caused mostly by human activities


  • For R packages and map data, use the following:
# Load required packages ---------------------------------------------------
library(tidyverse)   
library(sf)          # simple-features objects for spatial data
library(tigris)      # fast access to Census TIGER/Line shapefiles
library(ggthemes)    
library(gganimate)   # turn static ggplots into animations
library(showtext)    # Google-font rendering inside R graphics

# Register Google fonts so they render the same on any machine -------------
font_add_google("Roboto",       "roboto")  # body text
font_add_google("Alegreya Sans","aleg")    # title text
showtext_auto()                            # activate showtext globally

# 1) Download   US state & county geometries (simplified “cartographic” 
#    boundaries) at 1:5 million scale, then shift Alaska, Hawai‘i & PR closer ---

us_states <- states(
  cb         = TRUE,   # use generalized cartographic boundary files
  resolution = "5m",   # 1 : 5 000 000
  year       = 2024    # 2024 vintage
) |> 
  shift_geometry()     # moves AK/HI & PR into a CONUS-friendly layout

us_counties <- counties(
  state      = NULL,   # all states + DC + territories
  cb         = TRUE,
  resolution = "5m",
  year       = 2024
) |> 
  shift_geometry()

# 2) Build a bounding box that covers only the 50 states -------------------
#    (exclude DC + territories so the map crops nicely)

main_states <- states(cb = TRUE, resolution = "5m", year = 2024)

keep <- state.abb                        # built-in vec of 50 state abbr.
main_states <- main_states |>
  filter(STUSPS %in% keep) |>
  shift_geometry()

bb <- sf::st_bbox(main_states)           # xmin / ymin / xmax / ymax vector

# 3) Quick static map -------------------------------------------------------
ggplot() +
  geom_sf(data = us_counties,
          fill  = NA,        # counties as thin outlines
          color = "grey80",
          linewidth = 0.01) +
  geom_sf(data = us_states,
          fill  = NA,        # thicker state borders
          color = "black") +
  coord_sf(                   # crop to 50-state bounding box
    xlim   = c(bb["xmin"], bb["xmax"]),
    ylim   = c(bb["ymin"], bb["ymax"]),
    expand = FALSE
  ) +
  theme_map()                 # clean, border-free background
Click to Check the Answer! 
library(tidyverse)
library(sf)
library(tigris)
library(ggthemes)
library(gganimate)
library(showtext)
font_add_google("Roboto", "roboto")
font_add_google("Lato", "lato")
font_add_google("Poppins", "poppins")
font_add_google("Nunito", "nunito")
font_add_google("Annie Use Your Telescope", "annie")
font_add_google("Pacifico", "pacifico")
font_add_google("Alegreya Sans", "aleg")
showtext_auto()  # Enables Google fonts in R plots.

# 1) pull every U.S. state/county as a simplified (“cartographic boundary”) file:
us_states <- states(
  cb         = TRUE,        # cartographic boundary (generalized)
  resolution = "5m",        # 1:5 million scale (also "500k" or "20m")
  year       = 2024         # vintage
  ) |> 
  shift_geometry()

us_counties <- counties(
  state = NULL, 
  cb = TRUE, 
  resolution = "5m", 
  year = 2024) |> 
  shift_geometry()


# 1. get the box:
# pull everything that states() gives us
main_states <- states(cb = TRUE, resolution = "5m", year = 2024)

# keep only the 50 US states
keep <- c(state.abb)
main_states <- main_states |> 
  filter(STUSPS %in% keep) |> 
  shift_geometry()
  

bb <- sf::st_bbox(main_states)

# climate data for human belief
climate_opinion <- read_csv(
  'https://bcdanl.github.io/data/climate-change-public-opinion.csv')

climate_opinion_human <- climate_opinion |> 
  filter(variable %in% c("human")) |> 
  select(GeoID, GeoName,
         starts_with("x")) |> 
  pivot_longer(cols = x2014:x2024,
               names_to = "year",
               values_to = "belief") |> 
  mutate(year = str_replace_all(year, "x", ""),
         year = as.integer(year)) |> 
  rename(GEOID = GeoID)

unique(climate_opinion_human$year)


# us_counties6 <- rbind(us_counties |> mutate(year = 2014),
#                       us_counties |> mutate(year = 2016),
#                       us_counties |> mutate(year = 2018),
#                       us_counties |> mutate(year = 2020),
#                       us_counties |> mutate(year = 2022),
#                       us_counties |> mutate(year = 2024)) |> 
#   left_join(climate_opinion_human) |> 
#   filter(STUSPS != "PR")

us_counties_join <- us_counties |> 
  left_join(climate_opinion_human) |> 
  filter(STUSPS != "PR") |> 
  filter(year %in% seq(2014,2024, 2))

# facet
ggplot() +
  geom_sf(data = us_counties_join,
          aes(geometry = geometry, fill = belief/100),
          color = NA, size = 0.01) +
  geom_sf(data = us_states |> filter(STUSPS != "PR"),
          fill = NA, color = "grey", size = .25) +
  coord_sf(
    xlim = c(bb["xmin"], bb["xmax"]),
    ylim = c(bb["ymin"], bb["ymax"]),
    expand = FALSE
  ) +
  scale_fill_viridis_c(
    option = "inferno",
    limits = c(min((us_counties6$belief)/100), max((us_counties6$belief)/100)), 
    breaks = seq(.4, 1, .1),
    labels = scales::percent_format(accuracy = 1),
    guide = guide_colorbar(
      title           = "% Yes",
      title.position  = "left",
      barwidth        = unit(17.5, "cm"),
      # barheight       = unit(0.5, "cm"),
      direction       = "horizontal",
      nbin            = 500,
      raster          = FALSE,
      ticks           = TRUE,
      draw.llim       = FALSE,
      reverse         = FALSE
    )
  ) +
  facet_wrap(~year) +
  labs(
    title = "Is global warming caused mostly by human activities?",
    caption = "Source: Yale Climate Change Communication, 2025"
  ) +
  theme_map() +
  theme(plot.title = element_text(hjust = .5,
                                  margin = margin(0,0,0,0),
                                  family = 'aleg',
                                  face = 'bold',
                                  size = rel(2.75)),
        plot.subtitle = element_text(hjust = .75,
                                     vjust = -87.5,
                                     family = 'roboto',
                                     face = 'bold',
                                     size = rel(1.75)),
        plot.caption = element_text(family = 'roboto',
                                    face = 'italic'),
        strip.background = element_rect(fill = NA,
                                         color = NA),
        strip.text = element_text(
          family = 'roboto',
          face = 'bold',
          size = rel(1.75)),
        legend.position = 'top',
        legend.justification = 'center',
        legend.title = element_text(family = 'roboto', 
                                    face = 'bold',
                                    vjust = .75)
  )



Question 3

  • Animate the maps for the even-numbered years used in Question 2:

  • To view this animation bigger, click this.

  • For animate(), use the following parameter setting:

ggplot() +
  ...
 labs(subtitle = "Year: {floor(frame_time)}") +
  ...

anim <- animate(p_anim, 
                nframes = 20, 
                fps = 5, 
                width = 8, height = 5,
                units = "in",
                res = 900)  # resolution can be reduced

# To save animation as PATHNAME_FOR_YOUR_FILE.gif
anim_save("PATHNAME_FOR_YOUR_FILE.gif", anim)

  • You can open “PATHNAME_FOR_YOUR_FILE.gif” in any web browser (Chrome is recommended).

  • transition_time(year) treats year as a continuous numeric variable.

  • After you call animate(), gganimate:

  1. Spans the full numeric range of the variable — here 2014 → 2024, a width of 10 years.
  2. Slices that span into nframes equally spaced instants.
    • Here nframes = 20, so the frame-to-frame step is \[ \begin{align} \Delta = \frac{2024 - 2014}{20 - 1} \;=\; \frac{10}{19} \;\approx\; 0.526 \text{ years}. \end{align} \]
  3. Assigns each frame its own “virtual” time stamp:
    • \(t_0 = 2014\), \(t_1 = 2014.526\), \(t_2 = 2015.053\), …
Click to Check the Answer! 
p <- ggplot() +
  geom_sf(data = us_counties_join,
          aes(geometry = geometry, fill = belief/100),
          color = NA, size = 0.01) +
  geom_sf(data = us_states |> filter(STUSPS != "PR"),
          fill = NA, color = "grey", size = .25) +
  coord_sf(
    xlim = c(bb["xmin"], bb["xmax"]),
    ylim = c(bb["ymin"], bb["ymax"]),
    expand = FALSE
  ) +
  scale_fill_viridis_c(
    option = "inferno",
    limits = c(min((us_counties6$belief)/100), max((us_counties6$belief)/100)), 
    breaks = seq(.4, 1, .1),
    labels = scales::percent_format(accuracy = 1),
    guide = guide_colorbar(
      title           = "% Yes",
      title.position  = "left",
      barwidth        = unit(17.5, "cm"),
      # barheight       = unit(0.5, "cm"),
      direction       = "horizontal",
      nbin            = 500,
      raster          = FALSE,
      ticks           = TRUE,
      draw.llim       = FALSE,
      reverse         = FALSE
    )
  ) +
  labs(
    title = "Is global warming caused mostly by human activities?",
    subtitle = "Year: {floor(frame_time)}",
    caption = "Source: Yale Climate Change Communication, 2025"
  ) +
  theme_map() +
  theme(plot.title = element_text(hjust = .5,
                                  margin = margin(0,0,-10,0),
                                  family = 'aleg',
                                  face = 'bold',
                                  size = rel(2.75)),
        plot.subtitle = element_text(hjust = .75,
                                     vjust = -87.5,
                                     family = 'roboto',
                                     face = 'bold',
                                     size = rel(1.75)),
        plot.caption = element_text(family = 'roboto',
                                    face = 'italic'),
        legend.position = 'top',
        legend.justification = 'center',
        legend.title = element_text(family = 'roboto', 
                                    face = 'bold',
                                    vjust = .75)
  )
# p

p_anim <- p +
  transition_time(year) +        # animate through year
  ease_aes('linear')

anim <- animate(p_anim, 
                # nframes = length(years)*18, 
                nframes = 20, 
                fps = 5, 
                width = 8, height = 5,
                units = "in",
                res = 900)
anim_save("/Users/bchoe/Downloads/sf_cc_human2025.gif", anim)
Back to top