Ped Crash 02

library(sf)
library(dplyr)
library(ggplot2)
library(tigris)
library(readr)
library(RColorBrewer)
library(classInt)
library(leaflet)
library(htmlwidgets)

options(tigris_use_cache = TRUE)


setwd("C:/Users/mvx13/OneDrive - Texas State University/Papers/2026/2027_TRBAM/Data/OtherData")
library(readxl)

dat= read_excel("Pedestrian2017_2025a.xlsx")
dim(dat)

names(dat)

dat1= dat[, c(1, 24, 25, 28, 35)]
data_unique <- dat1[!duplicated(dat1$Crash_ID), ]
dim(data_unique)

crash= data_unique

crash_sf <- crash %>%
  filter(
    !is.na(Longitude),
    !is.na(Latitude),
    Longitude < 0,
    Latitude > 0
  ) %>%
  st_as_sf(
    coords = c("Longitude", "Latitude"),
    crs = 4326,
    remove = FALSE
  )

# ------------------------------------------------------------
# Texas state boundary
# ------------------------------------------------------------

tx_boundary <- states(cb = TRUE, year = 2020, class = "sf") %>%
  filter(STUSPS == "TX") %>%
  st_transform(3083) %>%
  st_make_valid()

# Project crash points
crash_tx <- crash_sf %>%
  st_transform(3083)

# Keep only crashes inside Texas
crash_tx <- st_filter(crash_tx, tx_boundary)

# ------------------------------------------------------------
# Create approximately 80,000 square grids over Texas
# ------------------------------------------------------------

target_grid_count <- 80000

# Texas area in square meters
tx_area_m2 <- as.numeric(st_area(tx_boundary))

# Cell size in meters
cell_size_m <- sqrt(tx_area_m2 / target_grid_count)

cell_size_m

# ------------------------------------------------------------
# Build square fishnet grid
# ------------------------------------------------------------

tx_grid <- st_make_grid(
  tx_boundary,
  cellsize = cell_size_m,
  square = TRUE
) %>%
  st_sf(grid_id = seq_along(.), geometry = .)

# Clip grid to Texas boundary
tx_grid <- st_intersection(tx_grid, tx_boundary) %>%
  select(grid_id, geometry)

# Check actual number of cells
nrow(tx_grid)

# ------------------------------------------------------------
# Spatially join crashes to grid cells
# ------------------------------------------------------------

crash_grid_join <- st_join(
  crash_tx,
  tx_grid,
  join = st_within,
  left = FALSE
)

# Count crashes per grid
grid_crash_counts <- crash_grid_join %>%
  st_drop_geometry() %>%
  count(grid_id, name = "crash_count")

# Join crash counts back to grid
tx_grid_crash_map <- tx_grid %>%
  left_join(grid_crash_counts, by = "grid_id") %>%
  mutate(
    crash_count = ifelse(is.na(crash_count), 0, crash_count)
  )

summary(tx_grid_crash_map$crash_count)

quantile(
  tx_grid_crash_map$crash_count,
  probs = c(0, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99, 1),
  na.rm = TRUE
)

tx_grid_crash_map %>%
  st_drop_geometry() %>%
  count(crash_count, sort = FALSE) %>%
  arrange(crash_count) %>%
  head(30)

# ------------------------------------------------------------
# Data-driven bins: zero separate + Jenks for positive counts
# ------------------------------------------------------------

positive_counts <- tx_grid_crash_map$crash_count[
  !is.na(tx_grid_crash_map$crash_count) &
    tx_grid_crash_map$crash_count > 0
]

n_classes <- 6

jenks_breaks <- classIntervals(
  positive_counts,
  n = n_classes,
  style = "jenks"
)$brks

# Clean breaks
breaks_final <- sort(unique(c(0, round(jenks_breaks), Inf)))

breaks_final

# ------------------------------------------------------------
# Create readable labels
# ------------------------------------------------------------

make_bin_labels <- function(breaks) {
  
  labels <- character(length(breaks) - 1)
  
  for (i in seq_along(labels)) {
    
    lower <- breaks[i]
    upper <- breaks[i + 1]
    
    if (lower == 0 && upper == 1) {
      labels[i] <- "0"
    } else if (is.infinite(upper)) {
      labels[i] <- paste0(lower + 1, "+")
    } else {
      labels[i] <- paste0(lower + 1, "–", upper)
    }
  }
  
  labels
}

bin_labels <- make_bin_labels(breaks_final)

bin_labels

tx_grid_crash_map <- tx_grid_crash_map %>%
  mutate(
    crash_count_bin = cut(
      crash_count,
      breaks = breaks_final,
      labels = bin_labels,
      include.lowest = TRUE,
      right = TRUE
    )
  )

# Check number of grid cells in each bin
tx_grid_crash_map %>%
  st_drop_geometry() %>%
  count(crash_count_bin)

display.brewer.pal(9, "YlOrRd")

display.brewer.pal(9, "OrRd")

display.brewer.pal(9, "Reds")

palette_name <- "YlOrRd"

n_bins <- length(levels(tx_grid_crash_map$crash_count_bin))

map_colors <- colorRampPalette(
  brewer.pal(9, palette_name)
)(n_bins)

# Make zero-crash grid cells light gray
map_colors[1] <- "gray92"

# ------------------------------------------------------------
# Static crash-count grid map
# ------------------------------------------------------------

