避難場所への経路をRで探索してみよう

wget https://download.geofabrik.de/asia/japan/chubu-latest.osm.pbf

docker run -t -v "${PWD}:/data" osrm/osrm-backend osrm-extract \
  -p /opt/foot.lua /data/chubu-latest.osm.pbf
docker run -t -v "${PWD}:/data" osrm/osrm-backend osrm-partition \
  /data/chubu-latest.osrm
docker run -t -v "${PWD}:/data" osrm/osrm-backend osrm-customize \
  /data/chubu-latest.osrm

docker run -t -i -p 5000:5000 -v "${PWD}:/data" \
  osrm/osrm-backend osrm-routed \
  --algorithm mld /data/chubu-latest.osrm

> library(osrm)
> 
> options(osrm.server = "http://localhost:5000/",
+         osrm.profile = "foot")

> # 金沢駅の位置
> kanazawa_eki <- c(136.648175, 36.578044)
> # 兼六園の位置
> kenrokuen <- c(136.662603, 36.562085)

> library(leaflet)
> leaflet() |>
+   addTiles() |>
+   setView(lng = (kanazawa_eki[1] + kenrokuen[1]) / 2,
+           lat = (kanazawa_eki[2] + kenrokuen[2]) / 2,
+           zoom = 14) |>
+   addCircleMarkers(lng = kanazawa_eki[1], lat = kanazawa_eki[2],
+                    label = "金沢駅", color = "blue", opacity = 0.9,
+                    labelOptions = labelOptions(noHide = TRUE)) |>
+   addCircleMarkers(lng = kenrokuen[1], lat = kenrokuen[2],
+                    label = "兼六園", color = "red", opacity = 0.9,
+                    labelOptions = labelOptions(noHide = TRUE))

> (route <- osrmRoute(src = kanazawa_eki, dst = kenrokuen))
Simple feature collection with 1 feature and 4 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 136.6483 ymin: 36.56208 xmax: 136.6628 ymax: 36.57815
Geodetic CRS:  WGS 84
        src dst duration distance                       geometry
src_dst src dst 33.58167   2.7222 LINESTRING (136.6483 36.577...

> leaflet() |>
+   addTiles() |>
+   setView(lng = (kanazawa_eki[1] + kenrokuen[1]) / 2,
+           lat = (kanazawa_eki[2] + kenrokuen[2]) / 2,
+           zoom = 14) |>
+   addCircleMarkers(lng = kanazawa_eki[1], lat = kanazawa_eki[2],
+                    label = "金沢駅", color = "blue", opacity = 0.9,
+                    labelOptions = labelOptions(noHide = TRUE)) |>
+   addCircleMarkers(lng = kenrokuen[1], lat = kenrokuen[2],
+                    label = "兼六園", color = "red", opacity = 0.9,
+                    labelOptions = labelOptions(noHide = TRUE)) |>
+   addPolylines(data = route, col = "purple", weight = 4, opacity = 0.9)

> library(tidyverse)
> library(sf)
> 
> # Data
> data_dir <- "data"
> data_file <- "ishikawa-central-evacuation-space.csv"
> evac_data <- read_csv(file.path(data_dir, data_file))

> eki_df <- data.frame(lng = kanazawa_eki[1],
+                      lat = kanazawa_eki[2])
> evac_df <- evac_data |>
+   dplyr::select(lng = `経度`, lon = `緯度`)
> results <- osrmTable(src = eki_df,
+                      dst = evac_df,
+                      measure = "duration")
> results_df <- data.frame(id = 1:nrow(evac_data),
+                          duration = results$durations[1, ]) |>
+   dplyr::arrange(duration)

> evac_routes <- purrr::map_df(
+   1:5,
+   \(i)
+   osrmRoute(src = kanazawa_eki,
+             dst = evac_df[results_df$id[i], ]) |>
+     dplyr::mutate(rank = i)
+ )

> # カラーパレットの作成関数
> pal_fun <- colorNumeric("viridis", domain = 1:6)
> 
> leaflet() |>
+   addTiles() |>
+   setView(lng = kanazawa_eki[1], lat = kanazawa_eki[2],
+           zoom = 16) |>
+   addCircleMarkers(data = evac_data, 
+                    lng = ~`経度`, lat = ~`緯度`,
+                    label = ~`名称`,
+                    color = "green", opacity = 0.5) |>
+   addPolylines(data = evac_routes,
+                color = pal_fun(1:5),
+                weight = 6, opacity = 0.8) |>
+   addLegend(data = evac_routes,
+             colors = pal_fun(1:5),
+             labels = paste0("避難経路", 1:5),
+             opacity = 1,
+             position = "topright")

> evac_sf <- evac_data |>
+   sf::st_as_sf(coords = c("経度", "緯度"), crs = 4326) # WGS84
> evac_routes <- evac_routes |>
+   dplyr::mutate(rank = ordered(rank))
> 
> # 金沢駅中心にほぼ2km四方
> lng <- kanazawa_eki[1] + c(-0.0125, 0.0125)
> lat <- kanazawa_eki[2] + c(-0.00833, 0.00833)
> bb <- st_bbox(c(xmin = lng[1], xmax = lng[2],
+                 ymin = lat[1], ymax = lat[2]), crs = 4326) # WGS84

> library(tmap)
> library(tmap.mapgl)
> 
> tmap_mode("maplibre")
> tm_shape(evac_sf, bbox = bb) +
+   tm_maplibre() +
+   tm_symbols(fill = "red", fill_alpha = 0.67)

> tm_shape(evac_sf, bbox = bb) +
+   tm_maplibre(pitch = 60) +
+   tm_symbols(fill = "red", fill_alpha = 0.67) +
+   tm_shape(evac_routes) +
+   tm_lines(col = "rank", lwd = 5)