---
title: "MilWALKeeWalks"
output: 
  html_document:
    toc: true
    toc_depth: 5
    toc_float:
      collapsed: false
      smooth_scroll: true
editor_options:
  chunk_output_type: console
---


# Input Data & Configuration

## Libraries
```{r libs, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
date()
rm(list=ls())
library(tidyverse)
library(ggmap)
library(sf)
library(osrm)
library(smoothr)
library(ggnewscale)
library(RColorBrewer)
library(magick)
library(rsvg)
library(tidycensus)
```


## Load TOPS data
```{r loadTOPS, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
load(file = "data/TOPS/TOPS_data.Rda")
load(file = "data/TOPS/vuln_roles.Rda")
load(file = "data/TOPS/retrieve_date.Rda")
load(file = "data/TOPS/injury_severity.Rda")
```

## filter to county
```{r filterTOPS, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
focus_county <- "MILWAUKEE"
TOPS_data_filtered <- TOPS_data %>% filter(CNTYNAME == focus_county)
```

## identify start and end dates
```{r startenddates, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
year_min <- min(year(TOPS_data_filtered$date))
year_max <- max(year(TOPS_data_filtered$date))
month_end <- 9
date_max <- ymd(paste(year_max, month_end + 1,1)) - 1
```

## intro charts
```{r introCharts, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
ggplot() +
  geom_col(data = TOPS_data_filtered %>%
#             filter(MUNINAME %in% "MILWAUKEE") %>%
             filter(! year %in% c(year_max)) %>%
             filter(! is.na(vulnerable_role)) %>%
             group_by(month, vulnerable_role) %>% 
             summarize(total = n()),
           aes(x = month,
               y = total/((year_max - 1) - year_min + 1),
               fill = vulnerable_role),
           position = position_dodge()) +
  geom_line(data = TOPS_data_filtered %>%
#             filter(MUNINAME %in% "MILWAUKEE") %>%
              filter(year == year_max) %>%
              filter(!is.na(vulnerable_role)) %>%
              group_by(month, vulnerable_role) %>% 
              summarize(total = n()),
            aes(x = month,
                y = total,
                color = vulnerable_role,
                group = vulnerable_role),
            linewidth = 1) +
  scale_y_continuous(expand = expansion(mult = c(0,0.1))) +
  scale_fill_manual(values = c("sienna3", "deepskyblue3")) +
  scale_color_manual(values = c("sienna4", "deepskyblue4")) +
  labs(title = paste0("Car crashes involving pedestrians & bicyclists"),
       subtitle = paste0(str_to_title(focus_county), " County"),
       x = "Month",
       y = "Crashes per month",
       fill = paste0("Average\n", year_min, " - ", year_max - 1),
       color = year_max,
       caption = paste0("crash data from UW TOPS lab - retrieved ",
                        strftime(retrieve_date, format = "%m/%Y"),
                        "\nper direction of the WisDOT Bureau of Transportation Safety")) +
      theme(plot.caption = element_text(color = "grey"))
ggsave(filename = paste0("figures/MilWALKee_Walks/", "month_role.png"),
       device = png,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

ggplot() +
  geom_col(data = TOPS_data_filtered %>%
#             filter(MUNINAME %in% "MILWAUKEE") %>%
             filter(vulnerable_role == "Pedestrian",
                    !is.na(ped_age)) %>%
             filter(! year %in% c(year_max)) %>%
             mutate(age = ifelse(ped_age <= 18, "child", "adult")) %>%
             group_by(month, age) %>% 
             summarize(total = n()/((year_max - 1) - year_min + 1)),
           aes(x = month,
               y = total,
               fill = age),
           position = position_dodge()) +
    geom_line(data = TOPS_data_filtered %>%
#                filter(MUNINAME %in% "MILWAUKEE") %>%
                filter(year == year_max) %>%
                filter(vulnerable_role == "Pedestrian",
                       !is.na(ped_age)) %>%
                mutate(age = ifelse(ped_age <= 18, "child", "adult")) %>%
                group_by(month, age, year) %>% 
                summarize(total = n()),
              aes(x = month,
                  y = total,
                  color = age,
                  group = age),
              linewidth = 1) +
  scale_y_continuous(expand = expansion(mult = c(0,0.1))) +
  scale_fill_manual(values = c("deeppink1", "darkgoldenrod1")) +
  scale_color_manual(values = c("deeppink3", "darkgoldenrod3")) +
  labs(title = paste0("Car crashes involving pedestrians"),
       subtitle = paste0(str_to_title(focus_county), " County"),
       x = "Month",
       y = "Crashes per month",
       fill = paste0("Average\n", year_min, " - ", year_max - 1),
       color = year_max,
       caption = paste0("crash data from UW TOPS lab - retrieved ",
                        strftime(retrieve_date, format = "%m/%Y"),
                        "\nper direction of the WisDOT Bureau of Transportation Safety")) +
      theme(plot.caption = element_text(color = "grey"))
ggsave(filename = paste0("figures/MilWALKee_Walks/", "month_age.png"),
       device = png,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

ggplot(data = TOPS_data_filtered %>%
  filter(
    #vulnerable_role == "Pedestrian",
         month(date) <= month_end,
         ) %>%
#    filter(ped_age < 18,
#          MUNINAME %in% "MILWAUKEE") %>%
  group_by(year) %>% 
  summarize(total = n())) +
  geom_col(aes(x = year,
               y = total),
           fill = "lightblue4") +
  scale_y_continuous(expand = expansion(mult = c(0,0.1))) +
  labs(title = paste0("Car crashes involving pedestrians & bicyclists"),
       subtitle = paste0(str_to_title(focus_county), " County | ", "January - ", format(ymd(paste(2025, month_end, 01)), "%B")),
       x = NULL,
       y = "Crashes per year",
       caption = paste0("crash data from UW TOPS lab - retrieved ",
                        strftime(retrieve_date, format = "%m/%Y"),
                        "\nper direction of the WisDOT Bureau of Transportation Safety")) +
      theme(plot.caption = element_text(color = "grey"))
ggsave(filename = paste0("figures/MilWALKee_Walks/", "vuln_years.png"),
       device = png,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

ggplot(data = TOPS_data_filtered %>%
  filter(
    #vulnerable_role == "Pedestrian",
         month(date) <= month_end,
         ped_inj %in% c("K", "A"),
         ) %>%
  group_by(year, ped_inj_name) %>% 
  summarize(total = n())) +
  geom_col(aes(x = year,
               y = total,
               fill = ped_inj_name),
           position = position_dodge()) +
  scale_y_continuous(expand = expansion(mult = c(0,0.1))) +
  scale_fill_manual(values = setNames(injury_severity$color, injury_severity$InjSevName), name = "Injury severity") +
labs(title = paste0("Car crashes involving pedestrians & bicyclists - fatal and serious injuries"),
     subtitle = paste0(str_to_title(focus_county), " County | ", "January - ", format(ymd(paste(2025, month_end, 01)), "%B")),
     x = NULL,
     y = "Crashes per year",
     caption = paste0("crash data from UW TOPS lab - retrieved ",
                      strftime(retrieve_date, format = "%m/%Y"),
                      "\nper direction of the WisDOT Bureau of Transportation Safety")) +
      theme(plot.caption = element_text(color = "grey"))
ggsave(filename = paste0("figures/MilWALKee_Walks/", "vuln_years_severe.png"),
       device = png,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

ggplot(data = TOPS_data_filtered %>%
         filter(
           #vulnerable_role == "Pedestrian",
           ped_inj %in% c("K", "A"),
                month(date) <= month_end) %>%
         mutate(mke_city = factor(ifelse(MUNINAME %in% "MILWAUKEE", "City of Milwaukee", "Other municipalities"),
                                  levels = c("Other municipalities", "City of Milwaukee"))) %>%
         group_by(year, mke_city) %>% 
         summarize(total = n())) +
  geom_col(aes(x = year,
               y = total,
               fill = mke_city),
           position = position_dodge()) +
  scale_y_continuous(expand = expansion(mult = c(0,0.1))) +
  labs(title = paste0("Car crashes involving pedestrians - fatal and severe injuries"),
       subtitle = paste0(str_to_title(focus_county), " County | ", "January - ", format(ymd(paste(2025, month_end, 01)), "%B")),
       x = NULL,
       y = "Crashes",
       fill = NULL,
       caption = paste0("crash data from UW TOPS lab - retrieved ",
                        strftime(retrieve_date, format = "%m/%Y"),
                        "\nper direction of the WisDOT Bureau of Transportation Safety")) +
      theme(plot.caption = element_text(color = "grey"))
ggsave(filename = paste0("figures/MilWALKee_Walks/", "ped_years_MKEcity.png"),
       device = png,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

ggplot(data = TOPS_data_filtered %>%
  filter(
    #vulnerable_role == "Pedestrian",
      #   month(date) <= 8,
         ped_inj %in% c("K", "A"),
         ) %>%
  group_by(year, month) %>% 
  summarize(total = n()) %>%
    arrange(year, month) %>%
    group_by(year) %>%
    mutate(cumcrashes = cumsum(total))) +
  geom_line(aes(x = month,
               y = cumcrashes,
               group = year,
               color = year)) +
  scale_y_continuous(expand = expansion(mult = c(0,0.1))) +
  scale_color_brewer(palette = "Set1") +
  scale_fill_manual(values = setNames(injury_severity$color, injury_severity$InjSevName), name = "Injury severity") +
labs(title = paste0("Car crashes involving pedestrians & bicyclists - fatal and serious injuries"),
     subtitle = paste0(str_to_title(focus_county), " County"),
     x = NULL,
     y = "Cumulative crashes",
     color = NULL,
     caption = paste0("crash data from UW TOPS lab - retrieved ",
                      strftime(retrieve_date, format = "%m/%Y"),
                      "\nper direction of the WisDOT Bureau of Transportation Safety")) +
      theme(plot.caption = element_text(color = "grey"))
ggsave(filename = paste0("figures/MilWALKee_Walks/", "vuln_years_cumulative_severe.png"),
       device = png,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

```

