Geocodingで自前のGISデータを作る
住所や郵便番号などの情報に対して、その地理座標（緯度・経度など）を付加する作業をGeocodingといいます。一般的には、地理座標データベースの中から、対応する住所を見つけ出す作業、つまり文字列のマッチングを行います。
文字列に関する作業ですから、各言語特有の処理が必要になります。例えば、英語の場合は略称の処理（St -> Street, W -> Westなど）がキーになります。本稿では、日本語のGeocodingを紹介します。
対象データ
まずはGeocodingの対象となるデータとして、今回は賃貸の家賃データを用います。これに地理座標を与えて、賃料の分布を視覚的に表現する（家賃マップの作成！）ことを試みたいと思います。対象として、東急目黒線沿線を選びました。都心に出やすく住みやすい人気のエリアのひとつです。
都営三田線に乗り入れており、三田・日吉キャンパスをつなげる慶應大学生のバイタルラインでもあります。実際的な理由としては、東京都と神奈川県にまたがりつつも、駅数が13と少なく例として扱いやすいことから選びました。
東急目黒線沿線の賃貸物件の情報を、住宅情報サイトSUUMOより取得しました。HTMLソースを処理してデータを抽出する、Webscrapingと呼ばれる方法を使うのですが、これについては別に機会があれば紹介したいと思います。
以下の作業は統計処理ソフトRを用います。まず、取得したデータを参照してみますと、このよう形式になっています。
> head(x[, c('rent', 'address', 'type', 'madori', 'year.direc')])
rent address type madori year.direc
1 10.6万円 神奈川県川崎市中原区下小田中５ アパート 2DK 築11年/南
2 10.6万円 神奈川県川崎市中原区下小田中５ アパート 2DK 築11年/南
3 10.6万円 神奈川県川崎市中原区上丸子八幡町 マンション 2DK 築25年/南西
4 10.6万円 神奈川県川崎市中原区今井西町 アパート 1LDK 築2年/西
5 10.6万円 神奈川県川崎市中原区今井西町 アパート 1LDK 築2年/西
6 10.6万円 神奈川県川崎市中原区今井西町 アパート 1LDK 築2年/西
住所はaddress変数に格納されており、どうやら町丁名レベルまで（〇〇丁目など）が記載されていて、街区（△△番地など）は略されていることがわかります。この文字列データに対して、対応する緯度・経度の情報をマッチングしていくのがGeocodingの目的です。
地理座標データベース
日本の地理座標データベースは、国土交通省国土政策局国土情報課のページから、ダウンロードすることができます。左記のリンクで利用約款を読んでから、「位置情報ダウンロードサービスへ」のリンクへ進みます。
次のようなページに行き着きます。ここでは「都道府県単位」を選びます。
都道府県選択です。今回の対象には東京都と神奈川県の住所が含まれるので、その２つを選択します。
ダウンロードするファイルを選択します。各都道府県ごとに2種類のファイルが用意されています。「街区」レベルのファイルには番地までのデータが、「大字・町丁目」レベルのファイルには〇〇丁目までのデータが格納されています。
「街区」レベルのほうが「大字・町丁目」レベルよりも細かく、ファイルサイズも非常に大きくなっています。今回の対象となる住所には番地までの情報は含まれていないので、「大字・町丁目」レベルを選択します。
このあと、利用約款に同意するとダウンロードリンクが表示されます。
ダウンロードしたZipファイルを解凍するといくつかのファイルがありますが、CSVファイルを使います。これをRで読み込んで、データを参照した結果が以下です。
> head(y)
都道府県コード 都道府県名 市区町村コード 市区町村名 大字町丁目コード
1 13 東京都 13101 千代田区 1.3101e+11
2 13 東京都 13101 千代田区 1.3101e+11
3 13 東京都 13101 千代田区 1.3101e+11
4 13 東京都 13101 千代田区 1.3101e+11
5 13 東京都 13101 千代田区 1.3101e+11
6 13 東京都 13101 千代田区 1.3101e+11
大字町丁目名 緯度 経度 原典資料コード 大字.字.丁目区分コード
1 一ツ橋一丁目 35.69163 139.7567 1 3
2 一ツ橋二丁目 35.69295 139.7573 1 3
3 一番町 35.68772 139.7397 1 1
4 永田町一丁目 35.67633 139.7457 1 3
5 永田町二丁目 35.67571 139.7405 1 3
6 猿楽町一丁目 35.69847 139.7599 1 3
「都道府県名」「市区町村名」「大字町丁目名」をくっつけて住所変数をつくり、また緯度・経度の変数名を変更しておきます。
> y$address <- paste(y[,2], y[,4], y[, 6], sep="")
> names(y)[7:8] <- c('lat', 'lon')
> head(y[, c('address', 'lat', 'lon')])
address lat lon
1 東京都千代田区一ツ橋一丁目 35.69163 139.7567
2 東京都千代田区一ツ橋二丁目 35.69295 139.7573
3 東京都千代田区一番町 35.68772 139.7397
4 東京都千代田区永田町一丁目 35.67633 139.7457
5 東京都千代田区永田町二丁目 35.67571 139.7405
6 東京都千代田区猿楽町一丁目 35.69847 139.7599
本題：住所文字列のマッチング
ここからが本題。ケースバイケースの手探りの作業です。まず、このままの状態でどれくらいの住所が一致しているかを調べてみます。xが賃貸データ、yが緯度・経度データです。
> table(x$address %in% y$address)
FALSE TRUE
24682 4718
上のコードで、xのうちのどれだけのケースがyの中に現れるかを数えています。残念ながら、このままでは6分の1弱しかマッチしません。理由は明白で、xのデータでは町丁名の表記が全角数字（例：神奈川県川崎市中原区下小田中5）となっているのに対して、yのデータでは漢数字になっているからです（例：東京都千代田区一ツ橋一丁目）。
この形式を一致させる必要があるのですが、数値と漢数字の変換はそれなりに面倒です。まず具体的にxに現れる漢数字の種類を調べてみましょう。
> p <- regexpr('[0-9]+$', x$address)
> suji <- substring(x$address, p)
> suji[p < 0] <- ''
> table(suji)
suji
5 3 4 1 2 6 7 8
4728 2555 4624 2569 6133 6420 1726 583 62
コード中のregexpr()は正規表現を用いて文字列パターンを探す関数で、’[0-9]+$’ は「末尾に現れる0から9の最長列」を示します。結果によると、データ内に現れる数字はいずれもひと桁ですから、これらを〇〇丁目と漢数字に置き換えれば今回の場合は用が足ります。
kanji <- c('一', '二', '三', '四', '五', '六', '七', '八', '九')
nums <- c('1', '2', '3', '4', '5', '6', '7', '8', '9')
for (i in 1:length(kanji)) {
x$address <- sub(nums[i], paste(kanji[i], '丁目', sep=''), x$address)
}
改めて住所の一致具合を調べてみると、下記のように随分改善されて不一致はわずかに11件となりました。残り少ないのでその内訳を調べてしまいましょう。すると、不一致の住所の種類は４件だけです（同じ住所に複数の物件が入っているため）。
> table(x$address %in% y$address)
FALSE TRUE
11 29389
> table(x$address[!(x$address %in% y$address)])
神奈川県川崎市中原区小杉
5
神奈川県川崎市幸区小倉
3
神奈川県川崎市幸区鹿島田三丁目
1
神奈川県川崎市高津区子母口子母口富士見台
2
なぜこれらがyのデータに現れないのかをひとつずつ調べていきます。
まず、中原区小杉については、ぴったり一致するものがないようです。これが小杉町のことだとして、何丁目かを指定する必要があるようです。今回のケースでは、これは不確定ですから欠損（マッチ不能）として扱うことにします。
> y$address[grep('神奈川県川崎市中原区小杉', y$address)]
[1] "神奈川県川崎市中原区小杉御殿町一丁目"
[2] "神奈川県川崎市中原区小杉御殿町二丁目"
[3] "神奈川県川崎市中原区小杉陣屋町一丁目"
[4] "神奈川県川崎市中原区小杉陣屋町二丁目"
[5] "神奈川県川崎市中原区小杉町一丁目"
[6] "神奈川県川崎市中原区小杉町二丁目"
[7] "神奈川県川崎市中原区小杉町三丁目"
幸区小倉も同様に、何丁目かを指定する必要があるようです。一から五丁目の平均値をとる手もありますが、今回は欠損扱いにすることにします。
> y$address[grep('神奈川県川崎市幸区小倉', y$address)]
[1] "神奈川県川崎市幸区小倉五丁目" "神奈川県川崎市幸区小倉二丁目"
[3] "神奈川県川崎市幸区小倉四丁目" "神奈川県川崎市幸区小倉一丁目"
[5] "神奈川県川崎市幸区小倉三丁目"
“神奈川県川崎市高津区子母口子母口富士見台”という不思議な住所ですが、これはどうやら入力ミスのようです。「子母口子母口」を「子母口」に修正します。
今回は綺麗に整えられた情報サイトからのデータなのでこれはレアケースとなっていますが、多くのGeocodingの実践では、こういった入力時のエラーや形式の非統一への対処が主たるタスクになります。
> y$address[grep('神奈川県川崎市高津区子母口', y$address)]
[1] "神奈川県川崎市高津区子母口"
[2] "神奈川県川崎市高津区子母口富士見台"
> x$address <- sub('子母口子母口', '子母口', x$address)
“神奈川県川崎市幸区鹿島田三丁目” に関しては、データベース上では「幸区鹿島田」として記載されているようです。もしかして三丁目しかないのかと思いましたが、調べてみたらそうではないらしいです。理由はわかりませんが、「三丁目」を削除してマッチングすることにします。
> y$address[grep('神奈川県川崎市幸区鹿島田', y$address)]
[1] "神奈川県川崎市幸区鹿島田"
> x$address <- sub('神奈川県川崎市幸区鹿島田三丁目', '神奈川県川崎市幸区鹿島田', x$address)
結果、8件の欠損が残りましたが、まずまずの割合だと思います。address変数でxとyを結合すればgeocoding完成です。
> table(x$address[!(x$address %in% y$address)])
神奈川県川崎市中原区小杉 神奈川県川崎市幸区小倉
5 3
> x <- merge(x, y)
家賃マップを作成する
早速作ったデータをもとに、目黒線沿線、目黒・日吉間の家賃マップを描いてみましょう（コードは下にまとめて載せます）。比較のため、「間取り１K、アパートまたはマンション、築15年以下」のものに限定しました（5261件）。また、参考のために目黒線の各駅を重ねてプロットしました（駅の地理座標はWikipediaを参照）。
この図では、プロットの色で賃料の高さを表しています。赤に近いほど高く、青に近いほど安い物件です。当然ですが、都心に近づくほど高いですね。
真ん中に物件のない線状の空間がありますが、多摩川でしょう。これより南が神奈川県川崎市で、7万円台以下の物件がぐっと増えます。人気の武蔵小杉や元住吉がこのあたりです。
この図は、データから賃金を局所平均によってスムーズにした結果を表示しています。また、周辺がごちゃごちゃしていないので、駅名もプロットしてあります。特に神奈川県に入ってからは、駅から離れるほど賃料が安くなる傾向が出ています。都内ではそうでもないですが、このあたりは目黒線の他にもたくさんの路線があるためだと思われます。
奥沢の少し上にぽこっと賃料の高くなるエリアがありますが、自由が丘がこのあたりです。このから分岐して目黒線の北側に東横線が伸びており、そのためか南側に比べて賃料が高くなっています。図の左上・左下に赤くなっているエリアがありますが、このあたりにはほぼデータがないので、これらは信頼性のない推計値です。
結び
本稿では、Geocodingの基礎的な方法を紹介しました。日本語の住所の文字列に対応する地理座標情報をマッチングすることが主眼ですから、応用の幅は広いと思います。国土交通省のデータは無料で公開されているので、何かの際には活用されてみてはいかがでしょうか。
＜グラフ用のメインコード＞
x$rent <- as.numeric(gsub('[万円]', '', x$rent))
x$year <- as.numeric(gsub('[^0-9]', '', x$year.direc))
x$menseki <- as.numeric(gsub('m2', '', x$menseki))
z <- read.table('station.txt', header=F)
flg <- x$madori %in% c('1K') &
x$type %in% c('アパート', 'マンション') &
x$year <= 15
w <- x[which(flg),]
library(gplots)
library(plotrix)
rng <- range(w$rent)
p <- 0.05
pt <- quantile(w$rent, c(p, 1-p))
lev <- c(floor(rng[1]), seq(floor(pt[1]), ceiling(pt[2]), by=0.5), ceiling(rng[2]))
nlev <- length(lev)
col <- rep('', nrow(w))
col.set <- rich.colors(nlev-1)
col.labels <- rep('', length=length(col.set))
for (i in 2:nlev) {
col[w$rent > lev[i-1] & w$rent <= lev[i]] <- col.set[i-1]
if (i==2) {
col.labels[i-1] <- paste(' ', '~', lev[i])
} else if (i==nlev) {
col.labels[i-1] <- paste(lev[i-1], '~', ' ')
} else {
col.labels[i-1] <- paste(lev[i-1], '~', lev[i])
}
}
png('plot.png', width=600, height=480)
par(mar=c(5, 4, 4, 6))
plot(w$lon, w$lat, col=col, pch=16, xlab='longitude', ylab='latitude', cex=1.2)
points(z$V4, z$V3, type='o', pch=15, lwd=4, col='azure2', cex=1.8)
points(z$V4, z$V3, type='o', pch=15, lwd=3, cex=1.5)
color.legend(par()$usr[2] + 0.003, par()$usr[3],
par()$usr[2] + 0.013, par()$usr[4],
col.labels, col.set, gradient='y', align='rb', cex=0.9)
dev.off()
source('./filled.contour3.R')
gd <- 100
fit <- loess(rent ~ lon + lat, data=w, span=0.5)
lon <- seq(min(w$lon), max(w$lon), length=gd)
lat <- seq(min(w$lat), max(w$lat), length=gd)
D <- expand.grid(lon=lon, lat=lat)
pre <- matrix(predict(fit, D), nrow=gd, ncol=gd)
png('cont.png', width=600, height=480)
filled.contour3(lon, lat, pre, color=rich.colors,
xlab='longitude', ylab='latitude')
points(z$V4, z$V3, type='c', pch=15, lwd=4, col='azure2', cex=1.8)
points(z$V4, z$V3, type='c', pch=15, lwd=3, cex=1.5)
text(z$V4, z$V3, z$V1)
dev.off()
filled.contour3 <-
function (x = seq(0, 1, length.out = nrow(z)),
y = seq(0, 1, length.out = ncol(z)), z, xlim = range(x, finite = TRUE),
ylim = range(y, finite = TRUE), zlim = range(z, finite = TRUE),
levels = pretty(zlim, nlevels), nlevels = 20, color.palette = cm.colors,
col = color.palette(length(levels) - 1), plot.title, plot.axes,
key.title, key.axes, asp = NA, xaxs = "i", yaxs = "i", las = 1,
axes = TRUE, frame.plot = axes,mar, ...)
{
# modification by Ian Taylor of the filled.contour function
# to remove the key and facilitate overplotting with contour()
# further modified by Carey McGilliard and Bridget Ferris
# to allow multiple plots on one page
if (missing(z)) {
if (!missing(x)) {
if (is.list(x)) {
z <- x$z
y <- x$y
x <- x$x
}
else {
z <- x
x <- seq.int(0, 1, length.out = nrow(z))
}
}
else stop("no 'z' matrix specified")
}
else if (is.list(x)) {
y <- x$y
x <- x$x
}
if (any(diff(x) <= 0) || any(diff(y) <= 0))
stop("increasing 'x' and 'y' values expected")
# mar.orig <- (par.orig <- par(c("mar", "las", "mfrow")))$mar
# on.exit(par(par.orig))
# w <- (3 + mar.orig[2]) * par("csi") * 2.54
# par(las = las)
# mar <- mar.orig
plot.new()
# par(mar=mar)
plot.window(xlim, ylim, "", xaxs = xaxs, yaxs = yaxs, asp = asp)
if (!is.matrix(z) || nrow(z) <= 1 || ncol(z) <= 1)
stop("no proper 'z' matrix specified")
if (!is.double(z))
storage.mode(z) <- "double"
.filled.contour(as.double(x), as.double(y), z, as.double(levels),
col = col)
if (missing(plot.axes)) {
if (axes) {
title(main = "", xlab = "", ylab = "")
Axis(x, side = 1)
Axis(y, side = 2)
}
}
else plot.axes
if (frame.plot)
box()
if (missing(plot.title))
title(...)
else plot.title
invisible()
}
filled.legend <-
function (x = seq(0, 1, length.out = nrow(z)), y = seq(0, 1,
length.out = ncol(z)), z, xlim = range(x, finite = TRUE),
ylim = range(y, finite = TRUE), zlim = range(z, finite = TRUE),
levels = pretty(zlim, nlevels), nlevels = 20, color.palette = cm.colors,
col = color.palette(length(levels) - 1), plot.title, plot.axes,
key.title, key.axes, asp = NA, xaxs = "i", yaxs = "i", las = 1,
axes = TRUE, frame.plot = axes, ...)
{
# modification of filled.contour by Carey McGilliard and Bridget Ferris
# designed to just plot the legend
if (missing(z)) {
if (!missing(x)) {
if (is.list(x)) {
z <- x$z
y <- x$y
x <- x$x
}
else {
z <- x
x <- seq.int(0, 1, length.out = nrow(z))
}
}
else stop("no 'z' matrix specified")
}
else if (is.list(x)) {
y <- x$y
x <- x$x
}
if (any(diff(x) <= 0) || any(diff(y) <= 0))
stop("increasing 'x' and 'y' values expected")
# mar.orig <- (par.orig <- par(c("mar", "las", "mfrow")))$mar
# on.exit(par(par.orig))
# w <- (3 + mar.orig[2L]) * par("csi") * 2.54
#layout(matrix(c(2, 1), ncol = 2L), widths = c(1, lcm(w)))
# par(las = las)
# mar <- mar.orig
# mar[4L] <- mar[2L]
# mar[2L] <- 1
# par(mar = mar)
# plot.new()
plot.window(xlim = c(0, 1), ylim = range(levels), xaxs = "i",
yaxs = "i")
rect(0, levels[-length(levels)], 1, levels[-1L], col = col)
if (missing(key.axes)) {
if (axes)
axis(4)
}
else key.axes
box()
}
#
# if (!missing(key.title))
# key.title
# mar <- mar.orig
# mar[4L] <- 1
# par(mar = mar)
# plot.new()
# plot.window(xlim, ylim, "", xaxs = xaxs, yaxs = yaxs, asp = asp)
# if (!is.matrix(z) || nrow(z) <= 1L || ncol(z) <= 1L)
# stop("no proper 'z' matrix specified")
# if (!is.double(z))
# storage.mode(z) <- "double"
# .Internal(filledcontour(as.double(x), as.double(y), z, as.double(levels),
# col = col))
# if (missing(plot.axes)) {
# if (axes) {
# title(main = "", xlab = "", ylab = "")
# Axis(x, side = 1)
# Axis(y, side = 2)
# }
# }
# else plot.axes
# if (frame.plot)
# box()
# if (missing(plot.title))
# title(...)
# else plot.title
# invisible()
#}
CodeIQコード銀行にあなたのコードを預けてみませんか？
- CodeIQコード銀行ではあなたのコードを財産と考えます。
- お預かりいただいたコードは、CodeIQコード銀行がしっかり評価し、フィードバックいたします。
- 当コード銀行にお預けいただいたコードは、企業がみてスカウトをかける可能性があります。
- 転職したい方や将来転職することを考えている方で、今の自分のスキルレベルを知りたい方はぜひ挑戦してみてください。
- 企業からスカウトがきたら困る人は挑戦しないでください。
興味を持った方はこちらからチャレンジを！
エール大学経済学部博士課程。政治経済学・統計学が専門。日本の新聞記事のテキスト分析とマーケティング行動を研究している。将来の夢は囲碁電王戦に出場すること。