Scenarios

Author

R Félix

Work day vs. weekend

With the same method presented in Accessibility, let’s compare the population that can be reach from city center on a peak hour and a Sunday evening, using public transit and max 1 transfer.

Public transit - services available

First, we can have a look at the available routes on both time periods. We use tidytransit::filter_feed_by_date() function, to filter our GTFS by day and by hour.

Code 
# Example with Carris
carris_different_trips_all = length(unique(gtfs_carris$trips$trip_id)) # 53028

# wednesday morning
carris_wednesday = gtfs_carris |> filter_feed_by_date("2025-10-01")
carris_different_trips_wednesday = length(unique(carris_wednesday$trips$trip_id)) # 11504
carris_wednesday_peakh = gtfs_carris |> filter_feed_by_date("2025-10-01",
                                                            min_departure_time = "08:00:00",
                                                            max_arrival_time = "09:00:00")
carris_different_trips_wednesday_h = length(unique(carris_wednesday_peakh$trips$trip_id)) # 1352

# sunday evening
carris_sunday = gtfs_carris |> filter_feed_by_date("2025-10-05")
carris_different_trips_sunday = length(unique(carris_sunday$trips$trip_id)) # 6268
carris_sunday_peakh = gtfs_carris |> filter_feed_by_date("2025-10-05",
                                                          min_departure_time = "22:00:00",
                                                          max_arrival_time = "23:00:00")
carris_different_trips_sunday_h = length(unique(carris_sunday_peakh$trips$trip_id)) # 264

# comparision
carris_different_trips_sunday/carris_different_trips_wednesday # 54%
carris_different_trips_sunday_h/carris_different_trips_wednesday_h # 19%

For instance, for Carris (bus operator), on a Sunday there are 54% of the departures of a Wednesday, and on a Sunday evening (22h) there are 19% of the departures that are on a Wednesday morning (8h).

Transfers needed to reach all stops

By considering the PTransit stops as destinations (and not population), we can also estimate how many bus stops are reachable, in total and by number of transfers needed, for each our.

These estimates are based on GTFS data only, and do not rely on r5r.

summary(gtfs_carris)
tidygtfs object
files        agency, stops, routes, trips, stop_times, calendar, calendar_dates, shapes, transfers
agency       Carris
service      from 2025-09-01 to 2025-12-01
uses         stop_times (no frequencies)
# routes       175
# trips      53028
# stop_ids    2341
# stop_names  1139
# shapes       307
# We will use the Cais Sodré bus stop name
from_station = "Cais Sodré"
stop_times_carris = filter_stop_times(gtfs_carris, "2025-10-01", "08:00:00", "09:00:00")
trip_duration_cais = travel_times(stop_times_carris, from_station, stop_dist_check = FALSE)

nrow(trip_duration_cais) # 989

summary(trip_duration_cais$travel_time/60) # in minutos
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00   19.85   29.87   29.72   40.02   62.30

Not counting the waiting time at Cais do Sodré it is possible to reach 989 stops in the network (42%) during rush hour—within the 1 hour time frame defined.

Code 
round(prop.table(table(trip_duration_cais$transfers))*100, 1) # transfers percentage
   0    1    2    3 
22.6 55.5 19.7  2.1

Only 22% of stops are directly accessible (without bus transfers). 78% are accessible with a maximum of 1 transfer, and 2.1% of stops require at least 4 buses to reach (3 transfers).

Your turn

Try to do the same math for a sunday evening.

Mobility Peak-hour

We use the same code as in Isochrones.

# define some parameters
mode = c("SUBWAY", "BUS") # TRANSIT, BUS, SUBWAY, RAIL, CAR, FERRY, WALK, BIKE, TRAM
mode_egress = "WALK"
max_walk_time = 10 # in minutes
max_trip_duration = 90 # in minutes
time_window = 30 # in minutes
time_intervals <- seq(0, 100, 10)
departure_datetime_HP = as.POSIXct("01-10-2025 8:00:00", format = "%d-%m-%Y %H:%M:%S") # quarta-feira

# calculate travel time matrix
ttm_zer_PT_wednesday = travel_time_matrix(r5r_network = r5r_lisboa,
                              origins = BAIXA,
                              destinations = POINTS,
                              mode = mode,
                              mode_egress = mode_egress,
                              departure_datetime = departure_datetime_HP,
                              max_walk_time = max_walk_time,
                              max_trip_duration = max_trip_duration,
                              time_window = time_window,
                              max_rides = 3, # max 1 transfer
                              verbose = FALSE)

