メモ：「道路を方角ごとに塗り分けると、その街のでき方がわかる :: デイリーポータルZ」を地理院地図ベクトルタイルとggplot2で

この記事を読んで、ベクトルタイルのよい使い方（ベクトルなので方向を取り出すことができる）だと思ったのでggplot2でもやってみた。もうすでに誰かやってそうだなと思いつつ...

元ネタになっている記事はこれ。

library(ggplot2)
#> Warning: replacing previous import 'vctrs::data_frame' by 'tibble::data_frame'
#> when loading 'dplyr'
library(patchwork)
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(dplyr, warn.conflicts = FALSE)

sec <- function(x) {
  1 / cos(x)
}

lonlat2xy <- function(lat_deg, lon_deg, zoom) {
  n <- 2^zoom
  
  x <- (n * (lon_deg + 180)) %/% 360
  lat_rad <- lat_deg * pi / 180
  y <- (n * (1 - log(tan(lat_rad) + sec(lat_rad)) / pi)) %/% 2
  
  list(x = x, y = y)
}

zoom_level <- 14
# 国技館らへん
xy_base <- lonlat2xy(35.697, 139.793, zoom_level)

# 上下左右も
xy <- expand.grid(x = xy_base$x + -1:1, y = xy_base$y + -1:1)

v <- purrr::pmap(xy, function(x, y) {
  url <- glue::glue("https://cyberjapandata.gsi.go.jp/xyz/experimental_bvmap/{zoom_level}/{x}/{y}.pbf")
  protolite::read_mvt_sf(url)
})

road <- purrr::map_dfr(v, "road")

ggplot(road) +
  geom_sf() +
  theme_minimal() +
  labs(caption = "地理院地図ベクトルタイルを加工して作成")

# 少し MULTILINESTRING が混じっているが少数なので無視
table(as.character(st_geometry_type(road)))
#> 
#>      LINESTRING MULTILINESTRING 
#>           33774               6

# こんな感じで角度を計算できるはず（あってる？）
m <- as.matrix(road$geometry[[1]])
atan2(m[2, 1] - m[1, 1], m[2, 2] - m[1, 2])
#> [1] 1.380848

# 実験
sfc <- st_sfc(purrr::map(2 * pi * 1:100 / 100, ~ st_linestring(rbind(c(0, 0), c(sin(.x), cos(.x))))))
d <- st_sf(geometry = sfc)
d %>% 
  mutate(
    angle = purrr::map_dbl(geometry, ~ {
      m <- as.matrix(.x)
      atan2(m[2, 1] - m[1, 1], m[2, 2] - m[1, 2]) / pi * 180
    }),
    direction = abs(angle %% 90 - 45)
  ) %>% 
  ggplot()+
  geom_sf(aes(colour = direction)) +
  theme_minimal() +
  scale_colour_viridis_c(option = "B")

# あってそうなので元の地図にも同じ計算を適用してみる
road <- road %>% 
  mutate(
    angle = purrr::map_dbl(geometry, ~ {
      m <- as.matrix(.x)
      atan2(m[2, 1] - m[1, 1], m[2, 2] - m[1, 2]) / pi * 180
    }),
    direction = abs(angle %% 90 - 45)
  )

ggplot(road) +
  geom_sf(aes(colour = direction)) +
  theme_minimal() +
  scale_colour_viridis_c(option = "B") +
  labs(caption = "地理院地図ベクトルタイルを加工して作成")

^{Created on 2020-08-15 by the reprex package (v0.3.0)}