library(tidyverse)
library(ggmap)
library(sf)
library(osrm)
library(smoothr)
library(ggnewscale)
library(RColorBrewer)
library(magick)
library(rsvg)
library(parallel)

## add data from WiscTransPortal Crash Data Retrieval Facility ----
## query:  SELECT * 
## FROM DTCRPRD.SUMMARY_COMBINED C 
## WHERE C.CRSHDATE BETWEEN TO_DATE('2022-JAN','YYYY-MM') AND 
## LAST_DAY(TO_DATE('2022-DEC','YYYY-MM')) AND 
## (C.BIKEFLAG = 'Y' OR C.PEDFLAG = 'Y') 
## ORDER BY C.DOCTNMBR

## Load TOPS data ----
## load TOPS data for the whole state (crashes involving bikes and pedestrians),
TOPS_data <- as.list(NULL)
for (file in list.files(path = "data/TOPS", pattern = "crash-data-download")) {
  message(paste("importing data from file: ", file))
  year <- substr(file, 21, 24)
  csv_run <- read_csv(file = paste0("data/TOPS/",file), col_types = cols(.default = "c"))
  TOPS_data[[file]] <- csv_run
}
rm(csv_run, file, year)
TOPS_data <- bind_rows(TOPS_data)

## clean up data
TOPS_data <- TOPS_data %>% 
  mutate(date = mdy(CRSHDATE), 
         age1 = as.double(AGE1),
         age2 = as.double(AGE2),
         latitude  = as.double(LATDECDG),
         longitude = as.double(LONDECDG)) %>% 
  mutate(month = month(date, label = TRUE),
         year = as.factor(year(date)))

# Injury Severy Index and Color -----
injury_severity <- data.frame(InjSevName = c("No apparent injury", "Possible Injury", "Suspected Minor Injury","Suspected Serious Injury","Fatality"), 
                              code = c("O", "C", "B", "A", "K"),
                              #                 color = c("#fafa6e", "#edc346", "#d88d2d", "#bd5721", "#9b1c1c"))
                              color = c("#fafa6e", "#edc346", "#d88d2d", "#d88d21", "#9b1c1c" ))

TOPS_data <- left_join(TOPS_data, injury_severity %>% select(InjSevName, code), join_by(INJSVR1 == code)) %>% 
  mutate(InjSevName = factor(InjSevName, levels = injury_severity$InjSevName)) %>%
  rename(InjSevName1 = InjSevName)
TOPS_data <- left_join(TOPS_data, injury_severity %>% select(InjSevName, code), join_by(INJSVR2 == code)) %>% 
  mutate(InjSevName = factor(InjSevName, levels = injury_severity$InjSevName)) %>%
  rename(InjSevName2 = InjSevName)

TOPS_data <- TOPS_data %>% mutate(ped_inj = ifelse(ROLE1 %in% c("BIKE", "PED"), 
                                      INJSVR1, 
                                      ifelse(ROLE2 %in% c("BIKE", "PED"), 
                                             INJSVR2,
                                             NA)))

TOPS_data <- left_join(TOPS_data, injury_severity %>% select(InjSevName, code), join_by(ped_inj == code)) %>% 
  mutate(InjSevName = factor(InjSevName, levels = injury_severity$InjSevName)) %>%
  rename(ped_inj_name = InjSevName)

# Race names
race <- data.frame(race_name = c("Asian", "Black", "Indian","Hispanic","White"), 
                   code = c("A", "B", "I", "H", "W"))

TOPS_data <- left_join(TOPS_data, race %>% select(race_name, code), join_by(RACE1 == code)) %>% rename(race_name1 = race_name)
TOPS_data <- left_join(TOPS_data, race %>% select(race_name, code), join_by(RACE2 == code)) %>% rename(race_name2 = race_name)

logo <- image_read(path = "other/BFW_Logo_180_x_200_transparent_background.png")


## set tile server info
# register stadia API key ----
register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, 36))
options(ggmap.file_drawer = "basemaps")
# dir.create(file_drawer(), recursive = TRUE, showWarnings = FALSE)
# saveRDS(list(), file_drawer("index.rds"))
readRDS(file_drawer("index.rds"))
file_drawer("index.rds")

## set parameters -----

focus_muni <- c("MILWAUKEE", "MADISON")
focus_inj <- c("A", "K")
focus_role <- c("BIKE", "PED")
focus_years <- c("2023")

## generate maps for focus city

for(muni in focus_muni) {
  # create bounding box around crashes that happen in city.
  muni_data <- TOPS_data %>% filter(MUNINAME %in% muni)
  bbox <- c(left = min(muni_data$longitude, na.rm = TRUE), 
            bottom = min(muni_data$latitude, na.rm = TRUE), 
            right = max(muni_data$longitude, na.rm = TRUE), 
            top = max(muni_data$latitude, na.rm = TRUE))
  
  #get basemap
  basemap <- get_stadiamap(bbox = bbox, zoom = 12, maptype = "stamen_toner_lite")
  
  # generate map
  ggmap(basemap) +
    labs(title = paste0("Crashes between pedestrians/bicyclists in ", str_to_title(muni)),
         subtitle = paste0("that result in a severe injury or fatality | ",
                           focus_years),
         caption = "data from Wisconsin DOT, UW TOPS Laboratory, and OpenStreetMap",
         x = NULL,
         y = NULL) +
    theme(axis.text=element_blank(),
          axis.ticks=element_blank()) +
    
    ## add bike lts
    #geom_sf(data = bike_lts[[county]],
    #       inherit.aes = FALSE,
    #       aes(color = lts)) +
    #scale_color_manual(values = bike_lts_scale$color, name = "Bike Level of Traffic Stress") +
    
    # add crash locations
    new_scale_fill() +
    geom_point(data = TOPS_data %>%
                 filter(ROLE1 %in% focus_role 
                        & INJSVR1 %in% focus_inj
                        #   & age1 < 18 
                        | ROLE2 %in% focus_role
                        & INJSVR2 %in% focus_inj
                        #  & age2 < 18
                 ) %>%
                 filter(longitude >= as.double(bbox[1]),
                        latitude >= as.double(bbox[2]),
                        longitude <= as.double(bbox[3]),
                        latitude <= as.double(bbox[4])) %>%
                 filter(year %in% focus_years),
               aes(x = longitude,
                   y = latitude,
                   fill = ped_inj_name),
               shape = 21,
               size = 2) +
    scale_fill_manual(values = injury_severity %>% filter(code %in% focus_inj) %>% pull(color), name = "Crash Severity") +
    annotation_raster(logo,
                      # Position adjustments here using plot_box$max/min/range
                      ymin = bbox['top'] - 0.25 * 0.16,
                      ymax = bbox['top'],
                      xmin = bbox['right'] + 0.25 * 0.05,
                      xmax = bbox['right'] + 0.25 * 0.20) +
    coord_sf(clip = "off")
  
  ggsave(file = paste0("figures/city_maps/",
                       str_to_title(muni),
                       ".pdf"),
         
         title = paste0(str_to_title(muni), " Pedestrian/Bike crashes"),
         device = pdf,
         height = 8.5,
         width = 11,
         units = "in",
         create.dir = TRUE)
}