Lat-Long to Map V01

setwd("C:/Users/mvx13/OneDrive - Texas State University/Hackathon_Aaqib/Papers/TRB2025_Papers/Ambulance_Paper/Model Run by AH/Data Preparation")
library(data.table)
library(ggplot2)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:data.table':
## 
##     between, first, last

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(hexbin)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ✔ readr     2.1.5

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::between()     masks data.table::between()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ dplyr::first()       masks data.table::first()
## ✖ lubridate::hour()    masks data.table::hour()
## ✖ lubridate::isoweek() masks data.table::isoweek()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ dplyr::last()        masks data.table::last()
## ✖ lubridate::mday()    masks data.table::mday()
## ✖ lubridate::minute()  masks data.table::minute()
## ✖ lubridate::month()   masks data.table::month()
## ✖ lubridate::quarter() masks data.table::quarter()
## ✖ lubridate::second()  masks data.table::second()
## ✖ purrr::transpose()   masks data.table::transpose()
## ✖ lubridate::wday()    masks data.table::wday()
## ✖ lubridate::week()    masks data.table::week()
## ✖ lubridate::yday()    masks data.table::yday()
## ✖ lubridate::year()    masks data.table::year()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(sf)

## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

library(ggspatial)
library(tidygeocoder)
library(maps)

## 
## Attaching package: 'maps'
## 
## The following object is masked from 'package:purrr':
## 
##     map

dat0 <- fread("Ambulance_LL.csv")
dim(dat0)

## [1] 3785    5

head(dat0)

##      Type          Crash_ID2 Latitude Longitude Crash_Sev_IDOri
##    <char>             <char>    <num>     <num>          <char>
## 1: single TX_2017_15536208_1 29.86953 -95.40760             PDO
## 2: single TX_2017_15543821_1 32.74830 -97.50395             PDO
## 3: single TX_2017_15545943_1 33.43791 -97.76280             PDO
## 4: single TX_2017_15573323_1 27.78452 -97.41616             PDO
## 5: single TX_2017_15595247_1 32.82985 -96.95882             PDO
## 6: single TX_2017_15607733_1 30.65527 -95.83986             PDO

dat0 = dat0 %>% mutate_if(is.character, as.factor)
dat01= subset(dat0, Type=="single")
dat02= subset(dat0, Type=="multiple")


state_map_data <- map('state', fill = TRUE, plot = FALSE) %>% 
 st_as_sf()

Single

# Assuming dat0 is your dataframe
dat01 <- dat01 %>%
  mutate(Crash_Sev_IDOri = factor(Crash_Sev_IDOri, levels = c('KA', 'BC', 'PDO')))

ggplot() +
  geom_sf(data = state_map_data) +
  stat_density2d(data = dat01, aes(x = Longitude, y = Latitude, fill = ..level..), geom = "polygon", alpha = 0.5) +
  geom_point(data = dat01, aes(x = Longitude, y = Latitude), alpha = 0.5, color = "#798E87") +
  coord_sf(xlim = c(-106.5, -93.6), ylim = c(25.5, 36.2)) +
  facet_wrap(~ Crash_Sev_IDOri, nrow = 1) +
  theme_bw(base_size = 14) +
  scale_fill_gradientn(colours = c("#CDC08C", "#9C964A"), name = "Density", na.value = NA) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

## Warning: The dot-dot notation (`..level..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(level)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Multiple

# Assuming dat0 is your dataframe
dat02 <- dat02 %>%
  mutate(Crash_Sev_IDOri = factor(Crash_Sev_IDOri, levels = c('KA', 'BC', 'PDO')))

ggplot() +
  geom_sf(data = state_map_data) +
  stat_density2d(data = dat02, aes(x = Longitude, y = Latitude, fill = ..level..), geom = "polygon", alpha = 0.5) +
  geom_point(data = dat02, aes(x = Longitude, y = Latitude), alpha = 0.5, color = "#798E87") +
  coord_sf(xlim = c(-106.5, -93.6), ylim = c(25.5, 36.2)) +
  facet_wrap(~ Crash_Sev_IDOri, nrow = 1) +
  theme_bw(base_size = 14) +
  scale_fill_gradientn(colours = c("#CDC08C", "#9C964A"), name = "Density", na.value = NA) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))