## Milwaukee maps

## Load API keys from StadiaMaps
```{r APIkeys, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
# register stadia API key ----
register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, 36))
```

## add county census data ----
```{r countycensus, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
census_api_key(key = substr(read_file(file = "api_keys/census_api_key"), 1, 40))
county_populations <- get_estimates(geography = "county", 
                                    year = 2022, 
                                    product = "population", 
                                    state = "Wisconsin", 
                                    geometry = TRUE) %>%
  filter(variable == "POPESTIMATE") %>%
  mutate(County = str_to_upper(str_replace(NAME, " County, Wisconsin", "")))
county_populations <- st_transform(county_populations, crs = 4326) %>% filter(County %in% focus_county)
```

## make hex grids and join with crash data ----
```{r hexgrids, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}

hexgrid <- rowid_to_column(st_transform(st_as_sf(st_make_grid(st_transform(county_populations, crs = 32616),
               cellsize = 3000,
               what = 'polygons',
               square = FALSE
               )), crs = 4326), "ID")
yearsforprior <- year_max - year_min - 1
hex_crashes <- st_join(hexgrid, 
                       st_as_sf(TOPS_data_filtered %>% 
                                  filter(vulnerable_role %in% "Pedestrian") %>% 
                                  #filter(ped_inj %in% c("A", "K", "B")) %>%
                                  filter(!is.na(latitude)), 
                                coords = c("longitude", "latitude"), 
                                crs = 4326), 
                       join = st_contains) %>%
  filter(!is.na(year)) %>%
  filter(date > (date_max - (365 * (yearsforprior + 1)))) %>%
  mutate(lastyear = ifelse((date <= date_max - 365),
                "prior",
                "lastyear")) %>%
  group_by(ID, lastyear) %>%
  summarise(count = n(), .groups = 'drop') %>%
  st_drop_geometry() %>%
  pivot_wider(id_cols = ID, names_from = lastyear, values_from = count) %>%
  mutate(across(-ID, ~ replace_na(., 0))) %>%
  mutate(total = rowSums(dplyr::select(., -ID), na.rm = TRUE))

hex_crashes <- st_as_sf(left_join(hexgrid, hex_crashes), crs = 4326)

hex_crashes <- hex_crashes %>%
  mutate(prioryearlyaverage = prior/yearsforprior) %>%
  mutate(lastyearchangepercent = (lastyear - prioryearlyaverage)/prioryearlyaverage,
         lastyearchangecrashes = lastyear - prioryearlyaverage)

hex_crashes_points <- st_centroid(hex_crashes)

#hex_small_crashes

hexgrid_small <- rowid_to_column(st_transform(st_as_sf(st_make_grid(st_transform(county_populations, crs = 32616),
               cellsize = 500,
               what = 'polygons',
               square = FALSE
               )), crs = 4326), "ID")

hex_crashes_small <- st_join(hexgrid_small, 
                             st_as_sf(TOPS_data_filtered %>% 
                                        filter(vulnerable_role %in% "Pedestrian") %>%
                                        #filter(ped_inj %in% c("A", "K", "B")) %>%
                                        filter(!is.na(latitude)), 
                                      coords = c("longitude", "latitude"), 
                                      crs = 4326), 
                             join = st_contains) %>%
  filter(!is.na(year)) %>%
  filter(date > (date_max - (365 * (yearsforprior + 1)))) %>%
  mutate(lastyear = ifelse((date <= date_max - 365),
                "prior",
                "lastyear")) %>%
  group_by(ID, lastyear) %>%
  summarise(count = n(), .groups = 'drop') %>%
  st_drop_geometry() %>%
  pivot_wider(id_cols = ID, names_from = lastyear, values_from = count) %>%
  mutate(across(-ID, ~ replace_na(., 0))) %>%
  mutate(total = rowSums(dplyr::select(., -ID), na.rm = TRUE))

hex_crashes_small <- st_as_sf(left_join(hexgrid_small, hex_crashes_small), crs = 4326)

hex_crashes_small <- hex_crashes_small %>%
  mutate(prioryearlyaverage = prior/yearsforprior) %>%
  mutate(lastyearchangepercent = (lastyear - prioryearlyaverage)/prioryearlyaverage,
         lastyearchangecrashes = lastyear - prioryearlyaverage)

hex_crashes_small_points <- st_centroid(hex_crashes_small)


```

