Skip to contents

From raw GPS to a realized GTFS feed

Source: vignettes/articles/pipeline.Rmd
pipeline.Rmd

This article follows transit observations from raw vehicle positions to a static GTFS feed of the service that actually ran, across four packages:

Step Package Functions
Parse GTFS-Realtime protobuf gtfsrealtime read_gtfsrt_positions()
Infer trips and stop visits from coordinates gps2gtfs g2g_stops_from_gtfs(), g2g_terminals_from_gtfs(), g2g_extract_trips_and_stop_times(), g2g_shapes_from_trips()
Assemble a compliant feed gtfsrt2static rt2s_events_from_stop_times(), rt2s_assemble()
Write, read and analyse GTFS gtfsio, gtfstools export_gtfs(), read_gtfs(), …

The simple example scaffolds a feed from one day of a bus’s GPS pings with no planned feed. The advanced example matches GTFS-Realtime positions to an operator’s planned feed and compares planned with realized times.

Setup

gps2gtfs is not on CRAN. Install it, and gtfsrt2static, from GitHub:

# install.packages("pak")
pak::pak(c("e-kotov/gps2gtfs", "e-kotov/gtfsrt2static", "gtfsrealtime", "gtfstools"))
library(gps2gtfs)
library(gtfsrt2static)
library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following object is masked from 'package:base':
#> 
#>     %notin%

Vehicle Positions are GPS pings

GTFS-Realtime Vehicle Positions report where each vehicle is. gtfsrealtime parses them, here from one bundled snapshot of the New York City bus feed:

vp <- gtfsrealtime::read_gtfsrt_positions(
  system.file("nyc-vehicle-positions.pb.bz2", package = "gtfsrealtime"),
  timezone = "America/New_York"
)
vp <- as.data.table(vp)
vp[1:3, .(vehicle_id, trip_id, route_id, latitude, longitude, timestamp, speed)]
#>       vehicle_id                             trip_id route_id latitude
#>           <char>                              <char>   <char>    <num>
#> 1: MTA NYCT_9771    MV_A6-Weekday-SDon-110100_M5_527       M4 40.78720
#> 2: MTA NYCT_8440   JA_A6-Weekday-SDon-110800_Q17_334      Q17 40.74109
#> 3: MTA NYCT_9770 GH_A6-Weekday-SDon-111400_BX238_640     BX28 40.87092
#>    longitude           timestamp speed
#>        <num>              <POSc> <num>
#> 1: -73.95419 2026-01-21 18:58:29    NA
#> 2: -73.78681 2026-01-21 18:58:24    NA
#> 3: -73.84709 2026-01-21 18:58:32    NA

These are the column names gps2gtfs expects, so the reader’s output goes into the extractor with no adapter. Columns it does not use, such as trip_id, route_id and stop_id, are passed through. This feed leaves speed empty. gps2gtfs detects dwell from pings with speed == 0, so with such a feed it warns and treats every ping as moving.

One snapshot holds one ping per vehicle:

nrow(vp)
#> [1] 3018
uniqueN(vp$vehicle_id)
#> [1] 3018

A trajectory needs many polls of the same vehicle. Archive the endpoint over a day with rt2s_collect(), then read the daily ZIP with the same function. The examples below use trajectories that ship with gps2gtfs instead.

Simple example: scaffold a feed from GPS pings

gps2gtfs bundles one day of GPS pings from a bus on route 654 between Kandy and Digana, Sri Lanka, with the route’s stops and terminals. There is no planned GTFS feed for it.

Extract trips, stop visits and shapes 

