日本の人口密度を可視化する: population lines

少し前(4月下旬ごろ?)に、redditで人口密度の高さを表現した地図が話題になりました。

www.reddit.com

この地図は、James Cheshire博士 (@spatialanalysis)が2014年に投稿した “Population Lines Print” が元となっていて、再現性のあるRコード、ヨーロッパに焦点を当てた地図が描かれた（Henrik Lindberg @hnrklndbrg ）ことで話題が広がっています（という印象。今週のR Weeklyでもいくつかの記事が掲載されました。

日本の人口密度の情報を世界地図から伺うことはできますが、スケールダウンしたものがあった方がわかりやすいです。…というわけで作成してみました。

f:id:u_ribo:20170503102859p:plain

政府統計の総合窓口 e-Statでも「人口」に関するデータはたくさん公開されていますが、都道府県や市区町村より細かな人口データとして、「平成22年国勢調査の地域メッシュ統計 1kmメッシュ男女別人口総数及び世帯総数 (2010/10/01)」を用いることにしました。手動でのダウンロードだと気が遠くなりそうな作業だったのでRSeleniumを使って自動化しました。

以下、Rコードです。メッシュを扱うためjpmeshパッケージを使います。集計方法や可視化のコードはHenrikさんのものを応用しています。

library(magrittr)
library(tidyverse)
library(jpmesh)

1次メッシュ (5kmメッシュ)に分割されたテキストファイル (カンマ区切りなので注意)を読み込みます。同一形式の複数ファイル読み込みにはpurrr::map()を使うと便利です。データフレームとしてまとめたら、purrr::set_names()により変数名を変更しておきます。

d <- list.files(pattern = ".txt$", full.names = TRUE) %>% 
    map(~read_csv(file = ., locale = locale(encoding = "cp932"), 
        skip = 1)) %>% reduce(bind_rows) %>% 
    set_names(c("meshcode8", "population", 
        "man", "woman", "total_setai"))

このようなデータです。1kmメッシュコード、人口総数、男性の人口総数、女性の人口総数、世帯総数が含まれます。

meshcode8	population	man	woman	total_setai
36225735	54	34	20	19
36225737	2	2	0	2
36225738	113	57	56	46
36225745	437	235	202	191
36225758	10	5	5	4
36225759	6	6	0	6

メッシュから緯度経度を割り当てたデータを用意します。これは可視化の際の座標として用います。

d.mesh <- d$meshcode8[1:nrow(d)] %>% purrr::map(jpmesh::meshcode_to_latlon) %>% 
    reduce(bind_rows) %>% dplyr::select(lat = lat_center, 
    lng = long_center) %>% mutate(meshcode8 = d$meshcode8[1:nrow(d)])

d %<>% left_join(d.mesh, by = "meshcode8") %>% 
    group_by(lat = round(lat, 1), lng = round(lng, 
        1)) %>% summarize(value = sum(population, 
    na.rm = TRUE)) %>% ungroup() %>% complete(lat, 
    lng)

最初に表示した図を描画します。

d %>% ggplot(aes(lng, lat + 5 * (value/max(value, 
    na.rm = TRUE)))) + geom_line(size = 0.4, 
    alpha = 0.65, color = "#5A3E37", aes(group = lat), 
    na.rm = TRUE) + ggthemes::theme_map() + 
    coord_equal(0.9)

特に人口密度の高い市区町村はどこでしょうか。以前書いた、緯度経度から該当する市区町村を判定するコードを使って確かめます。

qiita.com

d.mesh.bind <- d %>% arrange(desc(value)) %>% 
    head(10) %>% left_join(d.mesh %>% mutate(lat = round(lat, 
    1), lng = round(lng, 1)) %>% distinct(lat, 
    lng, .keep_all = TRUE))
d.density.top <- d.mesh.bind %>% slice_rows("meshcode8") %>% 
    by_slice(~find_city(dfs.jp, lon = .$lng, 
        lat = .$lat)) %>% unnest() %>% 
    left_join(d.mesh.bind, by = "meshcode8") %>% 
    select(city_name, value, lat, lng)

city_name	value
東京都新宿区	990465
東京都大田区	906408
東京都墨田区	752090
東京都杉並区	739990
大阪府大阪市北区	719928
埼玉県川口市	633078
東京都江戸川区	570923
神奈川県川崎市高津区	570452
大阪府大阪市住吉区	515699
大阪府東大阪市	489196