::: The final interactive map shows all the road crashes that took place in the 6-year analysis timeline of the study titled “The Effects of Car-free Zones on Road Crashes and their Displacement: A quasi-experimental study of Berlin’s Flaniermeile Friedrichstraße”, supervised by Prof. Dr. Christian Traxler, Hertie School. The map has been produced as part of the study for readers to understand the status of accidents when there was treatment and when the treatment was off.

Total crashes during treatments = 7, no-treatment = 41
Timeline:

Jan 2018–Aug 2020 - Pre-treatment-Phase-I
30 crashes in 2 years 7 months

Sep 2020–Oct 2021 - Treatment Phase-I (Only bikes, e-bikes and pedestrians allowed)
3 crashes in 1 year 2 months

Nov 2021–Jan 2023 - Pre-treatment Phase-II
5 crashes in 1 year 2 months

Feb 2023–Jun 2023 - Treatment Phase-II (Only bikes, e-bikes and pedestrians allowed)
4 crashes in 5 months

Jul 2023–Dec 2023 - Post-treatment Phase-II
6 crashes in 6 months

Note: The road crashes occuring on the ends of the Promenade are not included i.e., the buffer extends at the arms (east and west) only, not North and South. Hence, crashes occuring on Leipzigerstraße and Französischestraße do not count into the treated sample of crashes.

Upon clicking on each crash, one can see the type of vehicles involved. :::

library(readr)
## Warning: package 'readr' was built under R version 4.3.2
library(ggplot2)
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.2
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.3.3
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.3.3
library(sf)
## Warning: package 'sf' was built under R version 4.3.3
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(lubridate)
## Warning: package 'lubridate' was built under R version 4.3.3
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

# --- 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 25m Buffer ---
print("Creating 25m buffer...")
## [1] "Creating 25m buffer..."
FF_Buffer25m_proj <- st_buffer(FF_Treated_proj, dist = 25)

# --- 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_Buffer25m_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: 65"
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", "darkgreen", "darkred", "green","yellow")
)

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) 
##                               36                                7 
## Cars allowed (Nov 2021–Jan 2023)  Treatment-1 (Sep 2020–Oct 2021) 
##                                8                                5 
##  Treatment-2 (Feb 2023–Jun 2023) 
##                                9
# --- 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_Buffer25m_wgs84 <- st_transform(FF_Buffer25m_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_Buffer25m_wgs84, color = "blue", weight = 1, fillOpacity = 0.05, label = "25m 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", "darkgreen", "darkred", "green","yellow"),
    labels = c("Cars-allowed", "Treatment-1","Cars-allowed", "Treatment-2","Cars-allowed"),
    title = "Accident Period Status"
    )

# Display map
print("Displaying map...")
## [1] "Displaying map..."
m

::: Code generated with the assistance of ChatGPT and Gemini.

Data source: Statistische Ämter des Bundes und der Länder. (n.d.). Unfallatlas.
https://unfallatlas.statistikportal.de/

:::