result <- g2g_extract_trips_and_stop_times(
  gps_data = g2g_data_gps,
  terminals_data = g2g_data_terminals,
  stops_data = g2g_data_stops,
  terminals_buffer_radius = 50,
  stops_buffer_radius = 30,
  stops_extended_buffer_radius = 50,
  stop_direction_map = c("Kandy-Digana" = "BT01", "Digana-Kandy" = "BT02"),
  return_trajectory = TRUE,
  diagnostics_warn = FALSE
)
#> Starting Pipeline for extracting Trip and Bus Stop Data using backend: rust
#> [INFO] Removed 56 duplicated (vehicle_id, timestamp) row(s).
#> [INFO] Auto-detected local UTM projection EPSG:32644 (Zone 44N)
#> Pipeline finished successfully!
result$trips[, .(trip_id, start_terminal, end_terminal, start_time, end_time)]
#>    trip_id start_terminal end_terminal          start_time            end_time
#>      <int>         <char>       <char>              <POSc>              <POSc>
#> 1:       1           BT02         BT01 2022-07-01 06:52:06 2022-07-01 07:38:17
#> 2:       2           BT01         BT02 2022-07-01 07:56:29 2022-07-01 08:50:15
#> 3:       3           BT02         BT01 2022-07-01 10:46:59 2022-07-01 11:42:31
#> 4:       4           BT01         BT02 2022-07-01 12:49:00 2022-07-01 13:44:14
#> 5:       5           BT02         BT01 2022-07-01 15:19:18 2022-07-01 16:13:18
#> 6:       6           BT01         BT02 2022-07-01 16:51:39 2022-07-01 18:02:08

shapes <- g2g_shapes_from_trips(result$trajectory)

The gps2gtfs Get started article explains these arguments and the coverage diagnostics.

Convert to observed stop events

rt2s_events_from_stop_times() turns the stop visits into the observed stop events table that every gtfsrt2static assembler consumes. gps2gtfs models one route per run, so the route identity is stamped on. shape_ref_prefix must match the shape_id_prefix of g2g_shapes_from_trips() ("SHP_" by default) so each trip links to its driven shape.

events <- rt2s_events_from_stop_times(
  result$stop_times,
  route_ref = "654",
  source = "gps",
  shape_ref_prefix = "SHP_"
)
events[1:3, .(trip_ref, route_ref, shape_ref, service_date, stop_ref,
              stop_sequence, arrival_time, provenance)]
#> Key: <service_date, trip_ref, stop_sequence>
#>          trip_ref route_ref shape_ref service_date stop_ref stop_sequence
#>            <char>    <char>    <char>       <Date>   <char>         <int>
#> 1: g2g_20220701_1       654     SHP_1   2022-07-01      201             1
#> 2: g2g_20220701_1       654     SHP_1   2022-07-01      202             2
#> 3: g2g_20220701_1       654     SHP_1   2022-07-01      203             3
#>           arrival_time provenance
#>                 <POSc>     <char>
#> 1: 2022-07-01 06:57:49   observed
#> 2: 2022-07-01 07:00:12   observed
#> 3: 2022-07-01 07:04:34   observed

Each trip has no official identifier, so it gets a synthetic one, unique across service dates.

Scaffold the feed

With no planned feed, rt2s_assemble() scaffolds every required file. Agency metadata and stop coordinates must be supplied, and strict = TRUE makes any gap an error, so a feed that returns is free of known spec-required gaps. The bundled pings are stored in UTC, so UTC is the feed’s timezone here; with real data, give the agency’s own.

stops <- data.frame(
  stop_id = as.character(g2g_data_stops$stop_id),
  stop_name = g2g_data_stops$address,
  stop_lat = g2g_data_stops$latitude,
  stop_lon = g2g_data_stops$longitude
)

feed <- rt2s_assemble(
  events,
  agency = list(name = "Route 654 demo", url = "https://example.org",
                timezone = "UTC"),
  stops = stops,
  shapes = shapes,
  strict = TRUE
)
#> [INFO] route_type not given; scaffolding routes as 3 (bus).
names(feed)
#> [1] "agency"         "stops"          "routes"         "trips"         
#> [5] "stop_times"     "calendar_dates" "feed_info"      "shapes"
feed$trips[]
#>    route_id   service_id        trip_id direction_id shape_id
#>      <char>       <char>         <char>        <int>   <char>
#> 1:      654 SVC_20220701 g2g_20220701_1            1    SHP_1
#> 2:      654 SVC_20220701 g2g_20220701_2            0    SHP_2
#> 3:      654 SVC_20220701 g2g_20220701_3            1    SHP_3
#> 4:      654 SVC_20220701 g2g_20220701_4            0    SHP_4
#> 5:      654 SVC_20220701 g2g_20220701_5            1    SHP_5
#> 6:      654 SVC_20220701 g2g_20220701_6            0    SHP_6
rt2s_publishable(feed)$publishable
#> [1] TRUE

Write it and analyse it with gtfstools 

