scooter-ehsa.R

library(ggplot2)

tmp <- tempdir()

clone_call <- glue::glue("
cd {tmp}
git clone [email protected]:marcosci/scooter_scraping.git")

# clone repo into temp dir
system(clone_call)

# get all data files
scoot_files <- list.files(
  file.path(tmp, "scooter_scraping", "data"),
  full.names = TRUE,
  pattern = "^20"
)

# create a look up table based on file names
lookups <- tibble::tibble(fname_full = scoot_files)|>
  dplyr::mutate(
    fname = basename(fname_full),
    date_time = substr(fname, 12, 30),
    date_time = lubridate::ymd_hms(date_time),
    company = substr(fname, 32, nchar(fname) - 4)
  )

# read in all files
raw_df <- readr::read_csv(scoot_files, col_types = "ddc", id = "fname_full")

# join raw table to lookup tbale
all_df <- dplyr::left_join(raw_df, lookups, by = "fname_full") |>
  dplyr::select(id, date_time, company, x, y)

# filter outliers
scoots <- all_df |> 
  dplyr::filter(company == "yoio",
                y < 50,
                x > 7.78,
                y < 48.03) |> 
  sf::st_as_sf(coords = c("x", "y"), crs = 4326)

# plot all scooter locations
ggplot(scoots) +
  geom_sf(alpha = 0.05)

# make a 20 x 20 grid around all points
grd <- sf::st_make_grid(scoots, n = c(20, 20)) |> 
  sf::st_as_sf() |> 
  dplyr::mutate(grid_id = dplyr::row_number())

# join id to points
pnts <- sf::st_join(scoots, grd)

# count points ensuring all combinations of dat and location have a value
pnts_grid <- pnts |>
  dplyr::as_tibble() |>
  dplyr::mutate(day = lubridate::floor_date(date_time, "day"),
                day = as.Date(day)) |>
  dplyr::filter(!is.na(grid_id)) |>
  dplyr::count(day, grid_id, .drop = FALSE) |>
  tidyr::complete(day, grid_id, fill = list(n = 0))

# remove unnecessary grids
grd_subset <- dplyr::semi_join(grd, pnts_grid)

# cast as spacetime
scoot_spt <-
  sfdep::spacetime(.data = pnts_grid,
                   .geometry = grd_subset,
                   .loc_col = "grid_id",
                   .time_col = "day")

# check if it is a spacetime cube
sfdep::is_spacetime_cube(scoot_spt)

# conduct emerging hotspot analysis
ehsa <- sfdep::emerging_hotspot_analysis(scoot_spt, "n")

# get freiburg data
fb_fpath <- "https://geoportal.freiburg.de/wfs/abi_gliederung/abi_gliederung?REQUEST=GetFeature&SRSNAME=EPSG:25832&SERVICE=WFS&VERSION=2.0.0&TYPENAMES=stadtteile&OUTPUTFORMAT=geojson"

freiburg <- sf::read_sf(fb_fpath) |> 
  sf::st_transform(crs = 4326) |> 
  # remove inner boundaries
  sf::st_union()
  
# join together
dplyr::left_join(ehsa, grd_subset, by = c("location" = "grid_id")) |>
  sf::st_as_sf() |>
  ggplot() +
  geom_sf(aes(fill = classification), lwd = 0.1, color = "black") +
  geom_sf(data = freiburg,
          lwd = 0.2, 
          color = "black",
          fill = NA) + 
  # specify colors
  scale_fill_manual(
    values = c(
      "no pattern detected" = "white", 
      "new coldspot" = "#84d6e3", 
      "oscilating coldspot" = "#6ba5c7",
      "sporadic coldspot" = "#5084a3",
      "consecutive coldspot" = "#386682",
      "sporadic hotspot" = "#db6381", 
      "consecutive hotspot" = "#a3485f")
  ) + 
  theme_light()