##目的
オープンな地図データを使って、Rでコロプレスマップ（塗り分け地図）を作成します。シェープファイルをggplot2で読めるデータフレームに変換しggplot2で図化、背景にグーグルマップを重ねます。
##アウトプットイメージ
##ggmapライブラリでできることを確認する。
１）まずはggmapライブラリを使いジオコーディングを試します。
港区 町丁目名（虎ノ門4丁目） 番地名(1-8)まで入力すると、経緯度の結果が返されました。
library(ggmap)
LonLatData<-geocode("港区 虎ノ門 4丁目１−８",source="google")
２）経緯度を渡してどのレベルの住所が返されるかを確認します。
取得した経緯度をrevgeocode関数に渡すと、経緯度→住所変換ができます。geocodeコマンドの結果は、番地・建物レベルで合っていることが分かりました。
#緯度経度データから住所をグーグルから取得：revgeocodeコマンド
revgeocode(c(lon=LonLatData[1,1],lat=LonLatData[1,2]))
３）get_googlemap関数を使いマップを取得します。
get_googlemapは、パラメータの経緯度を中心にしたグーグルマップをオブジェクトとして返します。
#get_googlemap関数 Google Static Maps API ver.2にアクセスする
GmapData<-get_googlemap(center=c(lon=LonLatData[1,1],lat=LonLatData[1,2]),
zoom=16,size=c(640,640),scale=2,format="png8",
maptype="roadmap",language="jpn",sensor=FALSE,
messaging=FALSE,urlonly=FALSE,filename="ggmapTemp",color="color",force=FALSE)
４）ggmap関数でグーグルマップのオブジェクトを描画します。
#プロット
ggmap(GmapData)
５）座標をマーカーに登録し、マーカー付きのグーグルマップを描画します。
#マーカーを準備する
pt = data.frame(lon<-LonLatData[1,1],lat<-LonLatData[1,2])
#マーカー付きの地図をプロットする
GmapData<-get_googlemap(center=c(lon=LonLatData[1,1],lat=LonLatData[1,2])
zoom=16,size=c(640,640),scale=2,format="png8",
maptype="roadmap",language="jpn",sensor=FALSE,
messaging=FALSE,urlonly=FALSE,filename="ggmapTemp",
color="color",force=FALSE,markers = pt)
#プロット
ggmap(GmapData)
６）今度はマーカーではなく、中心をポイントにして描画します。ggplot2ライブラリのgeom_point関数を使います。
#マーカーなしのグーグルマップ（markers=pオプションを書かない。）
GmapData<-get_googlemap(center=c(lon=LonLatData[1,1],lat=LonLatData[1,2]),
zoom=16,size=c(640,640),scale=2,format="png8",
maptype="roadmap",language="jpn",sensor=FALSE,
messaging=FALSE,urlonly=FALSE,
filename="ggmapTemp",color="color",force=FALSE)
#geom_point関数を使い、ポイントとgooglemapを同時にプロット。
ggmap(GmapData)+geom_point(aes(x=LonLatData[1,1],y=LonLatData[1,2]),
data=pt,size=6,colour='red')
##シェープファイルをRで描画する
次に、Rでシェープファイルをプロットします。
データは、Estatから取得した東京都港区の平成22年（2010年）版小地域境界データ（世界測地系版）をダウンロードして使っています。世界測地系緯度経度・Shape形式をダウンロードします。
※世界測地系じゃないとグーグルマップとオーバーレイできません。（直角座標系版のシェープファイルををggmapでグーグルマップに重ね、グーグルマップの投影法を直角座標系に変えることは可能なのかもしれませんが現時点では調べていませんのでここでは世界測地系版を使います）。
###１） maptoolsライブラリでシェープファイルを読み込みます。
#Rでシェープファイルをプロットする
#シェープファイルが存在するフォルダにアクセスする
setwd("/Users/foo/Documents/GIS/shape/世界測地系_WC13103")
library(maptools) #readShapePolyを使う
#maptoolsを使って、shapefileを読み込む。
shpfile <- "h22ka13103.shp"
sh <- readShapePoly(shpfile)
plot(sh)
港区の町丁目境界がプロットされました。
###２）シェープファイルの属性テーブルの中身を確認します。
Rのエディタ上でオブジェクト名ｓｈ@と入力すると、オブジェクトの下位クラスに何があるかが分かります。ここから、data が属性テーブルであろうと想像できます。
sh@data$と入力して、属性テーブル内に存在する変数リストを確認できます。
KEY_CODEは、東京都のコード13と港区のコード103＋KEYCODE1と同じ6桁の小地域コードを合わせた１１桁コードです。ここでは地区のコードとして使います。
###３) カラーパレットを作成し、シェープファイルの属性テーブルに付いている人口の値をカラーレンジでプロットします。
属性テーブルのJINKOフィールドに地区毎の人口データがあります。このデータを使ってコロプレスマップを描画します。
colorRampPalette関数で白から赤にかけての128段階のカラーパレットを用意し、人口データを（各値ー最小値）/レンジで０〜１に基準化し、127段階に合わせます。
# Scale the total population to the palette
head(sh@data)
#シェープファイルの属性テーブルのうち、JINKOフィールドを、pop変数に格納する
#JINKOフィールドをpopに格納する。
pop <- sh@data$JINKO
# カラーパレットを白から赤のレンジで128段階に分けたパレットを用意する
p <- colorRampPalette(c("white", "red"))(128)
palette(p)
#色を各地区の人口と最小の人口との差をとり、レンジで割ったものに・。。
cols <- (pop - min(pop))/diff(range(pop))*127+1
#plot関数は、プロットする。
plot(sh, col=cols)
こんなマップができました。
##ggplot2を使い、グーグルマップの上にポリゴンを重ねます。
###１）シェープファイルのポリゴンをggmapが理解できるデータフレームに変換します。ggplotでグーグルマップ上にシェープファイルをオーバーレイします。
##fortifyコマンドggmapを使うために、shapefileのポリゴンをggmapが理解できる形式に変換する。
points<- fortify(sh,region='KEY_CODE')
#pointsの中身を確認する,
head(points)
pointsの中身はこんな感じでした。シェープファイルの形式ではあった属性テーブル内の変数は無くなり、経緯度と地域ID、group等図形情報として必要な情報だけになっています。
###２) 人口データと地区コードのみ、シェープファイルから別のデータに切り出します。
元々のダウンロードしたシェープファイルの属性に含まれていた人口データ(JINKOフィールド）と、地区コードとしてKEY_CODEフィールドだけ、シェープファイルオブジェクトから切り出します。
#shapeファイルの属性テーブルに含まれる人口データを描画するために、KEY_CODEと人口のみcsvにエクスポートする
head(sh@data[,c("KEY_CODE","JINKO")])
attrib_df<-sh@data[,c("KEY_CODE","JINKO")]
###３) ポリゴンのデータフレーム(points)と、人口フィールドのデータフレーム(attrib_df)を結合します。
ポリゴンのデータフレームの地区コードはidというフィールド、属性データフレームではKEY_CODEが、結合するキーになります。
#いったん作成したpointsデータフレームに、attrib_dfをマージする。
#pointsのデータフレームは地域コードは[id]というフィールド名になっている
points2 <-merge(points,attrib_df,by.x="id",by.y="KEY_CODE",all.x=TRUE)
points2の中身を確認すると、ちゃんとJINKOフィールドが結合されています。
###４）ggmapライブラリのqmap関数で簡単にグーグルマップを取得し、地区境界に重ねます。
qmqp("住所")で簡単に港区を中心としたグーグルマップを作成します。
港区のグーグルマップだけプロットし、その後でggplot2ライブラリのgeom_polygon関数でpointsポリゴンを描画します。地区のボーダーラインは黒に指定します。
minatoku<-qmap("Minato-ku,Tokyo",zoom=13)
plot(minatoku)
#オーバーレイする data=pointsというのは、shapeファイルから取得したもの
minatoku + geom_polygon(aes(x=long,y=lat,
group=group, alpha=0.25),data=points, fill='white')
+geom_polygon(aes(x=long,y=lat, group=group), data=points, color='black', fill=NA)
###５)グーグルマップ上に塗り分け図（コロプレスマップ）を重ねます。
geom_polygonで人口データがマージされているpoints2データフレームを指定します。aesのオプションで、fill＝JINKOとして塗り分けする値のフィールド名を指定すると完成です。
#googlemapと一緒にpolygonを描画する場合 geom_polygom でPlot
minatoku + geom_polygon(aes(x=long,y=lat, group=group,fill=JINKO),
data=points2, color='black') +scale_fill_gradient(low='white', high='red')
##参考にしたページ
- 本ページはコチラのページPlotting Choropleths from Shapefiles in R with ggmap – Toronto Neighbourhoods by Populationのやり方をほぼ忠実になぞってみました。ただし例題とした町丁目別の国勢調査地図は、Estatからダウンロードしたものを使いました。
##迷った所
- シェープファイルをポリゴン化した時に元のシェープファイルに存在していた属性変数は消えてしまいます。そのため、別のデータフレームに地区コードと人口変数のみ切り出して、ポリゴンのデータフレームに結合する必要がありました。
##未解決の問題
- シェープファイルの地区名が文字化けしています。地図のラベルに使えないのは残念ですが、多分解決法はあるでしょう。
- シェープファイルに投影法が設定されていた場合のグーグルマップのオーバーレイ方法