{osmdata}でエトセトラ(2/3) 極座標グラフの作成

{osmdata}で取得した道路ネットワークで遊ぶ第2回。今回は道路の角度を求め、極座標グラフを描いてみる。

道路データの細分化
角度・距離の計算
極座標グラフの作成

道路データの細分化

前回取得した道路ネットワークデータは1行に複数の区間が含まれているものがあり、そのままでは角度を計算できない。たとえば若洲には一周しているようなレコードがある。

# とある若洲の道路（結果は1行）
ways %>% 
  dplyr::filter(osm_id == 175401748)

Simple feature collection with 1 feature and 9 fields
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: 139.8333 ymin: 35.61354 xmax: 139.8392 ymax: 35.62948
geographic CRS: WGS 84
# A tibble: 1 x 10
  osm_id name  highway waterway aerialway barrier man_made z_order other_tags
* <chr>  <chr> <chr>   <chr>    <chr>     <chr>   <chr>      <int> <chr>     
1 17540… NA    service NA       NA        NA      NA             0 "\"golf\"…
# … with 1 more variable: geometry <LINESTRING [°]>

# 地図上で見ると複数のLINESTRINGの集合に見える
ways %>% 
  dplyr::filter(osm_id == 175401748) %>% 
  leaflet::leaflet(options = leaflet::leafletOptions(zoomControl = FALSE)) %>%
  leaflet::addTiles() %>% 
  leaflet::addPolylines()

f:id:tak95:20200901025744p:plain

そのためLINESTRINGを構成する座標から複数の区間に分割する。

# LINESTRINGを分割
ways.segment  <-  
  ways %>% 
  # LINESTRINGの区間をリストに格納
  dplyr::mutate(coords = purrr::map(geometry, function(x) {
    x %>% 
      # 座標を取得
      sf::st_coordinates() %>% 
      dplyr::as_tibble() %>% 
      # LINESTRINGの起点
      dplyr::select(X, Y) %>% 
      dplyr::rename(lng1 = X,
                    lat1 = Y) %>% 
      # LINESTRINGの終点
      dplyr::mutate(lng2 = dplyr::lead(lng1),
                    lat2 = dplyr::lead(lat1)) %>% 
      # 終点がない行を除く
      dplyr::filter(!is.na(lng2))
  })) %>% 
  # ジオメトリ情報を削除
  sf::st_set_geometry(NULL) %>% 
  # リストを展開
  tidyr::unnest(cols  = coords) %>% 
  # 区間のLINESTRINGを作成
  dplyr::mutate(geometry = purrr::pmap(
    list(lng1, lat1, lng2, lat2),
    function(lng1, lat1, lng2, lat2) {
      c(lng1, lat1, lng2, lat2) %>% 
        matrix(., ncol = 2, byrow = TRUE) %>% 
        sf::st_linestring()
    })
  ) %>% 
  dplyr::mutate(geometry = sf::st_as_sfc(geometry, crs = 4326)) %>% 
  sf::st_sf()

角度・距離の計算

道路を区間毎に分割できたので角度を求める。Pythonのosmnxモジュールのget_bearing()を参考*1に角度を求める関数を定義する。

get_bearing  <- function(lng1, lat1, lng2, lat2) {

  rad2deg <- function(rad) {(rad * 180) / (pi)}
  deg2rad <- function(deg) {(deg * pi) / (180)}

  # get latitudes and the difference in longitude, as radians
  lat1 <- deg2rad(lat1)
  lat2 <- deg2rad(lat2)
  diff_lng <- deg2rad(lng2 - lng1)

  # calculate initial bearing from -180 degrees to +180 degree
  x <-
    sin(diff_lng) * cos(lat2)
  y <-
    cos(lat1) * sin(lat2) - (sin(lat1) * cos(lat2) * cos(diff_lng))
  initial_bearing <- 
    atan2(x, y) %>% 
    rad2deg()
  
  # normalize initial bearing to 0-360 degrees to get compass bearing
  bearing <-
    (initial_bearing + 360) %% 360
  
  return(bearing)
}

区間に分割した道路ネットワークに適用し、角度および距離の列を追加する。

ways.bearing  <-
  ways.segment %>% 
  # 角度を列に追加
  dplyr::mutate(bearing =  purrr::pmap_dbl(
    list(lng1, lat1,  lng2, lat2),
    get_bearing
  )) %>% 
  # 平面直角座標系に変換
  sf::st_transform(6677) %>% 
  # 距離を列に追加
  dplyr::mutate(dist = purrr::map_dbl(geometry, sf::st_length)) %>% 
  sf::st_transform(4326)

極座標グラフの作成

各道路の角度が求まったので極座標グラフを描いてみる。

ways.bearing %>% 
  # ジオメトリを削除
  sf::st_set_geometry(NULL) %>% 
  # 15度ごとに角度を丸める
  dplyr::mutate(bearing = bearing %/% 15 * 15) %>% 
  # 丸めた角度毎に距離を合計
  dplyr::group_by(bearing) %>% 
  dplyr::summarise(dist = sum(dist)) %>%  
  dplyr::ungroup() %>% 
  # 極座標グラフの描画
  ggplot2::ggplot(ggplot2::aes(x = bearing, y = dist)) + 
  ggplot2::geom_bar(stat="identity", fill = RColorBrewer::brewer.pal(3, "Dark2")[1]) + 
  ggplot2::coord_polar(start = pi/180 * (-7.5)) +
  ggplot2::scale_x_continuous(breaks = seq(0, 360, 45)) +
  ggplot2::ggtitle(stringr::str_c(dplyr::coalesce(city$N03_003, ""), city$N03_004)) +
  ggthemes::theme_gdocs(base_family = "HiraKakuPro-W3") +
  ggplot2::scale_fill_brewer(palette = "Set1") +
  ggplot2::theme(axis.title.x = ggplot2::element_blank(),
                 axis.title.y = ggplot2::element_blank(),
                 axis.text.y = ggplot2::element_blank())