```{r previousyearstring, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}

previousyearstring <- paste0(format(date_max - 365 + 1, "%m/%Y"), " - ", format(date_max, "%m/%Y"))
```

```{r importprojects, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}

projects_2024 <- read_sf(dsn = "data/MilWALKeeWalks/projects_2024.kml")
projects_2023 <- read_sf(dsn = "data/MilWALKeeWalks/projects_2023.kml")
```

```{r MilwaukeeMaps, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
# get basemap
bbox <- st_bbox(county_populations)
bbox <- c(left = as.double(bbox[1]), 
          bottom = as.double(bbox[2]), 
          right = as.double(bbox[3]), 
          top = as.double(bbox[4]))
basemap <- get_stadiamap(bbox = bbox, zoom = 12, maptype = "stamen_toner_lite")

# generate map with bubbles
ggmap(basemap) +
  labs(title = paste0("Car crashes involving pedestrians"),
    subtitle = paste0(str_to_title(focus_county),
                      " County | ",
                      year_min,
                      " - ",
                      year_max),
    caption = paste0("crash data from UW TOPS lab - retrieved ",
                     strftime(retrieve_date, format = "%m/%Y"),
                     "\nper direction of the WisDOT Bureau of Transportation Safety",
                     "\nbasemap from StadiaMaps and OpenStreetMap Contributers"),
    x = NULL,
    y = NULL,
    size = paste0("Total crashes"),
    fill = paste0(previousyearstring, "\ncompared to\nprevious years")) +
  theme(axis.text=element_blank(),
        axis.ticks=element_blank(),
        plot.caption = element_text(color = "grey", size = 8)) +
  # add crash locations
  geom_sf(data = hex_crashes_points %>% filter(is.double(total), !is.na(total)),
          inherit.aes = FALSE,
          aes(size = total,
              fill = lastyearchangepercent),
          linewidth = 0,
          shape = 21,
          color = "black") +
  scale_size_area() +
  scale_fill_gradient2(
    low = "darkgreen",
    mid = "white",
    high = "red",
    midpoint = 0,
    limits = c(-2, 2),
    oob = scales::squish,
    labels = scales::percent
  )

ggsave(file = paste0("figures/MilWALKee_Walks/",
                     "milwaukee_map.png"),
       device = png,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

ggmap(basemap) +
  labs(title = paste0("Car crashes involving pedestrians"),
    subtitle = paste0(str_to_title(focus_county),
                      " County | ",
                      previousyearstring),
    caption = paste0("crash data from UW TOPS lab - retrieved ",
                     strftime(retrieve_date, format = "%m/%Y"),
                     "\nper direction of the WisDOT Bureau of Transportation Safety",
                     "\nbasemap from StadiaMaps and OpenStreetMap Contributers"),
    x = NULL,
    y = NULL,
    size = paste0("Total crashes"),
    fill = "Crashes in the last year\ncompared to previous average") +
  theme(axis.text=element_blank(),
        axis.ticks=element_blank(),
        plot.caption = element_text(color = "grey", size = 8)) +
  # add crash locations
  geom_sf(data = hex_crashes,
          inherit.aes = FALSE,
          aes(fill = lastyearchangecrashes),
          alpha = 0.5) +
  scale_size_area() +
  scale_fill_gradient2(
    low = "darkgreen",
    mid = "white",
    high = "red",
    midpoint = 0,
#    limits = c(-20, 20),
#    oob = scales::squish
  ) # + geom_sf(data = projects_2023, inherit.aes = FALSE)

ggsave(file = paste0("figures/MilWALKee_Walks/",
                     "milwaukee_map_crashchange.png"),
       device = png,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

highlighted_areas <- hex_crashes %>% 
  mutate(lastyearchange = (lastyear - prior/yearsforprior)) %>% 
  filter(prior > 100,
         (lastyearchange > 5 |
            lastyearchange < -2)) %>% pull(ID)
highlighted_areas <- c(62, 69, 78, 85)

ggmap(basemap) +
  labs(title = paste0("Car crashes involving pedestrians"),
    subtitle = paste0(str_to_title(focus_county),
                      " County | ",
                      year_min, 
                      " - ", 
                      year_max),
    caption = paste0("crash data from UW TOPS lab - retrieved ",
                     strftime(retrieve_date, format = "%m/%Y"),
                     "\nper direction of the WisDOT Bureau of Transportation Safety",
                     "\nbasemap from StadiaMaps and OpenStreetMap Contributers"),
    x = NULL,
    y = NULL,
    size = paste0("Total crashes"),
    fill = "last year\ncompared to previous") +
  theme(axis.text=element_blank(),
        axis.ticks=element_blank(),
        plot.caption = element_text(color = "grey", size = 8)) +
  # add crash locations
  geom_sf(data = hex_crashes %>% filter(ID %in% highlighted_areas),
          inherit.aes = FALSE,
          aes(fill = lastyearchangepercent),
          alpha = 0.5) +
  geom_sf_text(data = hex_crashes_points %>% filter(ID %in% highlighted_areas),
                inherit.aes = FALSE,
                aes(label = paste0(ifelse(lastyearchangepercent > 0, "+",""),round(lastyearchangepercent * 100, 0), "%")),
               size = 3) +
  scale_size_area() +
  scale_fill_gradient2(
    low = "darkgreen",
    mid = "white",
    high = "red",
    midpoint = 0,
    limits = c(-1, 1),
    oob = scales::squish,
    labels = scales::percent
  )

ggsave(file = paste0("figures/MilWALKee_Walks/",
                     "milwaukee_map_highlighted.png"),
       device = png,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

bbox <- st_bbox(st_union(hexgrid %>% filter(ID %in% highlighted_areas)))
bbox <- c(left = as.double(bbox[1]), 
          bottom = as.double(bbox[2]), 
          right = as.double(bbox[3]), 
          top = as.double(bbox[4]))
basemap <- get_stadiamap(bbox = bbox, zoom = 14, maptype = "stamen_toner_lite")

# Map of high increase areas
ggmap(basemap) +
  labs(title = paste0("Car crashes involving pedestrians"),
    subtitle = paste0(str_to_title(focus_county),
                      " County | ",
                      year_min, 
                      " - ", 
                      year_max),
    caption = paste0("crash data from UW TOPS lab - retrieved ",
                     strftime(retrieve_date, format = "%m/%Y"),
                     "\nper direction of the WisDOT Bureau of Transportation Safety",
                     "\nbasemap from StadiaMaps and OpenStreetMap Contributers"),
    x = NULL,
    y = NULL,
    size = paste0("Total crashes"),
    fill = paste0("last 12 months\ncompared to previous\n", yearsforprior, " years")) +
  theme(axis.text=element_blank(),
        axis.ticks=element_blank(),
        plot.caption = element_text(color = "grey", size = 8)) +
#add hexagons
  new_scale_fill() +
  geom_sf(data = hex_crashes_small %>%
            filter(!is.na(lastyearchangecrashes)),
          inherit.aes = FALSE,
          aes(fill = lastyearchangecrashes),
          alpha = 0.5) +
  scale_size_area() +
  scale_fill_gradient2(
    low = "darkgreen",
    mid = "white",
    high = "red",
    midpoint = 0,
#    limits = c(-10, 10),
#    oob = scales::squish,
  name = "Crashes last year\ncompared to previous average") +
  # add crash locations
  new_scale_fill() +
  geom_point(data = TOPS_data_filtered %>% 
               filter(date > date_max - 365) %>%
               #filter(ped_inj %in% c("K", "A", "B")) %>%
               filter(vulnerable_role %in% "Pedestrian") %>%
               filter(longitude >= as.double(bbox[1]),
                      latitude >= as.double(bbox[2]),
                      longitude <= as.double(bbox[3]),
                      latitude <= as.double(bbox[4])) %>%
               arrange(ped_inj_name),
             aes(x = longitude,
                 y = latitude,
                 fill = ped_inj_name),
             shape = 23,
             size = 3) +
  scale_fill_manual(values = setNames(injury_severity$color, injury_severity$InjSevName), name = paste0("Crashes ", previousyearstring))# + geom_sf(data = projects_2023, inherit.aes = FALSE)

ggsave(file = paste0("figures/MilWALKee_Walks/",
                     "milwaukee_map_zoomchange.png"),
       device = png,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

## compare crashes in area
nrow(TOPS_data_filtered %>% 
               filter(date > date_max - 365) %>%
               filter(ped_inj %in% c("K", "A")) %>%
               filter(vulnerable_role %in% "Pedestrian") %>%
               filter(longitude >= as.double(bbox[1]),
                      latitude >= as.double(bbox[2]),
                      longitude <= as.double(bbox[3]),
                      latitude <= as.double(bbox[4])))
nrow(TOPS_data_filtered %>% 
               filter(date > (date_max - 365 * (yearsforprior + 1))) %>%
               filter(ped_inj %in% c("K", "A")) %>%
               filter(vulnerable_role %in% "Pedestrian") %>%
               filter(longitude >= as.double(bbox[1]),
                      latitude >= as.double(bbox[2]),
                      longitude <= as.double(bbox[3]),
                      latitude <= as.double(bbox[4])))/(yearsforprior + 1)
##highland ave

bbox <- c(left = -87.967, 
          bottom = 43.043, 
          right = -87.944, 
          top = 43.051)
basemap <- get_stadiamap(bbox = bbox, zoom = 15, maptype = "stamen_toner_lite")
ggmap(basemap) +
  labs(title = paste0("Crashes between cars and pedestrians"),
    subtitle = paste0(str_to_title(focus_county),
                      " County | ",
                      year_min, 
                      " - ", 
                      year_max),
    caption = paste0("crash data from UW TOPS lab - retrieved ",
                     strftime(retrieve_date, format = "%m/%Y"),
                     "\nper direction of the WisDOT Bureau of Transportation Safety",
                     "\nbasemap from StadiaMaps and OpenStreetMap Contributers"),
    x = NULL,
    y = NULL,
    size = paste0("Total crashes"),
    fill = paste0("last 12 months\ncompared to previous\n", yearsforprior, " years")) +
  theme(axis.text=element_blank(),
        axis.ticks=element_blank(),
        plot.caption = element_text(color = "grey", size = 8)) +
  # add crash locations
  new_scale_fill() +
  geom_point(data = TOPS_data_filtered %>% 
               filter(year == 2025) %>%
               #filter(ped_inj %in% c("K", "A", "B")) %>%
               filter(vulnerable_role %in% "Pedestrian") %>%
               filter(longitude >= as.double(bbox[1]),
                      latitude >= as.double(bbox[2]),
                      longitude <= as.double(bbox[3]),
                      latitude <= as.double(bbox[4])) %>%
               arrange(ped_inj_name),
             aes(x = longitude,
                 y = latitude,
                 fill = ped_inj_name),
             shape = 23,
             size = 3) +
  scale_fill_manual(values = setNames(injury_severity$color, injury_severity$InjSevName), name = paste0("Crashes ", previousyearstring))

bbox <- c(left = -87.910, 
          bottom = 43.032, 
          right = -87.896, 
          top = 43.044)
basemap <- get_stadiamap(bbox = bbox, zoom = 15, maptype = "stamen_toner_lite")
ggmap(basemap) +
  labs(title = paste0("Crashes between cars and pedestrians"),
    subtitle = paste0(str_to_title(focus_county),
                      " County | ",
                      year_min, 
                      " - ", 
                      year_max),
    caption = paste0("crash data from UW TOPS lab - retrieved ",
                     strftime(retrieve_date, format = "%m/%Y"),
                     "\nper direction of the WisDOT Bureau of Transportation Safety",
                     "\nbasemap from StadiaMaps and OpenStreetMap Contributers"),
    x = NULL,
    y = NULL,
    size = paste0("Total crashes"),
    fill = paste0("last 12 months\ncompared to previous\n", yearsforprior, " years")) +
  theme(axis.text=element_blank(),
        axis.ticks=element_blank(),
        plot.caption = element_text(color = "grey", size = 8)) +
  # add crash locations
  new_scale_fill() +
  geom_point(data = TOPS_data_filtered %>% 
               filter(year != 2025) %>%
               #filter(ped_inj %in% c("K", "A", "B")) %>%
               filter(vulnerable_role %in% "Pedestrian") %>%
               filter(longitude >= as.double(bbox[1]),
                      latitude >= as.double(bbox[2]),
                      longitude <= as.double(bbox[3]),
                      latitude <= as.double(bbox[4])) %>%
               arrange(ped_inj_name),
             aes(x = longitude,
                 y = latitude,
                 fill = ped_inj_name),
             shape = 23,
             size = 3) +
  scale_fill_manual(values = setNames(injury_severity$color, injury_severity$InjSevName), name = paste0("Crashes ", previousyearstring))
```


