R:複数ポイントの一括ジオコーディング→ggplot2でマッピングして背景にGoogleMapを重ねる。

目的

1.ggmapライブラリのgeocode関数を使い、複数の住所のジオコーディングを一括で行います。

2.geom_pointでポイントをプロットし、ggmapでグーグルマップの背景に重ねます。店舗名のラベルをgeom_textでプロットします。

アウトプットイメージ

東京都港区　スターバックスの店舗分布

データ準備

今回は、スターバックスの店舗検索より東京都港区の店舗を検索した結果のhtmlファイルから、店舗名ファイルと店舗住所ファイルを2つ作成しました。きれいな住所データベースではなくwebから取得できるデータをそのまま使い、どのぐらいエラーなく使えるのかを検証してみました。

店舗名ファイル：shopname_minatoku.tsv

1	虎ノ門２丁目タワー店
2	虎ノ門駅前店
3	TSUTAYA神谷町駅前店
4	新橋駅前店
5	TSUTAYA新橋店
6	芝大門店
7	ゆりかもめ新橋駅店
...

店舗住所ファイル:id_store_address2.tsv

1	東京都港区　虎ノ門2-3-17　虎ノ門２丁目タワー
2	東京都港区　虎ノ門1-1-23　ウンピン虎ノ門ビル　1F
3	東京都港区　虎ノ門5-2-5　神谷町MTコート
4	東京都港区　新橋2-19-3　カシケイビル　1F
5	東京都港区　新橋4-11-1　TK新橋ビル
6	東京都港区　芝公園2-3-4　ホテルコンソレイユ芝・東京
7	東京都港区　東新橋1-1-8
...

データの読み込み

上記２つのデータを読み込みます。
普通にread.tableで何もオプションを付けずに読み込むと、住所や店舗名は文字列ではなくfactorとして読み込まれてしまうため、読み込んだ後にas.characterで型の変換をしています。str（データフレーム）で型を確認します。
factor型だと、ジオコーディングではErrになります。

# 1)ポイントデータの読み込み
library(ggmap)
setwd("~/Documents/GIS/Data/Point_data")

stb_address<-read.table("id_store_address2.tsv",header=F,sep="\t")
str(stb_address)
colnames(stb_address)<-c("id","address")
stb_address$address <- as.character(stb_address$address)
str(stb_address)

shopname<-read.table("shopname_minatoku.tsv",header=F,sep="\t")

colnames(shopname)<-c("id","shopname")
shopname$shopname<-as.character(shopname$shopname)
str(shopname)
head(stb_address)

一括ジオコーディング

ggmapライブラリのgeocode関数については、以前の記事(R:ggplot2でシェープファイルから人口マップを描画し背景にGoogleMapを重ねる。)でも書きました。住所を渡すと、経度、緯度の数値を戻すAPIで、ピンポイントの精度で結果が返ります。しかし住所の文字列により不正確な経緯度が返されることもあります。

# 結果を格納するデータフレーム
LLmatrix<-array("",dim=c(nrow(stb_address),3))
# ジオコーディングがエラーで返ってきた場合の、レコード番号と住所を保持するデータフレーム
err_id<-array("",dim=c(nrow(stb_address),2))
# k--エラーの個数をカウントするカウンタ
k<-0
for (i in 1:nrow(stb_address)){
  str_address<-stb_address[i,2]
  LonLatData<-geocode(str_address,source="google")
  if(is.na(LonLatData[1])){
    k<-k+1
    err_id[k,1]<- i
    err_id[k,2]<-str_address
  }
    LLmatrix[i,1]<-shopname[i,2] #店舗名
    LLmatrix[i,2]<-LonLatData[1,1] #経度
    LLmatrix[i,3]<-LonLatData[1,2] #緯度
}
k
err_id<-err_id[1:k,]
colnames(LLmatrix)<-c("shop_name","longitude","latitude")
colnames(err_id)<-c("id","address")
# LLmatrixがなぜかmatrixでdataframeではないので、データフレームにしておく
LLdf<-as.data.frame(LLmatrix)
err_id

