# --- 1. Load Data ---
semi_final_treated <- read.csv("accidents_2018_2023.csv")
FF_Treated <- st_read("Treated_Street_Section_FF.gpkg")
# --- 2. Prepare Spatial Data ---
tryCatch({
# Keep original columns by setting remove = FALSE
accidents_sf <- st_as_sf(semi_final_treated, coords = c("XGCSWGS84", "YGCSWGS84"), crs = 4326, remove = FALSE)
}, error = function(e) {
stop("Error converting accidents to sf object: ", e$message)
})
if (!all(st_is_valid(FF_Treated))) {
print("Making FF_Treated polygon valid...")
FF_Treated <- st_make_valid(FF_Treated)
}
## [1] "Making FF_Treated polygon valid..."
projected_crs_epsg <- 25833
print(paste("Target projected CRS for buffering/filtering: EPSG", projected_crs_epsg))
## [1] "Target projected CRS for buffering/filtering: EPSG 25833"
print("Transforming geometries to projected CRS...")
## [1] "Transforming geometries to projected CRS..."
FF_Treated_proj <- st_transform(FF_Treated, projected_crs_epsg)
# Make sure accidents_sf retains the ist* columns after conversion if remove=FALSE wasn't set
accidents_sf_proj <- st_transform(accidents_sf, projected_crs_epsg)
# --- 3. Create 10m Buffer ---
print("Creating 60m buffer...")
## [1] "Creating 60m buffer..."
FF_Buffer10m_proj <- st_buffer(FF_Treated_proj, dist = 60)
# --- 4. Filter Accidents ---
print("Filtering accidents within the 25m buffer zone...")
## [1] "Filtering accidents within the 25m buffer zone..."
accidents_in_area_proj <- st_filter(accidents_sf_proj, FF_Buffer10m_proj)
print(paste("Number of accidents within FF_Treated + 25m buffer:", nrow(accidents_in_area_proj)))
## [1] "Number of accidents within FF_Treated + 25m buffer: 77"
if (nrow(accidents_in_area_proj) == 0) {
warning("No accidents found within the specified polygon and 25m buffer.")
stop("Stopping execution as no accidents were found in the area.")
}
# --- 5. Define Periods and Categorize Accidents ---
print("Defining time periods and categorizing accidents...")
## [1] "Defining time periods and categorizing accidents..."
accidents_in_area_proj <- accidents_in_area_proj %>%
mutate(
UJAHR = as.numeric(UJAHR),
UMONAT = as.numeric(UMONAT),
AccidentDate = tryCatch(ymd(paste0(UJAHR, "-", UMONAT, "-01")), error = function(e) NA, warning = function(w) NA)
) %>%
filter(!is.na(AccidentDate))
periods <- data.frame(
StartDate = ymd(c("2018-01-01", "2020-09-01", "2021-11-01", "2023-02-01", "2023-07-01")),
EndDate = ymd(c("2020-08-31", "2021-10-31", "2023-01-31", "2023-06-30", "2023-12-31")),
Status = c("Cars allowed", "Treatment-1", "Cars allowed", "Treatment-2", "Cars allowed"),
PeriodLabel = c("Cars allowed (Jan 2018–Aug 2020)", "Treatment-1 (Sep 2020–Oct 2021)", "Cars allowed (Nov 2021–Jan 2023)", "Treatment-2 (Feb 2023–Jun 2023)", "Cars allowed (Jul 2023–Dec 2023)"),
PeriodColor = c("red", "black", "red", "black", "red")
)
accidents_in_area_proj <- accidents_in_area_proj %>%
mutate(
PeriodLabel = case_when(
AccidentDate >= periods$StartDate[1] & AccidentDate <= periods$EndDate[1] ~ periods$PeriodLabel[1],
AccidentDate >= periods$StartDate[2] & AccidentDate <= periods$EndDate[2] ~ periods$PeriodLabel[2],
AccidentDate >= periods$StartDate[3] & AccidentDate <= periods$EndDate[3] ~ periods$PeriodLabel[3],
AccidentDate >= periods$StartDate[4] & AccidentDate <= periods$EndDate[4] ~ periods$PeriodLabel[4],
AccidentDate >= periods$StartDate[5] & AccidentDate <= periods$EndDate[5] ~ periods$PeriodLabel[5],
TRUE ~ "Other Period"
),
PeriodColor = case_when(
AccidentDate >= periods$StartDate[1] & AccidentDate <= periods$EndDate[1] ~ periods$PeriodColor[1],
AccidentDate >= periods$StartDate[2] & AccidentDate <= periods$EndDate[2] ~ periods$PeriodColor[2],
AccidentDate >= periods$StartDate[3] & AccidentDate <= periods$EndDate[3] ~ periods$PeriodColor[3],
AccidentDate >= periods$StartDate[4] & AccidentDate <= periods$EndDate[4] ~ periods$PeriodColor[4],
AccidentDate >= periods$StartDate[5] & AccidentDate <= periods$EndDate[5] ~ periods$PeriodColor[5],
TRUE ~ "grey"
),
Status = case_when(
AccidentDate >= periods$StartDate[1] & AccidentDate <= periods$EndDate[1] ~ periods$Status[1],
AccidentDate >= periods$StartDate[2] & AccidentDate <= periods$EndDate[2] ~ periods$Status[2],
AccidentDate >= periods$StartDate[3] & AccidentDate <= periods$EndDate[3] ~ periods$Status[3],
AccidentDate >= periods$StartDate[4] & AccidentDate <= periods$EndDate[4] ~ periods$Status[4],
AccidentDate >= periods$StartDate[5] & AccidentDate <= periods$EndDate[5] ~ periods$Status[5],
TRUE ~ "Other"
)
)
print("Accident counts per defined period:")
## [1] "Accident counts per defined period:"
print(table(accidents_in_area_proj$PeriodLabel))
##
## Cars allowed (Jan 2018–Aug 2020) Cars allowed (Jul 2023–Dec 2023)
## 43 8
## Cars allowed (Nov 2021–Jan 2023) Treatment-1 (Sep 2020–Oct 2021)
## 10 6
## Treatment-2 (Feb 2023–Jun 2023)
## 10
# --- 6. Prepare Data for Leaflet (Transform back & Add Involved Parties) ---
print("Transforming results back to WGS84 (EPSG:4326)...")
## [1] "Transforming results back to WGS84 (EPSG:4326)..."
accidents_in_area_wgs84 <- st_transform(accidents_in_area_proj, 4326)
# *** NEW: Check for required columns and create InvolvedParties string ***
print("Generating involved parties information for popups...")
## [1] "Generating involved parties information for popups..."
required_cols <- c("IstRad", "IstFuss", "IstPKW", "IstSonstige", "IstKrad", "IstGkfz")
# Ensure columns exist, if not create them as 0 (or handle error)
for (col_name in required_cols) {
if (!col_name %in% names(accidents_in_area_wgs84)) {
warning(paste("Column", col_name, "not found. Assuming 0 for all entries."))
accidents_in_area_wgs84[[col_name]] <- 0 # Add column as 0 if missing
}
# Ensure columns are numeric and replace NA with 0
if(any(is.na(accidents_in_area_wgs84[[col_name]]))){
accidents_in_area_wgs84[[col_name]][is.na(accidents_in_area_wgs84[[col_name]])] <- 0
}
accidents_in_area_wgs84[[col_name]] <- as.numeric(accidents_in_area_wgs84[[col_name]])
}
accidents_in_area_wgs84 <- accidents_in_area_wgs84 %>%
mutate(
InvolvedParties = paste0(
"<b>Involved:</b>",
# Add a line for each party type where the corresponding column is 1
ifelse(IstRad == 1, "<br>- Bicycle", ""),
ifelse(IstFuss == 1, "<br>- Pedestrian", ""),
ifelse(IstPKW == 1, "<br>- Car", ""),
ifelse(IstKrad == 1, "<br>- Motorcycle", ""),
ifelse(IstGkfz == 1, "<br>- Truck/Comm.", ""),
ifelse(IstSonstige == 1, "<br>- Other", ""),
# Add a "None" message if all flags are 0
ifelse(IstRad != 1 & IstFuss != 1 & IstPKW != 1 & IstKrad != 1 & IstGkfz != 1 & IstSonstige != 1, "<br>None specified", "")
)
)
# You can check the generated strings: head(accidents_in_area_wgs84$InvolvedParties)
# Transform polygons back to WGS84
FF_Treated_wgs84 <- st_transform(FF_Treated_proj, 4326)
FF_Buffer10m_wgs84 <- st_transform(FF_Buffer10m_proj, 4326)
# --- 7. Create Leaflet Map ---
print("Creating Leaflet map with period grouping and involved parties...")
## [1] "Creating Leaflet map with period grouping and involved parties..."
map_center_lng <- tryCatch({st_coordinates(st_centroid(st_union(FF_Treated_wgs84)))[1]}, error = function(e) { mean(st_bbox(FF_Treated_wgs84)[c(1, 3)]) })
map_center_lat <- tryCatch({st_coordinates(st_centroid(st_union(FF_Treated_wgs84)))[2]}, error = function(e) { mean(st_bbox(FF_Treated_wgs84)[c(2, 4)]) })
m <- leaflet(data = accidents_in_area_wgs84, options = leafletOptions(minZoom = 12)) %>%
addProviderTiles(providers$CartoDB.Positron, group = "Map") %>%
setView(lng = map_center_lng, lat = map_center_lat, zoom = 16) %>%
addPolygons(data = FF_Treated_wgs84, color = "black", weight = 3, fillOpacity = 0.1, label = "Friedrichstraße Polygon", group = "Boundaries") %>%
addPolygons(data = FF_Buffer10m_wgs84, color = "blue", weight = 1, fillOpacity = 0.05, label = "10m Buffer", group = "Boundaries")
# Add markers
m <- m %>% addCircleMarkers(
radius = 4,
color = ~PeriodColor,
fillOpacity = 0.7,
stroke = FALSE,
group = ~PeriodLabel,
# *** UPDATED POPUP ***
popup = ~paste0(
"<b>ID:</b> ", OBJECTID,
"<br><b>Date:</b> ", format(AccidentDate, "%Y-%m"),
"<br><b>Period:</b> ", PeriodLabel,
"<br><b>Hour:</b> ", USTUNDE,
"<br><b>Category:</b> ", UKATEGORIE,
"<br><b>Type:</b> ", UART,
"<br>", # Add a line break before involved parties
InvolvedParties # Add the pre-formatted string from the new column
)
)
period_groups <- unique(accidents_in_area_proj$PeriodLabel)
period_groups <- c(period_groups[period_groups != "Other Period"], period_groups[period_groups == "Other Period"])
m <- m %>% addLayersControl(
baseGroups = c("Map"),
overlayGroups = c("Boundaries", period_groups),
options = layersControlOptions(collapsed = FALSE)
)
m <- m %>% addLegend(
position = "bottomright",
colors = c("red", "black"),
labels = c("Control Period", "Treated Period"),
title = "Accident Period Status"
)
# Display map
print("Displaying map...")
## [1] "Displaying map..."
m