Thursday, March 19, 2020

Creating an animated gif of public transport networks using GTFS data



Creating an animated gif of public transport networks using GTFS data and the gtfs2gps package

In this gist we show with a reproducible example how to create an animation of public transport networks using GTFS data in R. We use a few packages to do this. One of the core packages here is the new gtfs2gps. The gtfs2gps package converts public transport data in GTFS format to GPS-like records in a data.frame/data.table, which we will be using to create a .gif with the gganimate package.

Step 1 - Convert GTFS to GPS-like records

The first step is to process a GTFS.zip file. The function gtfs2gps{gtfs2gps} interpolates the space-time position of each vehicle in each trip considering the network distance and average speed between stops. The output is a data.table where each row represents the timestamp of each vehicle at a given spatial resolution. In this example, we use a sample of the public transport network of Sao Paulo (Brazil) mapped every 30 meters.


# load libraries
library(gtfs2gps)
library(geobr)
library(ggplot2)
library(gganimate)
library(data.table)
library(ggthemes)
library(sf)
library(viridis)
library(sfheaders)



###### 1. process public transport data  ------------------

# get local GTFS.zip file
gtfs_zip <- system.file("extdata/saopaulo.zip", package="gtfs2gps" )

# read gtfs
gtfs_dt <- gtfs2gps::read_gtfs(gtfs_zip)

# filter transport services on weekdays
gtfs_dt <- gtfs2gps::filter_week_days(gtfs_dt)

# get transport network as sf object
shapes_sf <- gtfs2gps::gtfs_shapes_as_sf(gtfs_dt)

# Convert GTFS data into a data.table with GPS-like records
gps_dt <- gtfs2gps::gtfs2gps(gtfs_dt, spatial_resolution = 30, progress = T, parallel = T)
head(gps_dt)

# subset time interval. Get only vehicle positions between 7am and 7:30am
gps_dt2 <- gps_dt[ between(departure_time, as.ITime("07:00:"), as.ITime("07:31"))]

# remove observations with 'too high speed' (in km/h)
gps_dt2 <- gps_dt2[speed < 20 , ]


# Convert "GPS" points into sf
gps_sf <- sfheaders::sf_point(gps_dt2, x = "shape_pt_lon" , y = "shape_pt_lat", keep = T)
sf::st_crs(gps_sf) <- 4326
head(gps_sf)

Step 2 - Create gif

Once you have created the GPS-like dataset of the public transport network, creating an animation gif is simple as this:


###### 2. gif with blank background [fully reproducible] ------------------

# download Sao Paulo border for spatial context
sao_paulo <- geobr::read_municipality(code_muni = 3550308)

anim <- ggplot() +
          geom_sf(data = sao_paulo, color='gray50', fill=NA) +
          geom_sf(data = shapes_sf, color='gray90', size=0.01) +
          geom_point(data = gps_dt2, aes(x = shape_pt_lon, y=shape_pt_lat, colour = trip_id), size=1.5, alpha = 0.5, show.legend = FALSE) +
          scale_colour_viridis( discrete = T ) +
          coord_sf(ylim = c(min(gps_dt2$shape_pt_lat), max(gps_dt2$shape_pt_lat)) ) +
          theme_map() +
  
          # 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')


# save gif
anim_save(animation = anim, "./tests_rafa/spo_gtfs_sample_20fps__GIST.gif", fps = 20)

Extra

In case you'd like to add more context info to the gif, you can use rasters as spatial tiles. I've added here some background info using a mapbox tile created by Kaue Braga for the Access to Opportunities Project. Thanks Kaue.


###### 3. gif tile background ------------------

# read map tile for spatial context
map_tiles <- readr::read_rds("L:/Proj_acess_oport/data/maptiles_crop/2019/mapbox/maptile_crop_mapbox_spo_2019.rds")

# convert gps data to UTM
gps_sf_utm <- sf::st_transform(gps_sf, 3857)
gps_dt_utm <- sfheaders::sf_to_df( gps_sf_utm, fill = TRUE )

anim2 <- ggplot() +
          geom_raster(data = map_tiles, aes(x, y, fill = hex), alpha = 1) +
          scale_fill_identity() +
          geom_sf(data= st_transform(sao_paulo, 3857), color='gray50', fill=NA) +
          geom_point(data = gps_dt_utm, aes(x = x, y=y, colour = trip_id), size=1.5, alpha = 0.6, show.legend = FALSE) +
          scale_colour_viridis(discrete = T) +
          coord_sf(ylim = c(min(gps_dt_utm$y), max(gps_dt_utm$y)) ) +
          theme_map() + 
          
          # 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')
        

# save gif
anim_save(animation = anim2, "./tests_rafa/spo_gtfs_sample_20fps_tile__GIST.gif", fps = 20)