• err_idより、結果がNullになってしまった（ジオコーディング失敗）住所については、手動で住所を直してジオコーディングします。

• 例えば、「芝公園2-3-4　ホテルコンソレイユ芝・東京」では失敗しますが、「芝公園2-3-4」では成功、ということがわかりました。建物名によりうまく返ってきたり失敗したりがあるようです(これは例です）。

• 住所データの加工ができた時点で、web掲載の住所は、全角数字があったり全角スペースがあったりすることがわかります。全角数値でもジオコーディングはできるようですが、何がエラーの原因になるか分からないので統一させたほうが良さそうです。

<参考>　ジオコーディング失敗の住所について手動で取得しなおし、結果データフレームの値を更新する。

# 経緯度が取得できなかったレコードのレコード番号と、住所表記を確認する
# err_id
# i<-as.numeric(err_id[1,1])
# #恐らく住所に建物名があるのが悪いと判断し、建物名を取る
# #自動化するには、err_id[1,2]の単語数を取得して、2個目まで（住所番地まで）でカットするようにする
# LonLatData<-geocode("東京都港区　東新橋1-5-2",source="google")
# 返された経緯度の住所を確認するrevgeocode
# revgeocode(c(lon=LonLatData[1,1],lat=LonLatData[1,2]))
# 結果のデータフレームの経緯度の値を更新する
# LLmatrix[i,2]<-LonLatData[1,1] #経度
# LLmatrix[i,3]<-LonLatData[1,2] #緯度

ggplot2で複数のポイントを同時にマッピングする

１）経度、緯度の型を、ggmapが返した同じ精度で数値に変換します。

ｓｔｒ(LLdf)としてみれば分かりますが、経度、緯度のフィールドがfactorになっています。
longitude_n,latitude_nというフィールドに数値変換したものを格納します。

sample.r

# 経緯度の表示を小数点以下最大10桁までに変える
options(digits=10) 
str(LLdf)
LLdf$longitude_n<-as.numeric(as.character(LLdf$longitude))
LLdf$latitude_n<-as.numeric(as.character(LLdf$latitude))
head(LLdf)

２） ggplot2でポイントをプロット

とても簡単にプロットができます。この段階ではグーグルの背景はありません。

# --ggplot2でポイントをプロットする
library(ggplot2)
g<-ggplot() +  geom_point(data=LLdf, aes(x=longitude_n, y=latitude_n), color="red")
plot(g)

経度・緯度のプロット

３）背景にgoogleマップを重ねます。

# 簡単に港区を中心とするggmapオブジェクトを取得
map <- get_map(location = "Minato-ku,Tokoy", zoom = 14,maptype="hybrid")
mp <- ggmap(map) +geom_point(aes(x = longitude_n, y = latitude_n), 
                                    data = LLdf, alpha = .9,
                                    color="red")
plot(mp)

-- zoomオプションで指定するマップのスケール（解像度）は、calc_zoom（）というggmapの関数で得ることもできます。

４）店名を重ねてプロットします。geom_text（）を使います。

# geom_textでhjustを0とすると、ポイントの真上にラベルを配置して見た目がよくない
# ラベルを少し右にずらす。hjustにマイナスの値を設定--ここははっきり言って試行錯誤

# check_overlap=TRUEすると、重なっているラベルを取る
mp<-mp+geom_text(data=LLdf,aes(x=longitude_n, y=latitude_n, label=shop_name), 
            family="Osaka",size = 5, vjust = 0, hjust = -0.05,color="Yellow",
            check_overlap = TRUE,
            inherit.aes=TRUE)
plot(mp)

完成マップ

備考

• ジオコーディング処理が一括で可能な件数について1回何個までという制限があるはずですが詳しくは御自身でご確認下さい。
• ジオコーディングの結果が正確かどうか、revgeocode関数で経緯度→住所に再度変換し、同じかどうかを確認できます。
• 準備したggmapオブジェクトの範囲からはみ出てしまうポイントがあると警告メッセージが表示されます。

参考にしたページ

• How to draw good looking maps in R(ちょっと古いですが役に立ちます）
• Plotting Russian AiRstRikes in SyRia
• ggplot2のドキュメント