世界を六角形で表したい

これはなんでしょう。

そう、日本ですね。正しくは日本列島を簡略化し、各都道府県を六角形 (hexagons) で表現した図です。日本列島がおさまってしまうネタ画像が出回るほどに面積の大きな北海道が他の都道府県と同じサイズで小さくなってしまっていたり、現実の都道府県の位置関係を反映していない箇所があったりしますが、パッと見て、日本だなとわかるものではないでしょうか(日本に馴染みのない国籍の人が見てもわからなそう)。

こういった六角形を使った地図はアメリカを例によく見られます。こんなのとか、こんなの。また、最近だとこうした地理空間関係を利用したグラフ表現もよく見られるようになりました。

ikashnitsky.github.io

github.com

北アメリカ大陸は形が整っていて、ずるい...。そう感じてしまいますが、まあ日本でもやってみようという気持ちでやりました。

六角形地図の作り方

まずは六角形にしたい要素、今回は都道府県、の位置関係をXとYの座標で定義します。

library(magrittr)
library(tidyverse)
library(ggthemes)
df.jp.prefs <- tibble::frame_data(
  ~x, ~y, ~id,
  15, 14, "HKD",
  14, 12, "AOM",
  15, 11, "IWT",
  14, 11, "AKT",
  14, 10, "MYG",
  13, 10, "YGT",
  14, 9, "FKS",
  13, 9, "IBR",
  12, 9, "NGT",
  14, 8, "GNM",
  13, 8, "SIT",
  12, 8, "TCG",
  11, 8, "TYM",
  10, 8, "ISK",
  14, 7, "CHB",
  13, 7, "TKY",
  12, 7, "YMN",
  11, 7, "NGN",
  10, 7, "FKI",
  9, 7, "KYT",
  8, 7, "HYO",
  7, 7, "TTR",
  6, 7, "SMN",
  13, 6, "KNG",
  12, 6, "SZO",
  11, 6, "AIC",
  10, 6, "GIF",
  9, 6, "SIG",
  8, 6, "OSK",
  7, 6, "OKY",
  6, 6, "HRS",
  5, 6, "YMG",
  10, 5, "MIE",
  9, 5, "NAR",
  9, 4, "WKY",
  7, 4, "KGW",
  6, 4, "EHM",
  7, 3, "TKS",
  6, 3, "KUC",
  4, 5, "FKO",
  3, 5, "SAG",
  2, 5, "NGS",
  3, 4, "OIT",
  2, 4, "KUM",
  3, 3, "MYZ",
  2, 3, "KGS",
  1, 1, "OKN"
)

ggplot(df.jp.prefs, aes(x = x, y = y, group = id)) +
      geom_point() +
      coord_fixed(ratio = 1) +
      theme_map()

この段階では点ですが、なんとなく日本列島ができています。この点を六角形のポリゴンにしていきます。今回はこの作業にhexamapmakerパッケージを用いました。このパッケージはまだCRANに登録されていないので、利用する際はGitHubからインストールすることになります。

hexamapmakerの関数を使うと、簡単に六角形を作成することができます。早速やってみましょう。この結果が冒頭の図になります。

library(hexamapmaker)
df.jp.prefs <- fix_shape(df.jp.prefs)
df.jp.prefs.hex <- make_polygons(df.jp.prefs)

(p <- ggplot(df.jp.prefs.hex, aes(x, y, group = id)) +
    geom_polygon(fill = "white", colour = "black", show.legend = FALSE) +
    coord_fixed(ratio = 1) +
    theme_map())

sfオブジェクトでもhex

せっかくなのでsfオブジェクトとして扱えるようにしておきましょう。

まずは六角形のポリゴンを作ります。ポリゴンはsfパッケージのst_polygon()を使って簡単に作れます。テキストベースでもRオブジェクトから作っても良いです。

先日、id: yutannihilation さんにpurrrlyrは（メンテされない可能性が大なので）あかん、と釘を刺されましたが、便利なのでpurrrlyrで都道府県ごとのポリゴンを作成します。sfパッケージのst_sfc()により、geom情報を格納したsfcオブジェクトとして扱えるようにしましょう。

