Tauentzienstraße
Yogada Joshi
library(sf) |> suppressMessages()## Warning: package 'sf' was built under R version 4.3.3library(dplyr) |> suppressMessages()## Warning: package 'dplyr' was built under R version 4.3.2library(units) |> suppressMessages()## Warning: package 'units' was built under R version 4.3.3library(giscoR) |> suppressMessages()## Warning: package 'giscoR' was built under R version 4.3.3library(tmap) |> suppressMessages()## Warning: package 'tmap' was built under R version 4.3.3library(tidyr) |> suppressMessages()## Warning: package 'tidyr' was built under R version 4.3.3library(here) |> suppressMessages()## Warning: package 'here' was built under R version 4.3.3library(sf) |> suppressMessages()
library(ggplot2) |> suppressMessages()
library(lubridate) |> suppressMessages()## Warning: package 'lubridate' was built under R version 4.3.3library(mapview) |> suppressMessages()## Warning: package 'mapview' was built under R version 4.3.3library(readr) |> suppressMessages()## Warning: package 'readr' was built under R version 4.3.2Basics
# after checking lanes data, plot the linestring
# step 3: Plotting
plot(st_geometry(Tauentzienstraße), col = "blue")
tmap_mode("view")## ℹ tmap mode set to "view".tm_shape(Tauentzienstraße) + tm_lines(col = "blue", lwd = 2)# project to flat CRS
Tauentzienstraße_proj <- st_transform(Tauentzienstraße,25833) # UTM Zone 33N (Berlin)# Lenght
st_length(Tauentzienstraße_proj)## Units: [m]
## [1] 380.463277 114.491661 98.043728 48.125112 42.575629 113.946236
## [7] 8.100131 26.792154 60.952711 50.273855 1.869662 55.460338
## [13] 23.845874 15.657152sum(st_length(Tauentzienstraße_proj)) # = 1040.598 [m], note units not transformed## 1040.598 [m]sum(st_length(Tauentzienstraße_proj))/1000 # = 1.040598 [km]## 1.040598 [m]Constructing the polygon
::: Each lane = 3.5 meters
Sidewalks = 1.5 meters on each side = 3.0 meters total
Buffer distance = (number of lanes × 3.5) + 3.0, divided by 2 (since buffer expands outward) :::
# Step 2: Clean and ensure 'lanes' column is numeric
Tauentzienstraße_proj <- Tauentzienstraße_proj %>%
mutate(lanes = as.numeric(lanes))
# Step 3: Calculate half-width buffer per feature
# Assuming 1.5m average lane width on either side of the street to form the polygon/area of each street to geolocate/spatially subset accidents
Tauentzienstraße_proj <- Tauentzienstraße_proj %>%
mutate(buffer_width = ((lanes * 3.5) + 3.0) / 2)
# Step 4: Apply variable buffer
Tauentzienstraße_buffer <- st_buffer(Tauentzienstraße_proj, dist = Tauentzienstraße_proj$buffer_width)
# Optional: Combine into a single polygon
Tauentzienstraße_union <- st_union(Tauentzienstraße_buffer)
# Step 5: Transform back to WGS84 if needed
Tauentzienstraße_union_wgs84 <- st_transform(Tauentzienstraße_union, 4326)
st_is_valid(Tauentzienstraße_union_wgs84)## [1] TRUEprint("Tauentzienstraße_union_wgs84 has been converted from a multi-linestring to one single polygon")## [1] "Tauentzienstraße_union_wgs84 has been converted from a multi-linestring to one single polygon"tmap_mode("view")## ℹ tmap mode set to "view".tm_shape(Tauentzienstraße_union_wgs84) +
tm_polygons(col = "blue", fill_alpha = 0.5)import accidents and transpose/spatially subset
# 74,838 (Jan 2018- Dec 2023) 6 years Berlin accidents data
accidents <- readRDS("accidents_2018_2023.rds")# Convert to sf object (assuming the CSV has WGS84 coordinates)
accidents_sf <- st_as_sf(accidents,
coords = c("XGCSWGS84", "YGCSWGS84"), # Longitude, Latitude
crs = 4326) # Set Coordinate Reference System (WGS84)
#accidents_sf is a geometrical spatial object now, different from "accidents" a statistical dataframe.# Ensure both layers are in the same CRS (WGS84)
accidents_sf <- st_transform(accidents_sf, st_crs(Tauentzienstraße_union_wgs84))
st_crs(accidents_sf) # Check CRS of accident points## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]st_crs(Tauentzienstraße_union_wgs84) # Check CRS of polygon## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]# asks R which points are within Tauentzienstraße_union_wgs84
street_matrix <- st_intersects(accidents_sf, Tauentzienstraße_union_wgs84, sparse = FALSE)
# returns point-wise matrix of matching both accident coordinates and the polygon area coordinates
head(street_matrix)## [,1]
## [1,] FALSE
## [2,] FALSE
## [3,] FALSE
## [4,] FALSE
## [5,] FALSE
## [6,] FALSE# Identify which points intersect with any polygon
intersecting_street_rows <- apply(street_matrix, 1, any)Tauentzienstraße_accidents <- accidents_sf[intersecting_street_rows, ]write_csv(st_drop_geometry(Tauentzienstraße_accidents), "accident_observations/Tauentzienstraße_accidents_2018_2023.csv")tmap_mode("view")## ℹ tmap mode set to "view".tm_shape(Tauentzienstraße_accidents) +
tm_polygons(col = "lightgreen", alpha = 0.5) +
tm_shape(Tauentzienstraße_accidents) +
tm_dots(col = "red", size = 0.5)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').[v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.