cucumber flesh

Rを中心としたデータ分析・統計解析らへんの話題をしていくだけ

地物範囲のタイル座標を得る

id:yutannihilationさんが、ggplot2::geom_sf()にOpen Street Map (OSM)のタイルを重ねるという面白い試みをしています。

yutani.rbind.io

OSMでは、画像データとして地図画像を配信していて、次のようなURLで参照されます。

http://[abc].tile.openstreetmap.org/zoom/x/y.png

ドメイン名の他、zoomxディレクトリがあり、ファイルはyという名前で置かれています。ここでzoomは地図のズームレベル、xとyはタイルの座標を示す、ピクセル座標と呼ばれる値です。OSMをはじめ、地図タイルの配信を行なっているサービスは同様にこの「XYZ地図タイル」形式を採用しています。

OSMに限らず、目的の地物を範囲に収めるタイルを用意するには、このピクセル座標を取得する必要があるわけですが、馴染みのある緯度経度からピクセル座標の変換は少し頭を使わなくてはいけません(私だけかもしれませんが)。今回の話題は、すでにyutannihilationさんの記事にも書かれていますが、緯度経度からピクセル座標を取得する処理をおさらいをするというものです(二番煎じです。本当にありがとうございました)。

基礎

緯度経度およびズームレベルからピクセル座標を求める処理は、OSMWikiに書かれている通りです。計算式の詳細は参考にあげるページに詳しいのでそちらをご覧ください。

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を表示した時に描画される領域と一致しています。

f:id:u_ribo:20180610094505p:plain

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

f:id:u_ribo:20180610095646p:plain

やりました!岡山県をカバーする範囲のタイル座標を重ねることができました。

ですが、右下のタイルなど、一部のタイルには「岡山県」が含まれていません。最後に、地物が含まれないタイルを除外する処理を加えておきましょう。これまでの処理を関数としてまとめ、一度に実行できるようにしておきます。

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)

f:id:u_ribo:20180610095629p:plain

ちょっと長いコードになってしまったので、もっと簡単な書き方があれば教えてください。

Enjoy!