## identify Halloween trick-or-treating days
```{r trickortreatdays, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
trickortreatdays <- data_frame(year = seq(year_min, year_max, 1))
trickortreatdays <- trickortreatdays %>%
  mutate(halloween = ymd(paste(year, "10, 31"))) %>%
  mutate(wday = wday(halloween, label = TRUE)) %>%
  mutate(satbefore = floor_date(halloween, "week", week_start = 6),
         sunbefore = floor_date(halloween, "week"))

trickortreatdays <- c(trickortreatdays$halloween, trickortreatdays$satbefore, trickortreatdays$sunbefore)
TOPS_data_filtered <- TOPS_data_filtered %>% mutate(trickortreat = ifelse(date %in% trickortreatdays, TRUE, FALSE))
TOPS_data <- TOPS_data %>% mutate(trickortreat = ifelse(date %in% trickortreatdays, TRUE, FALSE))

```

## Halloween
```{r Halloween, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
ggplot(data = TOPS_data_filtered %>% 
#         filter(ped_inj %in% c("K", "A", "B")) %>%
#  filter(ped_age <=18) %>%
  filter(vulnerable_role == "Pedestrian") %>%
  mutate(mday = mday(date)) %>%
  mutate(date_yearagnostic = ymd(paste("2025", month, mday))) %>%
  group_by(date_yearagnostic, year, trickortreat) %>% 
  summarize(total = n())) +
  geom_col(aes(x = date_yearagnostic,
               y = total,
               fill = trickortreat)) +
  scale_x_date(minor_breaks = "month", date_labels = "%b", expand = expansion(mult = c(0,0))) +
  scale_y_continuous(breaks = seq(0, 9, 3), labels = c(0, 3, 6, "")) +
  scale_fill_manual(values = c("black", "orange")) +
  facet_grid(year ~ .) +
  labs(title = paste0("Car crashes involving pedestrians - Halloween"),
       subtitle = paste0(str_to_title(focus_county), " County | ", year_min, " - ", year_max),
       x = NULL,
       y = "Crashes per day",
       fill = "Days with\nTrick-or-Treating",
       caption = paste0("crash data from UW TOPS lab - retrieved ",
                        strftime(retrieve_date, format = "%m/%Y"),
                        "\nper direction of the WisDOT Bureau of Transportation Safety")) +
      theme(plot.caption = element_text(color = "grey"))
ggsave(filename = paste0("figures/MilWALKee_Walks/", "halloween.png"),
       device = png,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

ggplot(data = TOPS_data_filtered %>% 
#         filter(ped_inj %in% c("K", "A", "B")) %>%
#  filter(ped_age <=18) %>%
  mutate(wday = wday(date, label = TRUE)) %>%
  filter(wday(date) %in% c(1, 7)) %>%
  filter(vulnerable_role == "Pedestrian") %>%
  group_by(date, wday, trickortreat) %>% 
  summarize(total = n()) %>%
  group_by(wday, trickortreat) %>%
  summarize(average = mean(total, na.rm = TRUE))) +
  geom_col(aes(x = wday,
               y = average,
               fill = trickortreat),
           position = position_dodge()) +
  scale_fill_manual(values = c("black", "orange")) +
  labs(title = paste0("Car crashes involving pedestrians - Halloween"),
       subtitle = paste0(str_to_title(focus_county), " County | ", year_min, " - ", year_max),
       x = NULL,
       y = "Crashes per day",
       fill = "Days with\nTrick-or-Treating",
       caption = paste0("crash data from UW TOPS lab - retrieved ",
                        strftime(retrieve_date, format = "%m/%Y"),
                        "\nper direction of the WisDOT Bureau of Transportation Safety")) +
      theme(plot.caption = element_text(color = "grey"))

ggsave(filename = paste0("figures/MilWALKee_Walks/", "halloween_wday.png"),
       device = png,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

```

