メモ：地理院地図ベクトルタイルをgeom_sf()でプロット

protoliteパッケージがmapbox vector tileを読めるようになったらしい、というツイートを見ました。

New version for #rstats protolite: Highly Optimized Protocol Buffer Serializers https://t.co/qEoFHkyUV9. This version adds support for reading mvt (Mapbox vector tiles) as 'sf' objects. Which other proto formats should we add?
— Jeroen Ooms (@opencpu) 2019年10月10日

なんか適当なやつで試してみよう、と思ってそういえば地理院地図のベクトルタイルがあるのを思い出したのでメモ。ちなみに、protoliteを使わなくても、ふつうにsf::read_sf()で読める気はします（GDALのコンパイル時のオプションによるかも）

ベクトルタイルは、以下のURL形式で配信されています。

https://cyberjapandata.gsi.go.jp/xyz/experimental_bvmap/{z}/{x}/{y}.pbf

適当に、ズームレベル7で緯度35°、経度139°の地点をタイル座標に変換します（コードはSlippy map tilenames - OpenStreetMap WikiのR版）。

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)
}

lonlat2xy(35, 139, 7)
#> $x
#> [1] 113
#> 
#> $y
#> [1] 50

得られたx、y、あとズームレベル（z）の数字を先程のURLに入れて、read_mvt_sf()でsfのオブジェクトとして読み込みます。

v <- protolite::read_mvt_sf("https://cyberjapandata.gsi.go.jp/xyz/experimental_bvmap/7/113/50.pbf")

ここでvはsfオブジェクトのリストになっており、さまざまなレイヤーを含んでいます（ズームレベルによってレイヤーの数は違う）。とりあえずざっくりプロットしてみるとこんな感じ。

library(ggplot2)

ggplot() +
  geom_sf(data = v$coastline) +
  geom_sf(data = v$waterarea, fill = "blue", alpha = 0.6) +
  geom_sf(data = v$river, colour = "blue") +
  geom_sf(data = v$railway, colour = "orange") +
  geom_sf(data = v$road, colour = "gray") +
  geom_sf(data = v$symbol) +
  # symbolには「甲府」みたいな市町村名がknjに入ってるものと、「松本空港」みたいな地物名がnameに入ってるデータが混在しています
  # 丁寧にやるならデータをsplitしたほうがいいですが、適当にプロットしたいだけなのでNAの警告が出るのは無視してknjもnameも全部プロットする
  geom_sf_text(data = v$symbol, aes(label = knj), hjust = 0, vjust = 1, nudge_x = 0.05) +
  geom_sf_text(data = v$symbol, aes(label = name), hjust = 0, vjust = 0, nudge_x = -0.05, colour = "red") +
  theme_minimal() +
  labs(caption = "地理院地図ベクトルタイルを加工して作成")
#> Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may
#> not give correct results for longitude/latitude data

#> Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may
#> not give correct results for longitude/latitude data
#> Warning: Removed 9 rows containing missing values (geom_text).
#> Warning: Removed 29 rows containing missing values (geom_text).

f:id:yutannihilation:20191016084703p:plain