madison-metro/madison-metro.R

library(tidyverse)
library(influxdbclient)
library(glue)
library(ggmap)
library(sf)

# parameters needed to make connection to Database
token <- substr(read_file(file = 'api_keys/influxdb_madison-metro'), 1, 88)
org <- "32b7fde0efd8a3b3"
bucket <- "metro_vehicles"
  

days <- 7

influx_connection <- InfluxDBClient$new(url = "https://influxdb.dendroalsia.net",
                                token = token,
                                org = org)
#---
# Fields to query
fields <- c("des", "spd", "pdist", "lon", "lat", "dly", "origtatripno")

# An empty list to store results for each field
results <- vector("list", length(fields))

# Loop through each field, get data, and coerce types if needed
for (i in seq_along(fields)) {
  field <- fields[i]
  
  query_string <- glue('from(bucket: "{bucket}") ',
                       '|> range(start: -{days}d) ',
                       '|> filter(fn: (r) => r["_measurement"] == "vehicle_data")',
                       '|> filter(fn: (r) => r["_field"] == "{field}")')
  
  data <- influx_connection$query(query_string)
  
  # Ensure the columns are coerced to consistent types 
  # (Optionally add coercion based on your expected types)
  data <- bind_rows(data) %>%
    mutate(value = as.character(`_value`),
           field = `_field`) %>% 
    select(time, rt, pid, vid, value, field)
  
  results[[i]] <- data
}

# Bind all results together
metro_raw <- bind_rows(results)

metro_raw <- pivot_wider(metro_raw, values_from = value, names_from = field) %>%
  distinct(pid, vid, lat, lon, spd, .keep_all = TRUE)

routes_categorized <- read_csv(file = "routes_categorized.csv", col_types = "cc")

metro_data <- metro_raw %>%
  mutate(time = with_tz(time, "America/Chicago"),
         spd = as.double(spd),
         pdist = as.double(pdist),
         lon = as.double(lon),
         lat = as.double(lat)) %>%
  mutate(date = date(time)) %>%
  group_by(pid, vid) %>%
  arrange(time) %>%
  mutate(pdist_lag = lag(pdist),
         time_lag = lag(time)) %>%
  ungroup() %>%
  mutate(spd_calc = case_when(pdist_lag > pdist ~ NA,
                              pdist_lag <= pdist ~ (pdist - pdist_lag)/as.double(difftime(time, time_lag, units = "hours"))/5280)) %>%
  left_join(routes_categorized, by = "pid")


bucket_feet <- 500

lat_round <- bucket_feet/364481.35
lon_round <- bucket_feet/267203.05

metro_summary <- metro_data %>%
  mutate(lat_bucket = round(lat / lat_round) * lat_round,
         lon_bucket = round(lon / lon_round) * lon_round) %>%
  group_by(rt, name, pid, lat_bucket, lon_bucket) %>%
  summarise(spd = median(spd, na.rm = TRUE),
            spd_calc = median(spd_calc, na.rm = TRUE),
            pdist = median(pdist),
            trip_count = length(unique(origtatripno)))
  
metro_data_sf <- st_as_sf(metro_data %>% filter(!is.na(lon), !is.na(lat)), coords = c("lon", "lat"), remove = FALSE)
metro_summary_sf <- st_as_sf(metro_summary %>% filter(!is.na(lon_bucket), !is.na(lat_bucket)), coords = c("lon_bucket", "lat_bucket"), remove = FALSE)

metro_segments <- metro_summary %>%
  group_by(rt, pid) %>%
  arrange(pdist) %>%
  mutate(lat_bucket_lag = lag(lat_bucket),
         lon_bucket_lag = lag(lon_bucket)) %>%
  filter(!is.na(lat_bucket) & !is.na(lon_bucket) & !is.na(lat_bucket_lag) & !is.na(lon_bucket_lag)) %>%
  mutate(
    geometry = pmap(list(lat_bucket, lon_bucket, lat_bucket_lag, lon_bucket_lag), 
                    ~st_linestring(matrix(c(..2, ..1, ..4, ..3), ncol = 2, byrow = TRUE)))
  ) %>%
  st_as_sf(sf_column_name = "geometry") %>%
  group_by(rt, name, lat_bucket, lon_bucket) %>%
  summarise(spd_calc = weighted.mean(spd_calc, trip_count))

