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

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

見えないRの関数のソースコードを読む

要約

• lookupパッケージで標準の関数定義ソースコードの出力機能を改善する
  • 総称関数や.C(), .Internal()などの関数で呼び出されるコードも出力
• prettycodeパッケージで関数定義のハイライトを有効にする
• prettycodeはRの起動時に読み込み、lookupは適宜、名前空間を指定してlookup::lookup()で実行、という運用にした

ソースコードの閲覧機能の向上とハイライト機能

最近ちょくちょく、Rのソースコードの読み方が変わっていくんではないかなと思っています。読み方というか出力方法というか。

百聞は一見に如かず。次の画像をご覧ください。この画像には通常のRでのソースコード出力と異なる点が2箇所あります。

この画像はターミナル上で起動したRで、head()のソースコードを表示している場面です。何かお気づきになられるでしょうか。手元にRを実行できる環境がある方は、同じようにソースコードを表示してみると良いかもしれません。

通常のRの実行結果と異なるのは次の点です。

• 関数のソースコードがハイライトされている
• 隠蔽される総称関数のクラスへの参照を行っている

これは通常のRでは行えない機能です。ではどうやっているのかというのが今回の話であり、今後のスタンダードになり得るのではないかと思っているものです。

Rでは、関数名をコンソールに打ち込んでエンターを押すとソースコードが出力されます。ですが、そのままではハイライトはされないし、head()のような総称関数では別途、getS3method("head", "default")を実行しないと隠蔽されたソースコードを読むことができません。

prettycodeおよびlookupはそれぞれ、関数のハイライト、ソースコードの出力機能を改善機能を提供します。prettycodeはlookupをインストールすると自動的にインストールされるので、試したい方はlookupをインストールしましょう。CRANには登録されていないのでGitHub経由で行う必要があります。

# githubinstall::githubinstall('lookup')
library(lookup)

lookupを実行するとprettycodeによるコードのハイライトと関数の出力機能改善が有効になります。

始めの画像のように、lookupはソースコードの出力を行うと、通常のRの挙動とは異なり目的のソースコードをきちんと表示してくれます。Rでは総称関数のようにある関数のメソッドとして定義された関数や.Internal(), .External(),.C(), .Call()によって呼び出される関数のソースを確認するには、逐次的に探索していく必要がありますがlookupを使うとこれらのソースコードの出力が簡単になります。

RStudioでlookupを読み込み、関数のソースコードを出力させると、ソース区画のタブでコードが表示されます。総称関数や.Callなどで呼び出されるコードについても同様にタブで出力されます。いずれもハイライトが有効化されており、オブジェクトがわかりやすくなっています。

運用方法

lookupのコード出力はGitHubに登録されているリポジトリのコードを読みにいくことで実現しています。そのためにGitHub APIを利用しており、認証のない状態では時間あたりの利用回数に制限があります。ですので、トークンを発行し設定して規制を緩和させることが推奨されています。

このような理由もあることから、私はlookupの読み込みは行わず、ソースを見たいときにlookup::lookup("function")のようにして必要な時にだけ有効にすることにしました。一方で通常の関数ソースコードの出力はハイライトさせたいのでprettycodeは.Rprofileに記述し、起動時に読み込むようにしました。

これまで隠蔽されたコードを読むには、面倒な手続きやいくつかの操作が必要でしたが、lookupを使うことで手間が大きく省けそうです。

Enjoy!