path <- tempfile(fileext = ".zip")
gtfsio::export_gtfs(feed, path)

gtfs <- gtfstools::read_gtfs(path)
merge(
  gtfstools::get_trip_duration(gtfs, unit = "min"),
  gtfstools::get_trip_speed(gtfs, unit = "km/h")[, .(trip_id, speed)],
  by = "trip_id"
)
#> Key: <trip_id>
#>           trip_id duration    speed
#>            <char>    <num>    <num>
#> 1: g2g_20220701_1 35.08333 27.08431
#> 2: g2g_20220701_2 44.40000 21.69604
#> 3: g2g_20220701_3 45.88333 20.88370
#> 4: g2g_20220701_4 45.66667 21.14421
#> 5: g2g_20220701_5 44.05000 21.83688
#> 6: g2g_20220701_6 56.61667 17.04455

get_trip_speed() measures each trip along its shape, which here is the road the bus actually drove. The same zip opens in tidytransit and other GTFS tools.

Advanced example: realized service against a planned feed

When the operator publishes a planned feed, and the Vehicle Positions carry its trip_id, the realized feed can keep the official identifiers. It then joins directly to the planned feed. The example uses a small two-stop route, so that every step fits on screen.

The planned feed

A planned feed with route R1, terminals A and B, two intermediate stops, and two trips: CS_1 from A to B at 08:00, and CS_2 back at 08:30.

lat <- seq(7.290, 7.320, length.out = 4)
lon <- seq(80.630, 80.660, length.out = 4)
planned_times <- c("08:00:00", "08:04:00", "08:08:00", "08:12:00",
                   "08:30:00", "08:34:00", "08:38:00", "08:42:00")

planned <- gtfsio::new_gtfs(list(
  agency = data.table(agency_id = "DT", agency_name = "Demo Transit",
                      agency_url = "https://example.org",
                      agency_timezone = "UTC"),
  routes = data.table(route_id = "R1", agency_id = "DT",
                      route_short_name = "1", route_type = 3L),
  calendar = data.table(service_id = "daily",
                        monday = 1L, tuesday = 1L, wednesday = 1L,
                        thursday = 1L, friday = 1L, saturday = 1L,
                        sunday = 1L,
                        start_date = 20260601L, end_date = 20260630L),
  stops = data.table(stop_id = c("A", "S1", "S2", "B"),
                     stop_name = c("Terminal A", "Stop 1", "Stop 2",
                                   "Terminal B"),
                     stop_lat = lat, stop_lon = lon),
  trips = data.table(route_id = "R1", service_id = "daily",
                     trip_id = c("CS_1", "CS_2"), direction_id = 0:1),
  stop_times = data.table(
    trip_id = rep(c("CS_1", "CS_2"), each = 4),
    arrival_time = planned_times,
    departure_time = planned_times,
    stop_id = c("A", "S1", "S2", "B", "B", "S2", "S1", "A"),
    stop_sequence = rep(1:4, 2)
  )
))

planned_zip <- tempfile(fileext = ".zip")
gtfsio::export_gtfs(planned, planned_zip)

Positions annotated with the planned trip

These are the parsed Vehicle Positions of the bus on 6 June 2026, polled every five minutes: one run from A to B, a layover ping without a trip, and the run back. The trip_id column is what read_gtfsrt_positions() returns when the feed annotates positions with the planned trip.

t0 <- as.POSIXct("2026-06-06 08:00:00", tz = "UTC")
positions <- data.table(
  vehicle_id = "7482",
  latitude = c(lat, 7.321, rev(lat)),
  longitude = c(lon, 80.661, rev(lon)),
  timestamp = t0 + c(0, 300, 600, 900, 1200, 1800, 2100, 2400, 2700),
  speed = c(0, 20, 20, 0, 0, 0, 20, 20, 0),
  trip_id = c(rep("CS_1", 4), NA, rep("CS_2", 4))
)

Stops and terminals from the planned feed

gps2gtfs derives both tables from the planned feed. The stops table holds the route’s non-terminal stops, labelled by the terminal their direction starts from, so no direction map is needed:

terminals <- g2g_terminals_from_gtfs(planned, route_id = "R1")
stops <- g2g_stops_from_gtfs(planned, route_id = "R1")
terminals
#>    terminal_id latitude longitude
#>         <char>    <num>     <num>
#> 1:           A     7.29     80.63
#> 2:           B     7.32     80.66
stops
#> Key: <direction, stop_id>
#>    stop_id latitude longitude direction
#>     <char>    <num>     <num>    <char>
#> 1:      S1     7.30     80.64         A
#> 2:      S2     7.31     80.65         A
#> 3:      S1     7.30     80.64         B
#> 4:      S2     7.31     80.65         B

Extract along the planned trips

trip_col = "trip_id" takes the trips from the annotation instead of inferring them from terminal crossings (the fast path). The annotation is kept as provided_trip_id.

result <- g2g_extract_trips_and_stop_times(
  gps_data = positions,
  terminals_data = terminals,
  stops_data = stops,
  terminals_buffer_radius = 100,
  stops_buffer_radius = 100,
  stops_extended_buffer_radius = 150,
  trip_col = "trip_id",
  diagnostics_warn = FALSE
)
#> Starting Pipeline for extracting Trip and Bus Stop Data using backend: rust
#> [INFO] Auto-detected local UTM projection EPSG:32644 (Zone 44N)
#> [INFO] Using supplied trip identities from column(s) 'trip_id' (fast path).
#> Pipeline finished successfully!
result$stop_times[, .(trip_id, provided_trip_id, stop_id, arrival_time)]
#>    trip_id provided_trip_id stop_id        arrival_time
#>      <int>           <char>  <char>              <POSc>
#> 1:       1             CS_1      S1 2026-06-06 08:05:00
#> 2:       1             CS_1      S2 2026-06-06 08:10:00
#> 3:       2             CS_2      S2 2026-06-06 08:35:00
#> 4:       2             CS_2      S1 2026-06-06 08:40:00

Assemble against the baseline

rt2s_events_from_stop_times() uses a provided_trip_id column as the trip identity when one is present. Given the planned feed as baseline, rt2s_assemble() keeps its agency, routes and stops, matches the observed trips by trip_id, and replaces stop_times.txt with the observed times for the service day:

events <- rt2s_events_from_stop_times(result$stop_times)
unique(events$trip_ref)
#> [1] "CS_1" "CS_2"

realized <- rt2s_assemble(events, baseline = planned_zip)
realized$trips
#>    route_id   service_id trip_id direction_id
#>      <char>       <char>  <char>        <int>
#> 1:       R1 SVC_20260606    CS_1            0
#> 2:       R1 SVC_20260606    CS_2            1
realized$stop_times
#>    trip_id arrival_time departure_time stop_id stop_sequence
#>     <char>       <char>         <char>  <char>         <int>
#> 1:    CS_1     08:05:00       08:05:00      S1             1
#> 2:    CS_1     08:10:00       08:10:00      S2             2
#> 3:    CS_2     08:35:00       08:35:00      S2             1
#> 4:    CS_2     08:40:00       08:40:00      S1             2

Compare planned and realized

Both feeds share trip and stop identifiers, so gtfstools reads both and a join gives the delay at each stop. Join on trip_id and stop_id, not on stop_sequence: for gps2gtfs input the realized stop_sequence counts the visits in observed order, while the planned one numbers every stop of the pattern, including the terminals. (In a loop that serves a stop twice, the stop alone does not identify the visit.)

realized_zip <- tempfile(fileext = ".zip")
gtfsio::export_gtfs(realized, realized_zip)

secs <- function(path) {
  g <- gtfstools::read_gtfs(path)
  gtfstools::convert_time_to_seconds(g, file = "stop_times")$stop_times
}
comparison <- merge(
  secs(planned_zip)[, .(trip_id, stop_id, planned = arrival_time_secs)],
  secs(realized_zip)[, .(trip_id, stop_id, realized = arrival_time_secs)],
  by = c("trip_id", "stop_id")
)
comparison[, delay_min := (realized - planned) / 60]
comparison
#> Key: <trip_id, stop_id>
#>    trip_id stop_id planned realized delay_min
#>     <char>  <char>   <int>    <int>     <num>
#> 1:    CS_1      S1   29040    29100         1
#> 2:    CS_1      S2   29280    29400         2
#> 3:    CS_2      S1   31080    31200         2
#> 4:    CS_2      S2   30840    30900         1

Each run was one minute late at its first intermediate stop and two at its second.

Where next