library(sf)

make_hex <- function(d) {
  res <- d %>% 
    mutate(geom = sf::st_polygon(list(rbind(c(min(d$x), min(d$y) + 0.577),
                                            c(min(d$x), min(d$y) + 0.577 + 1),
                                            c(mean(d$x), max(d$y)),
                                            c(max(d$x), min(d$y) + 0.577 + 1),
                                            c(max(d$x), min(d$y) + 0.577),
                                            c(mean(d$x), min(d$y)),
                                            c(min(d$x), min(d$y) + 0.577)
    )))) %>% 
    magrittr::use_series(geom)  
  return(res)
}

sfdf.jp.hex <- sf::st_sf(id = sort(df.jp.prefs$id),
               geometry = df.jp.prefs.hex %>% 
                 purrrlyr::slice_rows("id") %>% 
                 purrrlyr::by_slice(make_hex) %>% magrittr::use_series(.out) %>% 
                 sf::st_sfc())
sfdf.jp.hex %>% head()
## Simple feature collection with 6 features and 1 field
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 11 ymin: 5.154 xmax: 29 ymax: 19.924
## epsg (SRID):    NA
## proj4string:    NA
##    id                       geometry
## 1 AIC POLYGON ((21 8.885, 21 9.88...
## 2 AKT POLYGON ((26 16.77, 26 17.7...
## 3 AOM POLYGON ((27 18.347, 27 19....
## 4 CHB POLYGON ((26 10.462, 26 11....
## 5 EHM POLYGON ((11 5.731, 11 6.73...
## 6 FKI POLYGON ((18 10.462, 18 11....

# ラベル用のポリゴン重心データ
sfdf.jp.hex.centroid <- sfdf.jp.hex %>% mutate(x = map_dbl(geometry, ~st_centroid(.x)[[1]]), 
                  y = map_dbl(geometry, ~st_centroid(.x)[[2]]))

ggplot() + 
  geom_sf(data = sfdf.jp.hex) +
  geom_text(data = sfdf.jp.hex.centroid, 
            aes(label = id, 
                          x = x, y = y, alpha = 0.5), family = "IPAexGothic", 
            size = 3,
            show.legend = FALSE)

Enjoy!

地名情報を地図に反映させる

Google MapsやOpenStreetMapに慣れてくると、地名や地物のラベルが地図上に表示されていることに違和感がなくなり、感謝の気持ちが薄れてしまうなと感じる今日この頃です。みなさまいかがお過ごしでしょうか。

今日は、ラベルつきの地図への感謝の気持ちを思い出すべく、国土地理院が試験的に公開している「地名情報のベクトルタイル」を利用して、ラベルつきの地図を作ってみようと思います。

国土地理院提供の地名ベクトルタイル

地名情報のベクトルタイルは次のようなgeojsonで提供されます。

https://cyberjapandata.gsi.go.jp/xyz/experimental_nrpt/{z}/{x}/{y}.geojson

上記のテンプレートURLは地名情報（居住地名）であり、提供されるタイルの種類によってテンプレートURLが異なります。

取得したい範囲のデータは、xyzの座標を変更して取得します。ここで指定するxyzの座標は、タイル座標と呼ばれる多くのウェブ地図で利用されている方式での座標の表現方法です。二次元で表現される緯度経度に縮尺の概念を取り入れているのが特徴となります。

タイル座標は 国土地理院のサービス から調べることもできますが、変換するコードを書いておくと便利です。このページを参考に緯度経度をタイル座標に変換するRコードを書きました。

latlon2tile <- function(lon, lat, z) {
    x = trunc((lon/180 + 1) * 2^z/2)
    y = trunc(((-log(tan((45 + lat/2) * pi/180)) + 
        pi) * 2^z/(2 * pi)))
    return(list(x = x, y = y, zoom = z))
}

この関数を使って緯度経度とズームレベルからタイル座標を取得しましょう。

