地物範囲のタイル座標を得る
id:yutannihilationさんが、ggplot2::geom_sf()
にOpen Street Map (OSM)のタイルを重ねるという面白い試みをしています。
OSMでは、画像データとして地図画像を配信していて、次のようなURLで参照されます。
ドメイン名の他、zoom
、x
のディレクトリがあり、ファイルはy
という名前で置かれています。ここでzoom
は地図のズームレベル、xとyはタイルの座標を示す、ピクセル座標と呼ばれる値です。OSMをはじめ、地図タイルの配信を行なっているサービスは同様にこの「XYZ地図タイル」形式を採用しています。
OSMに限らず、目的の地物を範囲に収めるタイルを用意するには、このピクセル座標を取得する必要があるわけですが、馴染みのある緯度経度からピクセル座標の変換は少し頭を使わなくてはいけません(私だけかもしれませんが)。今回の話題は、すでにyutannihilationさんの記事にも書かれていますが、緯度経度からピクセル座標を取得する処理をおさらいをするというものです(二番煎じです。本当にありがとうございました)。
基礎
緯度経度およびズームレベルからピクセル座標を求める処理は、OSMのWikiに書かれている通りです。計算式の詳細は参考にあげるページに詳しいのでそちらをご覧ください。
library(magrittr)
Wikiに書かれているコードをちょっと修正して、緯度経度とズームレベルを与えた時にピクセル座標とズームレベルをリストで返す関数を定義します。
deg2num <- function(longitude, latitude, zoom){ lat_rad <- latitude * pi /180 n <- 2.0 ^ zoom xtile <- floor((longitude + 180.0) / 360.0 * n) ytile = floor((1.0 - log(tan(lat_rad) + (1 / cos(lat_rad))) / pi) / 2.0 * n) return(list(xtile = xtile, ytile = ytile, zoom = zoom)) }
実行結果です。
deg2num(140.113889, 35.613056, 10)
## $xtile
## [1] 910
##
## $ytile
## [1] 403
##
## $zoom
## [1] 10
deg2num(140.113889, 35.613056, 17)
## $xtile
## [1] 116549
##
## $ytile
## [1] 51643
##
## $zoom
## [1] 17
タイル座標が判明したので、次は一枚のタイルがカバーする領域の座標を求めることにします。これには、タイルの一角の座標を取得したあとで、周辺に隣接するタイルの座標を求めるという方法を採用します。まずは先ほどのタイル座標の領域に含まれる緯度経度座標の北東となる座標を取得する処理です。
num2deg <- function(xtile, ytile, zoom) { n = 2.0 ^ zoom longitude = xtile / n * 360.0 - 180.0 lat_rad = atan(sinh(pi * (1 - 2 * ytile / n))) latitude = lat_rad / 0.01745329 return(list(longitude = longitude, latitude = latitude, zoom = zoom)) }
deg2num(140.113889, 35.613056, 10) %>% purrr::pmap( ~ num2deg(..1, ..2, ..3) )
## [[1]]
## [[1]]$longitude
## [1] 139.921875
##
## [[1]]$latitude
## [1] 35.7465174211
##
## [[1]]$zoom
## [1] 10
同様の処理を隣接する全てのタイル座標で実行することで、元のタイル座標の領域を取得できます。これはポリゴンとして扱えるようにしておきます。
tile_bbox <- function(xtile, ytile, zoom) { ytile <- c(ytile, ytile + 1, ytile + 1, ytile, ytile) xtile <- c(xtile, xtile, xtile + 1, xtile + 1, xtile) purrr::map2_dfr( xtile, ytile, ~ num2deg(.x, .y, zoom)) %>% dplyr::select(-zoom) %>% as.matrix() %>% list() %>% sf::st_polygon() %>% sf::st_sfc(crs = 4326) }
deg2num(140.113889, 35.613056, 10) %>% purrr::pmap( ~ tile_bbox(..1, ..2, ..3) )
## [[1]]
## Geometry set for 1 feature
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 139.921875 ymin: 35.4606750714 xmax: 140.2734375 ymax: 35.7465174211
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
確認のめに一度描画してみましょう。北緯 35.613056、東経 140.113889をズームレベル10のタイル座標で示すと次のようになります。
library(mapview) mapview(deg2num(140.113889, 35.613056, 10) %>% purrr::pmap( ~ tile_bbox(..1, ..2, 10) ) %>% .[[1]])
本当かどうか怪しいので、さらに国土地理院のタイル座標確認ページと比較します。作成したポリゴンは、次のURLを表示した時に描画される領域と一致しています。
http://maps.gsi.go.jp/development/tileCoordCheck.html#10/35.6131/140.1139
作成したポリゴンは、次のURLを表示した時に描画される領域と一致しているようで一安心です。
地物を覆う領域のタイル座標
さてここからが本番です。一点の座標ではなく、地物を覆う領域のタイル座標にはどうすれば良いでしょう。答えは単に、対象の地物の矩形からタイルを作成することです。例として、「岡山県」を覆うタイルを描画してみましょう。
library(dplyr) library(sf) library(jpndistrict) df_33 <- jpn_pref(33) b <- df_33 %>% # 「岡山県」の矩形を作成し、その領域を求めます st_bbox() b
## xmin ymin xmax ymax
## 133.266694605 34.298390466 134.413218589 35.352861842
st_bbox()
により、地物の境界となる座標がわかったので、領域を覆うタイル座標を生成します。これにはexpand.grid()
を使っても良いですが、ここではtidyr::crossing()
を用いました。また、元のズームレベルも変数として残すようにしています。
df_tiles <- tidyr::crossing( x = seq( deg2num(b["xmin"], b["ymin"], z = 10)[["xtile"]], deg2num(b["xmax"], b["ymax"], z = 10)[["xtile"]] ), y = seq( deg2num(b["xmin"], b["ymin"], z = 10)[["ytile"]], deg2num(b["xmax"], b["ymax"], z = 10)[["ytile"]] ) ) %>% dplyr::mutate(zoom = 10) dplyr::glimpse(df_tiles)
## Observations: 20
## Variables: 3
## $ x <int> 891, 891, 891, 891, 891, 892, 892, 892, 892, 892, 893, 89...
## $ y <int> 404, 405, 406, 407, 408, 404, 405, 406, 407, 408, 404, 40...
## $ zoom <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 1...
「岡山県」を覆うタイル座標を得ることができました。このタイル座標を再びポリゴンに変換します。
df_tiles <- df_tiles %>% dplyr::mutate(geometry = purrr::pmap(., ~ tile_bbox(..1, ..2, ..3) %>% st_sf(crs = 4326))) %>% tidyr::unnest() %>% st_sf(crs = 4326)
おまけですが、各タイルの重心も変数として保存しておきましょう。
df_tiles <- df_tiles %>% dplyr::mutate(longitude = purrr::pmap_dbl(., ~ st_centroid(..4)[1]), latitude = purrr::pmap_dbl(., ~ st_centroid(..4)[2]))
library(ggplot2) ggplot() + geom_sf(data = df_33) + geom_sf(data = df_tiles, alpha = 0.5) + geom_text(data = df_tiles, aes(longitude, latitude, label = paste(zoom, x, y, sep = "/")), color = "red")
やりました!岡山県をカバーする範囲のタイル座標を重ねることができました。
ですが、右下のタイルなど、一部のタイルには「岡山県」が含まれていません。最後に、地物が含まれないタイルを除外する処理を加えておきましょう。これまでの処理を関数としてまとめ、一度に実行できるようにしておきます。
st_tiles <- function(x, zoom) { b <- st_bbox(x) zoom <- rlang::enquo(zoom) tidyr::crossing( x = seq( deg2num(b["xmin"], b["ymin"], z = !!zoom)[["xtile"]], deg2num(b["xmax"], b["ymax"], z = !!zoom)[["xtile"]] ), y = seq( deg2num(b["xmin"], b["ymin"], z = !!zoom)[["ytile"]], deg2num(b["xmax"], b["ymax"], z = !!zoom)[["ytile"]] ) ) %>% dplyr::mutate(zoom = !! zoom) %>% dplyr::mutate(geometry = purrr::pmap(., ~ tile_bbox(..1, ..2, ..3) %>% st_sf(crs = 4326))) %>% tidyr::unnest() %>% st_sf(crs = 4326) %>% dplyr::mutate(longitude = purrr::pmap_dbl(., ~ st_centroid(..4)[1]), latitude = purrr::pmap_dbl(., ~ st_centroid(..4)[2])) %>% dplyr::mutate(check = purrr::pmap(., ~ st_intersects(..6, !!rlang::eval_tidy(x), sparse = FALSE) %>% rowSums()) > 0) %>% dplyr::filter(check == TRUE) }
ズームレベルも11に変更します。
df_tiles_intersect <- st_tiles(df_33, 11) ggplot() + geom_sf(data = df_33) + geom_sf(data = df_tiles_intersect, alpha = 0.5) + geom_text(data = df_tiles_intersect, aes(longitude, latitude, label = paste(zoom, x, y, sep = "/")), color = "red", size = 1.8)
ちょっと長いコードになってしまったので、もっと簡単な書き方があれば教えてください。
Enjoy!