## Trucks
```{r Trucks, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
TOPS_data_filtered %>% 
  filter(ped_age < 18) %>%
  filter(nonvuln_vehicletype != "Other") %>%
  filter(vulnerable_role == "Pedestrian") %>%
  group_by(ped_inj_name, nonvuln_vehicletype) %>%
  summarise(total = n()) %>%
  group_by(nonvuln_vehicletype) %>%
  mutate(percent = total/sum(total),
         cumpercent = cumsum(total)/sum(total)) %>%
  ggplot(aes(x = nonvuln_vehicletype,
               y = percent,
               fill = ped_inj_name)) +
  geom_col() +
  geom_label(aes(label = paste0(round(percent * 100, 1), "%"),
                fill = ped_inj_name),
                position = position_stack(vjust = 0.5)) +
  scale_fill_manual(values = setNames(injury_severity$color, injury_severity$InjSevName), name = "Injury severity") +
  scale_y_continuous(labels = scales::percent) +
  labs(title = paste0("Car crashes involving pedestrians - Children"),
       subtitle = paste0(str_to_title(focus_county), " County | ", year_min, " - ", year_max),
       x = NULL,
       y = "Percent of crashes",
       caption = paste0("crash data from UW TOPS lab - retrieved ",
                        strftime(retrieve_date, format = "%m/%Y"),
                        "\nper direction of the WisDOT Bureau of Transportation Safety")) +
      theme(plot.caption = element_text(color = "grey"))
ggsave(filename = paste0("figures/MilWALKee_Walks/", "trucks_injsev.png"),
       device = png,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

```

