Generate animated transit maps using GTFS feeds. See https://sno.ws/transit-maps

gtfs_map_generator.R

##### Setup #####

# Load the necessary libraries
library(tidytransit) # Parse GTFS feeds
library(dplyr)
library(tidyr)
library(lubridate)
library(sf)
library(ggplot2)
library(gganimate)   # Animate the ggplot + tween between frames
library(ggspatial)   # Add map tiles
library(stinepack)   # Interpolate between stop times

# Set some simple parameters for the script
feed <- "CTA GTFS"  # Feed list can be found at tidytransit::feedlist
transit_type = 1  # See tidytransit::route_type_names

# Geographic area to plot. See boundingbox.klokantech.com for tool
bbox <- c(-87.736558, 41.807137, -87.556657, 41.943711)

# Trip departure date and times. Note that these are local time
dep_date = lubridate::today() - days(1) 
min_dep_time = "08:00:00"
max_arv_time = "09:00:00"

# Animation length scalar. Smaller = shorter animation
# A value of 1 means 1 minute of animation per hour depending on fps
min_to_hour_ratio = 0.25

# Animation FPS. Higher values = smoother anim. but bigger file
frames_per_second = 60



##### Get GTFS feed #####

# Download a GTFS feed and parse it with tidytransit
gtfs <- tidytransit::feedlist %>% 
  filter(t == feed) %>%
  pull(url_d) %>%
  read_gtfs()

# Keep only the routes (and corresponding trips and stops) specified
route_ids <- gtfs$routes %>%
  filter(route_type == transit_type) %>%
  pull(route_id)

trip_ids <- gtfs$trips %>%
  filter(route_id %in% route_ids) %>%
  pull(trip_id)

stop_ids <- gtfs$stop_times %>%
  filter(trip_id %in% trip_ids) %>%
  pull(stop_id) %>%
  unique()

# Load the static shapes of the routes and stops
gtfs_sf <- gtfs %>%
  tidytransit::gtfs_as_sf()

route_shapes <- gtfs_sf %>%
  get_route_geometry(route_ids) %>%
  left_join(gtfs$routes, by = "route_id") %>%
  mutate(route_color = paste0("#", route_color))

stop_shapes <- gtfs_sf$stops %>%
  filter(stop_id %in% stop_ids)



##### Get interpolated trips #####

# Each route is actually a series of points (held in the shapes
# section of the feed) that are joined into a line. We can use
# these points to define the waypoints that vehicles need to follow
route_waypoints <- gtfs$trips %>%
  filter(route_id %in% route_ids) %>%
  distinct(route_id, shape_id) %>%
  left_join(gtfs$shapes, by = "shape_id") %>%
  select(-shape_pt_sequence) %>%
  rename(
    lat = shape_pt_lat,
    lon = shape_pt_lon,
    dist = shape_dist_traveled
  )

# Create a data frame of stops with time stopped and geographic
# location. Only include trips within the time boundaries
# specified at top of script
stops_df <- gtfs$trips %>%
  filter(route_id %in% route_ids) %>%
  inner_join(
    filter_stop_times(gtfs, dep_date, min_dep_time, max_arv_time),
    by = "trip_id"
  ) %>%
  inner_join(gtfs$stops, by = "stop_id") %>%
  distinct(
    route_id, shape_id, trip_id, arrival_time,
    lat = stop_lat, lon = stop_lon, dist = shape_dist_traveled
  )

# For each trip, get ALL the waypoints along the route
# (basically the points between each stop for each trip)
waypoints_df <- stops_df %>%
  distinct(trip_id, shape_id) %>%
  inner_join(route_waypoints, by = "shape_id")

# Combine the known stop times with the unknown waypoint times,
# then fill in the missing waypoint data with time-series
# imputation based on the distance column
final_df <- stops_df %>%
  bind_rows(waypoints_df) %>%
  group_by(trip_id) %>%
  filter(sum(!is.na(arrival_time)) > 1) %>%
  mutate(
    time = as.POSIXct(
      as.character(arrival_time), format = "%H:%M:%S",
      tz = gtfs$agency$agency_timezone
    ) - days(1)
  ) %>%
  arrange(trip_id, dist) %>%
  group_by(trip_id, dist) %>%
  filter(row_number() == 1) %>%
  group_by(trip_id) %>%
  mutate(
    time = as.POSIXct(
      stinepack::na.stinterp(
        as.numeric(time),
        along = dist,
        na.rm = FALSE
      ),
      origin = "1970-01-01",
      tz = gtfs$agency$agency_timezone
    )
  ) %>%
  ungroup() %>%
  filter(!is.na(time))



##### Create animation #####

# If bounding box is not set, will use bbox of routes
if (!exists("bbox")) bbox <- st_bbox(route_shapes)

# Create a plot object to animate. Removing annotation_map_tile() will
# signficantly speed up drawing each frame
p <- ggplot(route_shapes) +
  annotation_map_tile(type = "cartolight", zoomin = 0, progress = "none") +
  geom_sf(data = route_shapes, aes(color = route_color), size = 1.5) +
  geom_sf(
    data = stop_shapes, stroke = 1.5,
    size = 4, shape = 21,
    color = "#000000", fill = "#ffffff"
  ) +
  coord_sf(xlim = c(bbox[1], bbox[3]), ylim = c(bbox[2], bbox[4])) +
  geom_point(
    data = final_df,
    aes(x = lon, y = lat, group = trip_id),
    size = 2.5,
    shape = 15
  ) +
  scale_color_identity() +
  transition_components(time) +
  ease_aes("sine-in-out") +
  theme_void() +
  labs(
    title = gtfs$agency$agency_name,
    subtitle = "{frame_time}") + 
  theme(
    legend.position = "none",
    plot.title = element_text(face = "bold", size = 24),
    plot.subtitle = element_text(size = 18)
  )

# Create the actual frames of animation
frames <- as.numeric(hms(max_arv_time) - hms(min_dep_time)) * min_to_hour_ratio
plot_mg <- animate(
  plot = p, nframes = frames, fps = frames_per_second,
  width = 800, height = 800,
  device = "png",
  renderer = file_renderer(dir = ".", overwrite = TRUE)
)

Hello Dan, great work on this gist, it has been very useful.
I work at a local government transport agency in Auckland, New Zealand. We would like to adapt this code to produce some visualizations for the public. I can't see a license for this code, may we use it to create derivative code and credit you as Dan Snow (http://sno.ws/)?

@TLDuhamel Glad to see people are getting some use from it! Consider it MIT licensed. Would love to see your results if you want to drop them in this thread once you're done!

Also note, I've been playing with a better version of this code for a bit with more features: actual stopping at stops, proper speed up + slow down, fixes for overlapping tracks, etc. I'll try to find some time in the next few weeks to write up a post about it.

@dfsnow I can now share the published animation, which was launched on Facebook: https://www.facebook.com/akltransport/videos/295156760023471/