summary(ttm_zer_PT_wednesday$travel_time_p50)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00   29.00   37.00   36.24   43.00   65.00
Code 
# add coordinates of destinations to travel time matrix
ttm_zer_PT_wednesday = ttm_zer_PT_wednesday |>
  left_join(POINTS, by = c("to_id" = "id"))

# interpolate estimates to get spatially smooth result
travel_times.interp <- with(na.omit(ttm_zer_PT_wednesday), interp(lon, lat, travel_time_p50)) |>
  with(cbind(travel_time=as.vector(z),  # Column-major order
             x=rep(x, times=length(y)),
             y=rep(y, each=length(x)))) |>
  as.data.frame() |> na.omit()

# find isochrone's bounding box to crop the map below
bb_x <- c(min(travel_times.interp$x), max(travel_times.interp$x))
bb_y <- c(min(travel_times.interp$y), max(travel_times.interp$y))
# plot
plotHP = ggplot(travel_times.interp) +
  geom_contour_filled(aes(x = x, y = y, z = travel_time), alpha = .7) +
  geom_sf(data = road_network_base, color = "gray55", lwd = 0.5, alpha = 0.4) +
  geom_sf(data = city_limit, fill = "transparent", color = "grey30") +
  geom_point(aes(x = lon, y = lat, color = 'Baixa'), data = BAIXA) +
  scale_fill_viridis_d(direction = -1, option = 'B') +
  scale_color_manual(values = c('Baixa' = 'black')) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  coord_sf(xlim = bb_x, ylim = bb_y) +
  labs(
    title = "Reaching from Baixa (Carris + Metro)",
    subtitle = "at 8am wednesday - 1 transf max",
    fill = "Trip duration \n[min]",
    color = ''
  ) +
  theme_minimal() +
  theme(axis.title = element_blank())

plotHP

Mobility Off peak-hour

departure_datetime_OFF = as.POSIXct("05-10-2025 22:00:00", format = "%d-%m-%Y %H:%M:%S") # sunday evening

# calculate travel time matrix
ttm_zer_PT_sunday = travel_time_matrix(r5r_network = r5r_lisboa,
                              origins = BAIXA,
                              destinations = POINTS,
                              mode = mode,
                              mode_egress = mode_egress,
                              departure_datetime = departure_datetime_OFF, # new date
                              max_walk_time = max_walk_time,
                              max_trip_duration = max_trip_duration,
                              time_window = time_window,
                              max_rides = 2, # max 1 transfer
                              verbose = FALSE)

summary(ttm_zer_PT_sunday$travel_time_p50)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00   32.00   41.00   40.81   50.00   79.00
Code 
# add coordinates of destinations to travel time matrix
ttm_zer_PT_sunday = ttm_zer_PT_sunday |>
  left_join(POINTS, by = c("to_id" = "id"))

# interpolate estimates to get spatially smooth result
travel_times.interp <- with(na.omit(ttm_zer_PT_sunday), interp(lon, lat, travel_time_p50)) |>
  with(cbind(travel_time=as.vector(z),  # Column-major order
             x=rep(x, times=length(y)),
             y=rep(y, each=length(x)))) |>
  as.data.frame() |> na.omit()

# find isochrone's bounding box to crop the map below
bb_x <- c(min(travel_times.interp$x), max(travel_times.interp$x))
bb_y <- c(min(travel_times.interp$y), max(travel_times.interp$y))
# plot
plotOFF = ggplot(travel_times.interp) +
  geom_contour_filled(aes(x = x, y = y, z = travel_time), alpha = .7) +
  geom_sf(data = road_network_base, color = "gray55", lwd = 0.5, alpha = 0.4) +
  geom_sf(data = city_limit, fill = "transparent", color = "grey30") +
  geom_point(aes(x = lon, y = lat, color = 'Baixa'), data = BAIXA) +
  scale_fill_viridis_d(direction = -1, option = 'B') +
  scale_color_manual(values = c('Baixa' = 'black')) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  coord_sf(xlim = bb_x, ylim = bb_y) +
  labs(
    title = "Reaching from Baixa (Carris + Metro)",
    subtitle = "at 10pm sunday - 1 transf max",
    fill = "Trip duration \n[min]",
    color = ''
  ) +
  theme_minimal() +
  theme(axis.title = element_blank())

plotOFF

Compare population

Based on Census 2021 (INE 2022), we can count how many people are reach.

Code 
poplisboa = sum(POINTS$residents) #
100* sum(ttm_zer_PT_wednesday$residents[ttm_zer_PT_wednesday$travel_time_p50 <= 15]) / poplisboa # 2.8%
100* sum(ttm_zer_PT_wednesday$residents[ttm_zer_PT_wednesday$travel_time_p50 <= 30]) / poplisboa # 4.4%
100* sum(ttm_zer_PT_wednesday$residents[ttm_zer_PT_wednesday$travel_time_p50 <= 45]) / poplisboa # 90.7%
100* sum(ttm_zer_PT_wednesday$residents[ttm_zer_PT_wednesday$travel_time_p50 <= 60]) / poplisboa # 98.9%