(xyz <- latlon2tile(133.938, 34.6636, 15))

$x
[1] 28575

$y
[1] 13016

$zoom
[1] 15

試験的に提供されている地名情報ベクトルタイルは住居表示住所以外はズームレベルが15で固定されています。ですので、z引数は固定します。

geojsonファイルの読み込みにはsf::read_sfを使い、sfオブジェクトとして扱えるようにしておきます。

library(sf)
# 居住地名
df.nrpt <- read_sf(paste0("https://cyberjapandata.gsi.go.jp/xyz/experimental_nrpt/", 
    paste(xyz[3], xyz[1], xyz[2], sep = "/"), 
    ".geojson"))
# 公共施設
df.nffpt <- read_sf(paste0("https://cyberjapandata.gsi.go.jp/xyz/experimental_pfpt/", 
    paste(xyz[3], xyz[1], xyz[2], sep = "/"), 
    ".geojson"))

今回指定したタイル座標は、自然地名のラベルが定義されていな箇所だったので、居住地、公共施設2つのラベルファイルを読み込んでおきました。

e-stat小地域ポリゴンデータ

次にベースとなる地図を用意します。市区町村単位のポリゴンであればjpndistrictパッケージが使えますが、今回はより狭い範囲のポリゴンが欲しいです。というわけでe-statの「地図で見る統計」ページから該当する小地域のポリゴンデータをダウンロードしてきます。今回利用したのは「平成２７年国勢調査（小地域） 岡山県 岡山市北区」および「平成２７年国勢調査（小地域） 岡山県 岡山市中区」になります。プログラムで取得したいという際は、以前書いたコードを修正すると大丈夫なはずです。

一部で話題の(?)e-Statが提供する国勢調査小地域のShapefileをRからダウンロードする関数を書きました (gifが高速過ぎる orz) コードは https://t.co/Dd2R58ov94 pic.twitter.com/SRo2H1Fk4W

— Uryu Shinya (@u_ribo) July 14, 2017

library(dplyr)
library(purrr)
# 小地域のデータを読み込みます
pref33 <- read_sf("~/Downloads/A002005212015DDSWC33102") %>% 
    rbind(read_sf("~/Downloads/A002005212015DDSWC33101"))  %>% 
# WGS84の投影座標系を緯度経度に変換しておきます
st_transform(crs = 4326) %>% # 列が多いので不要な列は削除してしまいます
# 小地域の名前が入った列だけを選びます
select(MOJI)

ダウンロードしてきた小地域データを、今回利用する国土地理院タイルの範囲だけに絞ります。ポイントとポリゴンの関係から、st_contains()を使って該当するポリゴンを求めます。

# 地名ポイントが含まれるポリゴンだけを抽出
which.row <- data_frame(contains = pref33$geometry %>% 
    map(~st_contains(.x, df.nrpt$geometry)) %>% 
    flatten() %>% map_lgl(~!identical(., 
    integer(0)))) %>% mutate(row_num = row_number()) %>% 
    filter(contains == TRUE) %>% use_series(row_num)
pref33mod <- pref33 %>% slice(which.row)

マッピング

さてデータが揃ったので地図を作りましょう。

ラベルとして用いる座標はgeometry列に保存されています。これをst_centroid()を使って緯度と経度それぞれ独立した列の値として利用可能にします。

df.nrpt %<>% select(name, kana) %>% mutate(lon = map_dbl(geometry, 
    ~st_centroid(.x)[[1]]), lat = map_dbl(geometry, 
    ~st_centroid(.x)[[2]]))
df.nffpt %<>% select(type, name = pfName) %>% 
    mutate(lon = map_dbl(geometry, ~st_centroid(.x)[[1]]), 
        lat = map_dbl(geometry, ~st_centroid(.x)[[2]]))

こんな感じでどうでしょうか。

最後に遊びごごろでfontawesomeを使って「地図記号」を再現してみます（適切な記号が選べていない感がありますが...）。

