septa_bus <- read_gtfs("google_bus.zip")
septa_rail <- read_gtfs("google_rail.zip")

bge_bus <- get_stop_frequency(
  septa_bus,
  start_time = "07:00:00",
  end_time = "09:00:00",
  service_ids = NULL,
  by_route = TRUE
) %>% filter(route_id == '23' | route_id == '16' | route_id == '53' | route_id == '56' | route_id == 'H' | route_id == 'XH')

septa_stops <- septa_bus[["stops"]]

bge_bus_with_stops <- merge(bge_bus, septa_stops, by = 'stop_id')

makePoints <- function(stop_lat, stop_lon) {
  return(st_point(c(stop_lat, stop_lon)))
}

geom_bus <- bge_bus_with_stops %>%
  select(stop_lat, stop_lon) %>% 
  pmap(makePoints) %>%
  st_as_sfc(crs = 4326) %>%
  {cbind(bge_bus_with_stops, geom = .)} %>%
  st_sf()


#plot headways
geom_bus <- geom_bus %>%
  arrange(mean_headway) %>%
  mutate(mean_headway=mean_headway/60)

geom_bus <- geom_bus %>% 
  filter(stop_lon < -75.14) %>% 
  filter(stop_lon > -75.16) %>% 
  filter(stop_lat > 40) %>% 
  filter(stop_lat < 40.02) %>%
  filter(mean_headway < 200)

Bus Headways Around BGE

Greater diameter correpsonds with greater headways or wait times.

pal <- colorFactor(
  palette = c("darkcyan", "orchid", "orange", "green", "midnightblue", "darkslategrey"),
  domain = geom_bus$route_id
)

# make interactive map
leaflet(geom_bus) %>%
  fitBounds(-75.167015,39.995993,-75.131325,40.021987) %>%
  addTiles() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircles(lng = ~stop_lon, lat = ~stop_lat, radius = ~mean_headway, weight = 1.5,
             color = ~pal(route_id), popup = ~paste("Route", route_id, "mean headway (min):", mean_headway)) %>%
  addLegend(position = "bottomleft", pal = pal, values = ~route_id,
            title = "")