Help us understand the problem. What is going on with this article?

R:ggplot2でシェープファイルから人口マップを描画し背景にGoogleMapを重ねる。

More than 3 years have passed since last update.

目的

オープンな地図データを使って、Rでコロプレスマップ(塗り分け地図)を作成します。シェープファイルをggplot2で読めるデータフレームに変換しggplot2で図化、背景にグーグルマップを重ねます。

アウトプットイメージ

image

ggmapライブラリでできることを確認する。

1)まずはggmapライブラリを使いジオコーディングを試します。

港区 町丁目名(虎ノ門4丁目) 番地名(1-8)まで入力すると、経緯度の結果が返されました。

library(ggmap)

LonLatData<-geocode("港区 虎ノ門 4丁目1−8",source="google")

image

2)経緯度を渡してどのレベルの住所が返されるかを確認します。

取得した経緯度をrevgeocode関数に渡すと、経緯度→住所変換ができます。geocodeコマンドの結果は、番地・建物レベルで合っていることが分かりました。

#緯度経度データから住所をグーグルから取得:revgeocodeコマンド
revgeocode(c(lon=LonLatData[1,1],lat=LonLatData[1,2]))

image

3)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)

4)ggmap関数でグーグルマップのオブジェクトを描画します。

#プロット
ggmap(GmapData)

image

5)座標をマーカーに登録し、マーカー付きのグーグルマップを描画します。

sample.R
#マーカーを準備する
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)

image

6)今度はマーカーではなく、中心をポイントにして描画します。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')

image

シェープファイルをRで描画する

次に、Rでシェープファイルをプロットします。
データは、Estatから取得した東京都港区の平成22年(2010年)版小地域境界データ(世界測地系版)をダウンロードして使っています。世界測地系緯度経度・Shape形式をダウンロードします。
※世界測地系じゃないとグーグルマップとオーバーレイできません。(直角座標系版のシェープファイルををggmapでグーグルマップに重ね、グーグルマップの投影法を直角座標系に変えることは可能なのかもしれませんが現時点では調べていませんのでここでは世界測地系版を使います)。

1) maptoolsライブラリでシェープファイルを読み込みます。

#Rでシェープファイルをプロットする
#シェープファイルが存在するフォルダにアクセスする
setwd("/Users/foo/Documents/GIS/shape/世界測地系_WC13103")
library(maptools) #readShapePolyを使う
#maptoolsを使って、shapefileを読み込む。
shpfile <- "h22ka13103.shp"
sh <- readShapePoly(shpfile)
plot(sh)

港区の町丁目境界がプロットされました。
image

2)シェープファイルの属性テーブルの中身を確認します。

Rのエディタ上でオブジェクト名sh@と入力すると、オブジェクトの下位クラスに何があるかが分かります。ここから、data が属性テーブルであろうと想像できます。

image

sh@data$と入力して、属性テーブル内に存在する変数リストを確認できます。
image

KEY_CODEは、東京都のコード13と港区のコード103+KEYCODE1と同じ6桁の小地域コードを合わせた11桁コードです。ここでは地区のコードとして使います。

3) カラーパレットを作成し、シェープファイルの属性テーブルに付いている人口の値をカラーレンジでプロットします。

属性テーブルのJINKOフィールドに地区毎の人口データがあります。このデータを使ってコロプレスマップを描画します。
colorRampPalette関数で白から赤にかけての128段階のカラーパレットを用意し、人口データを(各値ー最小値)/レンジで0〜1に基準化し、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)

こんなマップができました。
image

ggplot2を使い、グーグルマップの上にポリゴンを重ねます。

1)シェープファイルのポリゴンをggmapが理解できるデータフレームに変換します。ggplotでグーグルマップ上にシェープファイルをオーバーレイします。

##fortifyコマンドggmapを使うために、shapefileのポリゴンをggmapが理解できる形式に変換する。

points<- fortify(sh,region='KEY_CODE')
#pointsの中身を確認する,
head(points)

pointsの中身はこんな感じでした。シェープファイルの形式ではあった属性テーブル内の変数は無くなり、経緯度と地域ID、group等図形情報として必要な情報だけになっています。

image

2) 人口データと地区コードのみ、シェープファイルから別のデータに切り出します。

元々のダウンロードしたシェープファイルの属性に含まれていた人口データ(JINKOフィールド)と、地区コードとしてKEY_CODEフィールドだけ、シェープファイルオブジェクトから切り出します。

#shapeファイルの属性テーブルに含まれる人口データを描画するために、KEY_CODEと人口のみcsvにエクスポートする
head(sh@data[,c("KEY_CODE","JINKO")])
attrib_df<-sh@data[,c("KEY_CODE","JINKO")]

3) ポリゴンのデータフレーム(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フィールドが結合されています。

image

4)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)

image

5)グーグルマップ上に塗り分け図(コロプレスマップ)を重ねます。

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')

image

参考にしたページ

迷った所

  • シェープファイルをポリゴン化した時に元のシェープファイルに存在していた属性変数は消えてしまいます。そのため、別のデータフレームに地区コードと人口変数のみ切り出して、ポリゴンのデータフレームに結合する必要がありました。

未解決の問題

  • シェープファイルの地区名が文字化けしています。地図のラベルに使えないのは残念ですが、多分解決法はあるでしょう。
  • シェープファイルに投影法が設定されていた場合のグーグルマップのオーバーレイ方法
Why do not you register as a user and use Qiita more conveniently?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away