読者です 読者をやめる 読者になる 読者になる

まだ厨二病

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

日本の人口密度を可視化する: 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

せっかくなので位置を地図で確認してみましょう。

f:id:u_ribo:20170503104216p:plain

国勢調査は定期的に行われているので、経年変化を見てみるのも楽しそうです。

Enjoy!

広告を非表示にする