2015年11月1日日曜日

Rを使って郵便番号界の地図を作る(統計局のShape+日本郵便のCSV)

郵便番号がどんな感じで割り当てられているか、地図上で確認したかったんです。

GPSで取った地点情報がプライバシーに配慮云々で郵便番号になっていた、なんてことがありまして、その粒度がどんなものか確認したかったってのが発端です。

検索してみても無償のものは意外に見つからないもので、国際航業さんが売っていたりします↓

PAREA-Zip郵便番号界 | 国際航業株式会社

ノースリーブな(=無い袖は振れない)わたくしとしては、自分で作るしかないわけで。

郵便番号のリストは↓ここにありました。

郵便番号データダウンロード - 日本郵便

住所文字列と郵便番号が入っていればどれでもいいんですが、とりあえず、「住所の郵便番号(CSV形式)」の「読み仮名データの促音・拗音を小書きで表記するもの」というのをダウンロードしました。

次は地図です。こういうときに使う境界データとしては、統計局のe-StatでShapeファイルが公開されているので、これが便利。入手の方法は下記を参照。

e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ

今回は、我らがサラリーマンの聖地 新橋を有する港区を対象にしてみました。

とりあえず地図をプロットするところまでやってみましょう。

library(maptools)     # GIS系のパッケージ
# 座標系の情報を持つオブジェクトを作る
# longlatは緯度経度指定、WGS84は世界測地系であることを表す
pj <- CRS("+proj=longlat +datum=WGS84")

# シェープファイルの読み込み
map <- readShapePoly("h22ka13103.shp", proj4string=pj)

# とりあえずプロット
plot(map)

shapeをプロットしただけのもの

この地図上に追加で郵便番号をプロットすればよさそうですね。

読み込んだシェープファイルにはポリゴンの境界情報だけでなく、それぞれのポリゴンの属性の情報も含まれています。

こうやれば確認できます

head(map@data)

       AREA  (略)  KEN_NAME SITYO_NAME GST_NAME  (略)         MOJI
0 729237.30          東京都       <NA>     港区       元赤坂2丁目
1 148246.50          東京都       <NA>     港区       北青山1丁目
2  89514.52          東京都       <NA>     港区       元赤坂1丁目
3  93129.48          東京都       <NA>     港区         赤坂3丁目
4 151487.30          東京都       <NA>     港区         赤坂4丁目
5 199895.30          東京都       <NA>     港区       北青山2丁目

「MOJI」っていう列に住所が入っているようです。これだけにアクセスしたい場合は、「map$MOJI」と打てばOK。厳密には「map@data$MOJI」なんでしょうけど、省略できるようです。

これを郵便局のCSVファイルとマッチさせればよさそう。

# 郵便番号データの読み込み
post_data <- read.csv("13TOKYO.CSV", head=F)
post_data <- post_data[post_data$V8=="港区", ] # 港区だけ取り出し

# ○丁目の部分をカットして、MOJI2の列に格納
map$MOJI2 <- sub("[0-9].*", "", map$MOJI)

# MOJI2と郵便番号データの文字列がマッチしたら、
# 郵便番号を取り出して設定する
for ( i in 1:nrow(map@data) ) {
  # 複数マッチしたら、安直に1つ目を採用
  matched_index    <- grep(map$MOJI2[i], post_data$V9)[1]
  map$POST[i] <- post_data$V3[matched_index]
}

plot(map)
text( coordinates(map), as.character(map$POST), cex=0.5)

郵便番号をtext関数で追加プロットしたもの

一応、最低限のところまではできました。でも、文字が見にくいとか、地名も欲しいとか、色塗ったほうがキャッチーだよね(?)とか、ということで付け足してみたのが以下。

# ハイフンを入れる
post <- paste ( substr(map$POST, 1, 3), substr(map$POST, 4, 7), sep="-" )

# 地名&郵便番号の文字列を作る
labs <- paste(map$MOJI, post, sep="\n")

# 塗りわけのパレット
library(RColorBrewer) # カラーパレット
cols <- brewer.pal(12, "Set3")
cols <- rep(cols, 10) # 12色じゃ足りないので繰り返し

pdf("郵便番号界マップ(港区).pdf", family="Japan1GothicBBB", width=16.54, height=11.69)

# 郵便番号を因子型にすると、番号ごとに塗りわけできる
plot(map, col=cols[factor(map$POST)], border="gray", lwd=0.1)
text( coordinates(map), labs, cex=0.3)

dev.off()

上記のサンプルコードで作ったのが↓これです。

郵便番号界マップ(港区)(PDFファイル)


↓PDFファイルのスクリーンショット
郵便番号に、色と地名を追加したもの

なかなかいい感じになったんじゃないでしょうか。郵便番号ごとに塗り分けるために、因子型にするなんてのは、ちょっとしたtipsです。

が、ツッコミが入る前に白状しておかねばなりますまい。この郵便番号地図、正しくないところがあります。

住所のマッチングが結構甘いです。grepでマッチしたものが複数ある場合、単純に1つ目のものを採用しています。

例えば、↓こういうパターンに対してはうまく機能しています。

"1070052", "赤坂(次のビルを除く)"
"1076090", "赤坂赤坂アークヒルズ・アーク森ビル(地階・階層不明)"
"1076001", "赤坂赤坂アークヒルズ・アーク森ビル(1階)"
・・・

が、↓こんなやつだとうまくいきませんね。

"1050014", "芝(1~3丁目)"
"1080014", "芝(4、5丁目)"

芝1丁目、芝2丁目、・・・、芝5丁目を全部別の行にしてくれていたとしたら、まだプログラム側でなんとかしようなんて気も起きますが、「1~3」とか「4、5」とか、機械的に処理することを全く考慮していないデータなんですもの。今後の日本郵便のがんばりに期待したいところですね。(上から目線)

ということで、こんな間違いだらけのマップじゃ使えないよってかたは、国際航業さんから買ってください。




0 件のコメント:

コメントを投稿