cucumber flesh

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

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

mapping sf purrr dplyr ggrepel emojifont

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の「地図で見る統計」ページから該当する小地域のポリゴンデータをダウンロードしてきます。今回利用したのは「平成２７年国勢調査（小地域） 岡山県 岡山市北区」および「平成２７年国勢調査（小地域） 岡山県 岡山市中区」になります。プログラムで取得したいという際は、以前書いたコードを修正すると大丈夫なはずです。

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

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

f:id:u_ribo:20171012082602p:plain

最後に遊びごごろで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)

f:id:u_ribo:20171012082844p:plain

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!

クレジット

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