# get counts of routes
route_counts <- metro_data %>% group_by(pid, rt, des) %>% summarise(route_count = length(unique(origtatripno)))

# make charts
ggplot(data = metro_summary %>% filter(pid %in% (routes_categorized %>% filter(name %in% c("B_North", "B_South")) %>% pull (pid))),
       aes(x = pdist,
           y = spd_calc)) +
  geom_point() +
  geom_smooth() +
  facet_grid(name ~ .)

ggplot(data = metro_summary %>% filter(!is.na(name)),
       aes(x = name,
           y = spd_calc)) +
  geom_boxplot()

register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, 36))

bbox <- c(left = min(metro_data$lon, na.rm = TRUE),
          bottom = min(metro_data$lat, na.rm = TRUE),
          right = max(metro_data$lon, na.rm = TRUE),
          top = max(metro_data$lat, na.rm = TRUE))

#get basemap
basemap <- get_stadiamap(bbox = bbox, zoom = 13, maptype = "stamen_toner_lite")

quantile(metro_segments %>% filter(name %in% c("A_West")) %>% pull(spd_calc), c(0,0.25, 0.5, 0.75, 1), na.rm = TRUE)

for (route in unique(routes_categorized$name)){
  route_focus <- routes_categorized %>% filter(name == route) %>% pull(pid)
  ggmap(basemap) +
    labs(title = paste0("Metro Route Speed - ", route),
         subtitle = paste0("averaged between ",
                           sum(route_counts %>% filter(pid %in% route_focus) %>% pull(route_count)),
                           " bus trips - ",
                           min(date(metro_data$time)),
                           " to ",
                           max(date(metro_data$time))),
         x = NULL,
         y = NULL) +
    theme(axis.text=element_blank(),
          axis.ticks=element_blank(),
          plot.caption = element_text(color = "grey")) +
    geom_tile(data = metro_summary %>% 
                filter(name %in% route) %>% 
                group_by(lon_bucket, lat_bucket) %>% 
                summarise(spd_calc = weighted.mean(spd_calc, trip_count)),
            inherit.aes = FALSE,
            aes(x = lon_bucket,
                y = lat_bucket,
                fill = spd_calc,
                height = lat_round,
                width = lon_round)) +
    scale_fill_distiller(palette = "RdYlGn", direction = "reverse", limits = c(0,70), name = "Average speed or segment\n(calculated with locations, not reported speed)")
  ggsave(file = paste0("figures/",
                       route,
                       "_map.pdf"),
         title = paste0("Metro Route Speed - ", route),
         device = pdf,
         height = 8.5,
         width = 11,
         units = "in",
         create.dir = TRUE)

  ggplot(data = metro_data %>% filter(name %in% route)) +
    geom_boxplot(aes(x = date,
                    y = spd_calc,
                    group = date))
  ggsave(file = paste0("figures/",
                       route,
                       "_date.pdf"),
         title = paste0("Metro Route Speed - ", route),
         device = pdf,
         height = 8.5,
         width = 11,
         units = "in",
         create.dir = TRUE)
}


# ggmap(basemap) +
#   labs(title = paste0("Metro Route Speed - ", route),
#        subtitle = paste0("averaged between ",
#                          sum(route_counts %>% filter(pid %in% route_focus) %>% pull(route_count)),
#                          " bus trips - ",
#                          min(date(metro_data$time)),
#                          " to ",
#                          max(date(metro_data$time))),
#        x = NULL,
#        y = NULL) +
#   theme(axis.text=element_blank(),
#         axis.ticks=element_blank(),
#         plot.caption = element_text(color = "grey")) +
#   geom_sf(data = metro_segments %>% filter(name %in% route),
#           inherit.aes = FALSE,
#           aes(color = spd_calc),
#           linewidth = 1) +
#   scale_color_distiller(palette = "RdYlGn", direction = "reverse", limits = c(0,70), name = "Average speed or segment\n(calculated with locations, not reported speed)")