Animating public transport networks on 3D stacked maps

img { border: 5px solid #FFFFFF; }

install.packages('easypackages')
easypackages::packages('geobr', 'aopdata', 'gtfs2gps', 
                       'data.table', 'sf', 'viridis',
                       'magrittr', 'dplyr', 'ggnewscale',
                       'ggplot2', 'gganimate', 'av')



###### Functions to tilt sf ---------------------------------------------
# more info at https://www.urbandemographics.org/post/figures-map-layers-r/
rotate_data <- function(data, x_add = 0, y_add = 0) {
  
  shear_matrix <- function(){ matrix(c(2, 1.2, 0, 1), 2, 2) }
  
  rotate_matrix <- function(x){ 
    matrix(c(cos(x), sin(x), -sin(x), cos(x)), 2, 2) 
  }
  data %>% 
    dplyr::mutate(
      geometry = .$geometry * shear_matrix() * rotate_matrix(pi/20) + c(x_add, y_add)
    )
}

rotate_data_geom <- function(data, x_add = 0, y_add = 0) {
  shear_matrix <- function(){ matrix(c(2, 1.2, 0, 1), 2, 2) }
  
  rotate_matrix <- function(x) { 
    matrix(c(cos(x), sin(x), -sin(x), cos(x)), 2, 2) 
  }
  data %>% 
    dplyr::mutate(
      geom = .$geom * shear_matrix() * rotate_matrix(pi/20) + c(x_add, y_add)
    )
}


###### DATA: city boundary ---------------------------------------------
code_muni <- geobr::lookup_muni(name_muni = 'Fortaleza')$code_muni
city <- read_municipality(code_muni = code_muni)




###### DATA: population from aopdata ---------------------------------------------

#' more info at https://www.ipea.gov.br/acessooportunidades/en/
aop <- aopdata::read_access(city = 'fortaleza', mode="public_transport", geometry = T)

# select areas with population
aop <- subset(aop, P001 > 0)




###### DATA: pubic transport data ---------------------------------------------

# read GTFS data

# I'm using a local file, but you can use a sample GTFS by uncommenting the line below 
gtfs_file <-  './GTFS_fortaleza_20211020.zip'
#gtfs_file <- system.file("extdata/fortaleza.zip", package = "gtfs2gps")
gtfs_raw <- gtfs2gps::read_gtfs(gtfs_file)

# select trips btween 7am and 8am
gtfs <- gtfs2gps::filter_day_period(gtfs_raw, period_start="07:00:00", period_end = '08:00:00')

# select a sample of n routes
set.seed(42)
nroutes <- 80
myroutes <- sample(x = unique(gtfs$routes$route_id), size = nroutes, replace = F)
gtfs <- gtfs2gps::filter_by_route_id(gtfs, route_ids = myroutes)

# get route geometries
routes <- gtfs2gps::gtfs_shapes_as_sf(gtfs)

# convert GTFS to GPS-like data.table
gps <- gtfs2gps(gtfs, parallel = T, spatial_resolution = 200)
gps <- gps[ between(departure_time, as.ITime('07:00:00'), as.ITime('08:00:00'))]

# Convert "GPS" points of vehicle positions into sf
gps_points <- sfheaders::sf_point(gps, x = "shape_pt_lon" , y = "shape_pt_lat", keep = T)
sf::st_crs(gps_points) <- 4326

# rotate vehicle positions and back to points
gps_points_r <- rotate_data(gps_points, y_add = 0.3)
gps_points_df <- sfheaders::sf_to_df(gps_points_r, fill = T)

# remove missing
gps_points_df <- subset(gps_points_df, !is.na(departure_time))


###### Plot ---------------------------------------------

# annotate parameters
x = -80.92
clr="gray40"
sz = 4

# city boundary
temp1  <- ggplot() +
              geom_sf(data = city %>% rotate_data_geom(y_add = .03), color="gray90", fill="gray90", show.legend = FALSE) +
              annotate("text", label="City boundary", x=x, y=9.10, hjust = 0, color=clr, size=sz) +
              labs(caption = "Image by @UrbanDemog")

temp2 <- temp1 +
  
  # Population Income
  geom_sf(data = subset(aop, P001>0) %>% rotate_data(y_add = .1), aes(fill=R001 ), color=NA, show.legend = FALSE) +
  scale_fill_viridis_c(option = 'E') +
  annotate("text", label="Income", x=x, y= 9.18, hjust = 0, color=clr, size=sz) +
  
  # Employment accessibility
  new_scale_fill() +
  geom_sf(data = subset(aop, P001>0) %>% rotate_data(y_add = .2), aes(fill=CMATT60 ), color=NA, show.legend = FALSE) +
  scale_fill_viridis_c(option = 'inferno') +
  annotate("text", label="Job access", x=x, y= 9.25, hjust = 0, color=clr, size=sz) +

  # public transport routes
  geom_sf(data = routes %>% rotate_data(y_add = 0.3), color="gray80", alpha=.4, show.legend = FALSE) +
  annotate("text", label="Public transport \nNetwork", x=x, y= 9.35, hjust = 0, color=clr, size=sz) +

  # public transport trips
  new_scale_color() +
  geom_point(data = gps_points_df, aes(x = x, y=y, color = shape_id), size= .8, alpha = 0.5, show.legend = FALSE) +
  # annotate("text", label="Public transport \nvehicles", x=x, y= 9.35, hjust = 0, color=clr, size=sz) +
  scale_colour_viridis(discrete = T) +
  theme_void() +
  
  # gganimate specificatons
  labs(title="Time: {frame_time}") +
  transition_time(as.POSIXct(departure_time) + 10800) +  # need to fix issue with time zone
  shadow_wake(wake_length = 0.015, alpha = FALSE) +
  ease_aes('linear') +
  coord_sf(xlim = c(-81.5, -80.8)) +
  theme(plot.background = element_rect(fill="white", color="white"),
        plot.caption = element_text(color="gray30"))



###### Save ---------------------------------------------

# save gif
anim_save(animation = temp2, "stacked_map_gtfs2gps222.gif", fps = 10, width=550, height=400)

# save mp4
anim <- animate(temp2, duration = 10, fps = 10, renderer = av_renderer(), width=550, height=400)
anim_save(animation = anim, "stacked_map_gtfs2gps222.mp4", width=550, height=400)