cucumber flesh

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

geofacetで日本のデータが利用可能になりました

もう数ヶ月も前の話ですが、CRANに登録されているgeofacetパッケージが更新されて、バージョンが0.1.9となりました。geofacetパッケージは、ggplot2ベースの作図パッケージの一種ですが、グラフを地理的な空間関係と対応付けて描画するという特徴があります。例えば、次のような図を見たことがある人がいるのではないでしょうか。

f:id:u_ribo:20180622064350p:plain

この図はアメリカでの2000年から2016年にかけての失業率を州に分けて表示しています。ただ図を並べるだけでなく、各州の空間関係を考慮した配置をすることで各州の関係がわかりやすく可視化されていると思います。この図はgeofacetのサンプルコードとして掲載されているもので、下記を実行することで描画されます。

library(ggplot2)
library(geofacet)

ggplot(state_unemp, aes(year, rate)) +
  geom_line() +
  facet_geo(~ state, grid = "us_state_grid2") +
  scale_x_continuous(labels = function(x) paste0("'", substr(x, 3, 4))) +
  ylab("Unemployment Rate (%)")

geofacetには、こうした図を書くために、いくつかの国や地域に対応していますが、バージョン0.1.9で日本にも対応するようになりました!追加した当人として、geofacetで日本のデータを利用する例を紹介します。アメリカに比べると少し不恰好で、現実の大きさや配置を再現できていない部分もあるかと思いますが、ぜひお使いください!

github.com

続きを読む

ggplot2のsizeが意味するもの

付け焼き刃な知識故に、間違いございましたら指摘頂けますと幸いです。

ggplot2を使っていると、数種類の“size”を設定するタイミングがある。例えば、geom_point()で散布図を描画する際のポイントの大きさを指定するためのsize、または連続した線を描画するのに用いられるgeom_line()でのsize、あるいは文字を重ねるためのgeom_text()annotate()の中でのsizeである。また、少し異なる引数名であるがテーマ指定の関数theme_*()にも“base_size”という名の引数が用意されている。ここで肝心なことは、テキストとその他では、同じsizeであっても単位が異なるということだ。このことは次のStackOverflowの投稿にまとめられている。

stackoverflow.com

まとめて書いてしまうとこの通り。

  • レイヤーとして重ねるgeom_point()geom_line()などの要素のサイズはmm単位(ポイントの大きさ、線の太さなど)
  • 軸やタイトルなど、フォントを表示するテキスト系のサイズはpoint単位(1pointsはおよそ0.352778mm)
  • geom_point()のポイントの輪郭の太さを制御するstrokeの単位もpoint(ggplot2::.stroke

以下、この記事ではそのことを確かめてみたいと思う。

適当なggオブジェクトを用意し、それに手を加えて比較していくことにする。

library(ggplot2)
# 日本語の表示のためにフォントを指定
theme_set(theme_gray(base_family = "IPAexGothic"))

p <- 
  ggplot(mtcars, aes(wt, mpg))

まずはgeom_point()geom_line()でのsizeである。ここではsize = 10を指定した(既定では1.5)。

p_out <- 
  p + 
  geom_point(size = 10)

この結果をggsave()で出力した図を表示する。ポイントのサイズが10mmにしたので、検証のために出力時のサイズはその10倍となる10cmとした。

ggsave(filename = "out.png", 
       plot = p_out, 
       width = 10, 
       height = 10, 
       units = "cm")

f:id:u_ribo:20180611230853p:plain

確かにmmらしい。もっとも、これはソースコードにも示されている。

Size examples Should be specified with a numerical value (in millimetres), or from a variable source

github.com

次に、軸やタイトルテキストのsizeである。軸のラベルを除いて、size = 10を指定した。

p2 <- 
  p_out +
  # 適当な場所に文字「A」を表示
  annotate("text", x = 2, y = 15, size = 10, label = "A", color = "red") +
  ggtitle("文字の大きさはpointで指定")

p2 + 
  theme(title = element_text(size = 10),
         axis.title = element_text(size = 10),
         # グラフの両軸に用いられる数値の大きさ
         axis.text = element_text(size = 5))

f:id:u_ribo:20180611231008p:plain

フォントを扱っていて、なおかつ同じsize = 10を指定していてもannotate()が表示するフォントが大きい。これはannotate()がレイヤーの要素として評価されるためで、geom_point()のようにsizeはmm単位での指定となるためだ。

もし、タイトルや軸のサイズをmm単位にしたい時はpointをmmへ変換する必要がある。ggplot2::.ptはggplot2の内部で用いられている変数で、ドキュメントもある

p2 +
  theme(axis.title = element_text(size = ggplot2::.pt * 10))

f:id:u_ribo:20180611231028p:plain

最後にtheme_*()の引数“base_size”だが、こちらも単位はpointである。すなわち次の2つのコードによって出力される図は等しい。

ggplot2が生成するような描画オブジェクトの差分を確認するにはvdiffrが便利である。ここではwidget_diff()を用いて差分の検証を行うが、GitHubのように2つの画像を表示し、スライドさせながら差分を確認するwidget_slide()やテスト時に役立つ関数も備わっている。

qiita.com

# plotとaxisのラベルの大きさは、テーマのbase_sizeに対してI()で指定した値をかけた大きさとなる
p_default <- 
  p2 + 
  theme_gray(base_size = 11,
             base_family = "IPAexGothic")

p_custom <- 
  p2 +
  theme_replace(
    plot.title = element_text(size = 11 * I(1.2)),
    axis.title = element_text(size = 11 * I(0.8)))
# 2つのオブジェクトで差分がない場合は白い画面が表示される
vdiffr::widget_diff(
  p_default,
  p_custom
)

# 変更してみる
# フォントサイズが変更されたことで全体の配置も微妙にずれる
vdiffr::widget_diff(
  p_default,
  p_custom +
    theme(plot.title = element_text(size = 20))
)

f:id:u_ribo:20180611231716p:plain

学術誌に投稿する図では、軸のテキストサイズやらについて細かく指定しされているものもある。怒られないように、少しずつggplot2に慣れていこうと思う。 Enjoy!

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

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!