## seasonality
```{r seasonality, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
TOPS_data_filtered %>% 
#         filter(ped_inj %in% c("K", "A", "B")) %>%
  filter(ped_age >= 18) %>%
  filter(vulnerable_role == "Pedestrian",
         year != year_max,
         MUNINAME == "MILWAUKEE") %>%
  mutate(date_yearagnostic = ymd(paste("2025", month(date, label = FALSE), mday(date)))) %>%
  mutate(week = week(date_yearagnostic),
         weekdate = floor_date(date_yearagnostic, unit = "weeks", week_start = 3)) %>%
  group_by(weekdate) %>% 
  summarize(average = n()/(year_max - year_min + 1)) %>%
ggplot() +
  geom_col(aes(x = weekdate,
               y = average),
           fill = "darkorange3") +
  scale_x_date(minor_breaks = "month", date_labels = "%b", expand = expansion(mult = c(0,0))) +
  scale_y_continuous() +
  labs(title = paste0("Car crashes involving pedestrians - Adults"),
       subtitle = paste0("City of Milwaukee | ", year_min, " - ", year_max - 1),
       x = NULL,
       y = "Average crashes per week",
       fill = "Days with\nTrick-or-Treating",
       caption = paste0("crash data from UW TOPS lab - retrieved ",
                        strftime(retrieve_date, format = "%m/%Y"),
                        "\nper direction of the WisDOT Bureau of Transportation Safety")) +
      theme(plot.caption = element_text(color = "grey"))
ggsave(filename = paste0("figures/MilWALKee_Walks/", "MKEadultweek.png"),
       device = png,
       height = 8.5,
       width = 11,
       units = "in",
       create.dir = TRUE)

hex_crashes_small <- st_join(hexgrid_small, 
                             st_as_sf(TOPS_data_filtered %>% 
                                        filter(vulnerable_role %in% "Pedestrian",
                                               ped_age >= 18,
                                               MUNINAME == "MILWAUKEE") %>%
                                        #filter(ped_inj %in% c("A", "K", "B")) %>%
                                        filter(!is.na(latitude)), 
                                      coords = c("longitude", "latitude"), 
                                      crs = 4326), 
                             join = st_contains) %>%
  filter(!is.na(year)) %>%
  mutate(dateyearagnostic = ymd(paste("2025", month(date, label = FALSE), mday(date)))) %>%
  filter(dateyearagnostic >= ymd("2025-10-01"),
         dateyearagnostic <= ymd("2025-10-14")) %>%
  mutate(octoberfest = ifelse(dateyearagnostic <= ymd("2025-10-7"),
                "TRUE",
                "FALSE")) %>%
  group_by(ID, octoberfest) %>%
  summarise(count = n(), .groups = 'drop') %>%
  st_drop_geometry() %>%
  pivot_wider(id_cols = ID, names_from = octoberfest, values_from = count) %>%
  mutate(across(-ID, ~ replace_na(., 0))) %>%
  mutate(total = rowSums(dplyr::select(., -ID), na.rm = TRUE))

hex_crashes_small <- st_as_sf(left_join(hexgrid_small, hex_crashes_small), crs = 4326)

hex_crashes_small <- hex_crashes_small %>%
  mutate(octoberfestchange = `TRUE` - `FALSE`)
# Map of high increase areas
ggmap(basemap) +
  labs(title = paste0("Car crashes involving pedestrians"),
    subtitle = paste0(str_to_title(focus_county),
                      " County | ",
                      year_min, 
                      " - ", 
                      year_max),
    caption = paste0("crash data from UW TOPS lab - retrieved ",
                     strftime(retrieve_date, format = "%m/%Y"),
                     "\nper direction of the WisDOT Bureau of Transportation Safety",
                     "\nbasemap from StadiaMaps and OpenStreetMap Contributers"),
    x = NULL,
    y = NULL,
    size = paste0("Total crashes"),
    fill = paste0("last 12 months\ncompared to previous\n", yearsforprior, " years")) +
  theme(axis.text=element_blank(),
        axis.ticks=element_blank(),
        plot.caption = element_text(color = "grey", size = 8)) +
#add hexagons
  new_scale_fill() +
  geom_sf(data = hex_crashes_small %>%
            filter(!is.na(octoberfestchange)),
          inherit.aes = FALSE,
          aes(fill = octoberfestchange),
          alpha = 0.5) +
  scale_size_area() +
  scale_fill_gradient2(
    low = "darkgreen",
    mid = "white",
    high = "red",
    midpoint = 0,
#    limits = c(-10, 10),
#    oob = scales::squish,
  name = "Crashes last year\ncompared to previous average")

+
  # add crash locations
  new_scale_fill() +
  geom_point(data = TOPS_data_filtered %>% 
               #filter(date > max(TOPS_data_filtered$date) - 365) %>%
               #filter(ped_inj %in% c("K", "A", "B")) %>%
               filter(date > ymd("2024-09-20"),
                      date < ymd("2024-10-10")) %>%
               filter(vulnerable_role %in% "Pedestrian") %>%
               filter(longitude >= as.double(bbox[1]),
                      latitude >= as.double(bbox[2]),
                      longitude <= as.double(bbox[3]),
                      latitude <= as.double(bbox[4])) %>%
               arrange(ped_inj_name),
             aes(x = longitude,
                 y = latitude,
                 fill = ped_inj_name),
             shape = 23,
             size = 3) +
  scale_fill_manual(values = setNames(injury_severity$color, injury_severity$InjSevName), name = paste0("Crashes ", previousyearstring))# + geom_sf(data = projects_2023, inherit.aes = FALSE)
```