library(ggplot2)
library(ggrepel)
ggplot() + geom_sf(data = pref33mod, fill = "white") + 
    geom_text(data = df.nrpt, aes(label = name, 
        x = lon, y = lat, alpha = 0.5), family = "IPAexGothic", 
        size = 3) + geom_sf(data = df.nffpt, 
    aes(color = "red")) + geom_text_repel(data = df.nffpt, 
    aes(label = name, x = lon, y = lat), 
    size = 2, family = "IPAexGothic") + guides(alpha = FALSE, 
    color = FALSE)

library(emojifont)
load.fontawesome()
df.nffpt$label <- fontawesome(c("fa-dot-circle-o", "fa-plus-square",
                                "fa-times-circle", 
                    "fa-times-circle", "fa-university",
                    "fa-envelope", "fa-envelope"))

ggplot() + 
  geom_sf(data = pref33mod, fill = "white") +
  geom_text(data = df.nffpt, aes(label = label, x = lon, y = lat), 
            family = 'fontawesome-webfont') +
  geom_text(data = df.nrpt, aes(label = name, x = lon, y = lat,
                                alpha = 0.5), 
            family = "IPAexGothic", size = 3) + 
  geom_text_repel(data = df.nffpt, aes(label = name, x = lon, y = lat), 
             size = 2,
             family = "IPAexGothic") +
  guides(alpha = FALSE, color = FALSE)

Enjoy!

クレジット

地名情報（居住地名・自然地名・公共施設・住居表示住所）のベクトルタイルを加工して瓜生真也が作成しました。

Rラジオをやってみての感想

先日、年始に立てた目標の一つである「Rラジオ」をやることができた。自分一人では成り立たなくて、ゲストとして id:yutannihilation さんに参加してもらった。感謝しかない。謝謝（収録後に坦々麺を食べに言った🍜）。

uribo.hatenablog.com

口で喋るのは難しい。いや、自分の場合は文章すらまともに書けないのだけど、会話はもっと破滅的になる。何を言っているのか、会話の途中でわからなくなることがしばしばある(ん、これ矛盾してない?みたいな。文であればそれに気が付きにくいというだけかも)。

そんなせいもあり、Rラジオ (公式にはR Radio for the Rest of us)もところどころグダグダしているところがあったり、私が笑いすぎていたりしていて反省も多いのだが、出来としては、個人的には満足している。勉強会やslack、Twitterといったこれまでとは違う形で、RユーザがRの話題を気楽にできる環境としてRラジオが機能していけるということを感じたためだ。

また、これは反省会と称して公開している部分の録音が終わった後にユタニさんと話していた内容なのだが、私からの話題提供もすることができてとても良い。ブログに書くようなほどでもない些細なネタや旬なトピックス、議論したい話題をRユーザとできる。それが楽しい。

録音した自分の声を聞くのは嫌だなと思っていたけれど、聴き始めたらこのラジオ最高や！ってなった(声をもうちょい聴き取りやすくしたい)。また、Twitterでも反響があって嬉しい。がんばるぞい。

rlangradio.org

ちなみに、正式名称である「R Radio for the Rest of us」というのはユタニさん考案のもので、"for the Rest of us" の部分は、その昔あったアップルの標語を再利用しているもの。バリバリにRを使い回す人だけでなく、それ以外の、どんなRユーザでも取りこぼさないぞ！という気概を感じて即採用した。略すとRが三つでR3となる。良い。

最後の方で話しているのだけど、Rラジオは、誰がやってもOK、話したい人どうしで録音してネットにアップロードすれば完成、というものを目指している。

ゲストとの連絡の取り方であったり、運用方法、話す内容（コーナー的な）、質問の募集先など決まっていないことも多々あるのだが、ひとまず船は港を発した。大海原はRラジオという船を運んでくれるだろう。そう期待している。

次回は?

以前何人かに声をかけていたのでその人たちの記憶を思い出させたり、「友達の輪」方式を採用したり、いくつか候補はあるが決まっていない。

みなさんどうですか?

あと、会場としてお世話になったTokyoR Bar&Cafeは良いお店。みなさん是非。

やりました！ #rラジオ #rstats