100* sum(ttm_zer_PT_sunday$residents[ttm_zer_PT_sunday$travel_time_p50 <= 15]) / poplisboa # 1.3%
100* sum(ttm_zer_PT_sunday$residents[ttm_zer_PT_sunday$travel_time_p50 <= 30]) / poplisboa # 29.6%
100* sum(ttm_zer_PT_sunday$residents[ttm_zer_PT_sunday$travel_time_p50 <= 45]) / poplisboa # 73.5%
100* sum(ttm_zer_PT_sunday$residents[ttm_zer_PT_sunday$travel_time_p50 <= 60]) / poplisboa # 93.5%
Up to… Wednesday peak hour Sunday evening
30 min 42.3% 29.6%
45 min 90.7% 73.5%
60 min 98.9% 93.5%

New metro station or metro line

Let’s create a new metro line or a new station in your city.

New pink line (expansion in Lisbon)

Methodology

For that, an easy graphic interface software is the PTV Lines.

  1. Login into PTV Lines (upper right corner) with the credentials that were sent to you.
  2. Create a new scenario, and import an existing subway gtfs (.zip) - max 50 lines!1
  3. Edit that scenario by creating stops and lines.
  4. Export the GTFS file.
  5. Duplicate your r5r folder (baseline scenario) and delete the network.dat file inside. Rename it as ..._ver1 and replace the gtfs.zip by the one you just exported.
  6. Build a new network with that folder path.

Create a metro line with PTV Lines

First, create new station points (S) and name them. As many as you need.

Then, create a new line (L), name it something like “Expansion”, and create a new route (R). Name it 0 (you can edit later) and draw that line roughly by selecting an existing start station and connecting to your new stations.

Add trips to it (T), with a proper headway.

After defining the weekday and weekend trips, go to the Outbound route and “Create opposite direction”, so you can have both directions now. Rename the line routes, something as “Start station - End station”.

Extend a metro line with PTV Lines

This tutorial assumes that the stations for which the line will be extended to already exist (if not, see how to create new stations at Create a metro line with PTV Lines).

The extension is performed over a route. To access one, click on a transit line, and then on one of its routes.

On the map, hover the Start/End buttons and click it to enable the extension mode.

Once the extension mode is activated, your cursor will turn into a pen icon and you just have to click on the station you want to expand your line to.

Routing modes

Be aware of the several modes available on the button that is now displayed (automatic routing mode, street routing mode, straight line mode or exclusive PT routing mode). Make sure you use the one that best suits your use case.

Once you click on the station, you should see the new station in the route window, on the left part of the screen. You should adjust the run and dwell times accordingly.

Visualizing time diagram

Go to Connections, and click on a connecting station to visualize the time diagram of connections with other routes.

Generate GTFS

After all your edits, go back to the dashboard, and click on “Export supply”. Export as GTFS (.zip)

Your scenario will be saved in your dashboard. You can edit or duplicate them later.

Your task

Run some routing and accessibility analysis, and compare with the baseline scenario.

Consider comparing the estimated cost of your new infrastructure and the mobility gains.

See Pereira and Herszenhut (2023) Chapter 6 for an example.

If you do multi-modal analysis, you may need to use GTFShift::unify() to merge the new gtfs with other modes and create new transfers.

If you have more than 50 lines, consider filtering your GTFS by area or by mode, and use that one as the baseline scenario, and the improved one (PTV lines) as you comparison scenario.

Code example to filter a gtfs by mode 
# example to filter a gtfs by mode = TRAM
gtfs_all_modes = gtfstools::read_gtfs("data/gtfs_all.zip")
gtfs_tram = gtfstools::filter_by_route_type(gtfs_all_modes,
                                            route_type = 0) # see the documentation for the other modes' code
gtfstools::write_gtfs("data/gtfs_tram-zip")

References

INE. 2022. “Censos 2021- XVI Recenseamento Geral da População. VI Recenseamento Geral da Habitação.” Lisboa: Instituto National de Estatística. https://censos.ine.pt/xurl/pub/65586079.
Pereira, Rafael HM, and Daniel Herszenhut. 2023. Introduction to Urban Accessibility: A Practical Guide with r. Instituto de Pesquisa Econômica Aplicada (Ipea). https://ipeagit.github.io/intro_access_book/.

Footnotes

  1. Tip: you may need to filter by a mode (tram, or subway) to avoid having more than 50 lines.↩︎