ggplot() +
  geom_sf(
    data = tx_grid_crash_map,
    aes(fill = crash_count_bin),
    color = NA
  ) +
  scale_fill_manual(
    values = map_colors,
    na.value = "transparent",
    name = "Crash Count",
    drop = FALSE
  ) +
  labs(
    title = "Texas Pedestrian Crashes by Square Grid",
    subtitle = paste0(
      "Crash counts aggregated to approximately ",
      format(target_grid_count, big.mark = ","),
      " square grid cells"
    ),
    caption = "Data: Texas crash records; grid generated from Texas state boundary"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    plot.subtitle = element_text(size = 12),
    plot.caption = element_text(size = 9, color = "gray40"),
    legend.position = "right",
    legend.title = element_text(face = "bold"),
    legend.text = element_text(size = 9),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.title = element_blank(),
    axis.ticks = element_blank()
  )

tx_grid_crash_only <- tx_grid_crash_map %>%
  filter(crash_count > 0)

ggplot() +
  geom_sf(
    data = tx_boundary,
    fill = "gray97",
    color = "gray50",
    linewidth = 0.25
  ) +
  geom_sf(
    data = tx_grid_crash_only,
    aes(fill = crash_count_bin),
    color = NA
  ) +
  scale_fill_manual(
    values = map_colors,
    na.value = "transparent",
    name = "Crash Count",
    drop = FALSE
  ) +
  labs(
    title = "Texas Pedestrian Crash Concentrations by Square Grid",
    subtitle = "Only grid cells with at least one crash are shown",
    caption = "Data: Texas crash records; grid generated from Texas state boundary"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    plot.subtitle = element_text(size = 12),
    legend.position = "right",
    legend.title = element_text(face = "bold"),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    axis.title = element_blank(),
    axis.ticks = element_blank()
  )

ggplot() +
  geom_sf(
    data = tx_boundary,
    fill = "gray97",
    color = "gray50",
    linewidth = 0.25
  ) +
  geom_sf(
    data = tx_grid_crash_only,
    aes(fill = crash_count_bin),
    color = NA
  ) +
  geom_sf(
    data = crash_tx,
    size = 0.15,
    alpha = 0.20,
    color = "black"
  ) +
  scale_fill_manual(
    values = map_colors,
    na.value = "transparent",
    name = "Crash Count",
    drop = FALSE
  ) +
  labs(
    title = "Texas Pedestrian Crashes and Square-Grid Hotspots",
    subtitle = "Crash points overlaid on grid-level crash counts"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    plot.subtitle = element_text(size = 12),
    legend.position = "right",
    panel.grid = element_blank(),
    axis.text = element_blank(),
    axis.title = element_blank(),
    axis.ticks = element_blank()
  )

# ------------------------------------------------------------
# Prepare data for leaflet
# ------------------------------------------------------------

tx_grid_leaflet <- tx_grid_crash_map %>%
  st_transform(4326)

tx_grid_leaflet_crash_only <- tx_grid_leaflet %>%
  filter(crash_count > 0)

crash_leaflet <- crash_tx %>%
  st_transform(4326)

tx_boundary_leaflet <- tx_boundary %>%
  st_transform(4326)

tx_grid_leaflet_crash_only <- tx_grid_leaflet_crash_only %>%
  mutate(
    popup_text = paste0(
      "<b>Grid ID:</b> ", grid_id, "<br>",
      "<b>Crash Count:</b> ", crash_count, "<br>",
      "<b>Crash Count Bin:</b> ", crash_count_bin
    )
  )

pal_leaflet <- colorFactor(
  palette = map_colors,
  domain = levels(tx_grid_crash_map$crash_count_bin),
  na.color = "transparent"
)


# ------------------------------------------------------------
# Interactive leaflet map
# ------------------------------------------------------------

leaflet_grid_map <- leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  
  addPolygons(
    data = tx_boundary_leaflet,
    fill = FALSE,
    color = "gray40",
    weight = 1,
    opacity = 0.8,
    group = "Texas Boundary"
  ) %>%
  
  addPolygons(
    data = tx_grid_leaflet_crash_only,
    fillColor = ~pal_leaflet(crash_count_bin),
    fillOpacity = 0.75,
    color = "white",
    weight = 0.15,
    opacity = 0.7,
    popup = ~popup_text,
    label = ~paste0("Grid ", grid_id, ": ", crash_count, " crashes"),
    highlightOptions = highlightOptions(
      weight = 2,
      color = "#333333",
      fillOpacity = 0.90,
      bringToFront = TRUE
    ),
    group = "Crash Grid"
  ) %>%
  
  addCircleMarkers(
    data = crash_leaflet,
    lng = ~Longitude,
    lat = ~Latitude,
    radius = 1.5,
    stroke = FALSE,
    fillOpacity = 0.25,
    color = "black",
    popup = ~paste0(
      "<b>Crash ID:</b> ", Crash_ID, "<br>",
      "<b>Latitude:</b> ", Latitude, "<br>",
      "<b>Longitude:</b> ", Longitude
    ),
    group = "Crash Points"
  ) %>%
  
  addLegend(
    position = "bottomright",
    pal = pal_leaflet,
    values = levels(tx_grid_crash_map$crash_count_bin),
    title = "Crash Count",
    opacity = 0.8
  ) %>%
  
  addLayersControl(
    overlayGroups = c("Texas Boundary", "Crash Grid", "Crash Points"),
    options = layersControlOptions(collapsed = FALSE)
  )

leaflet_grid_map

# Save spatial grid
#st_write(
#  tx_grid_crash_map,
 # "tx_pedestrian_crashes_80000_grid.gpkg",
 # delete_dsn = TRUE
#)

# Save non-spatial table
#tx_grid_crash_map %>%
 # st_drop_geometry() %>%
 # write_csv("tx_pedestrian_crashes_80000_grid_counts.csv")