e-Stat(統計局)で公開されているShapeファイルを、Rで表示する

-
まずはShapeファイルを入手するところから。

総務省統計局がe-Statというサイトを運営していて、そこに全国各地域のShapeファイルが公開されています。

地図で見る統計（統計ＧＩＳ)


(1) 「平成２２年国勢調査（小地域） 2010/10/01」を選択
(2) 「男女別人口総数及び世帯総数」にチェック
(3) 「統計表各種データダウンロードへ」ボタンをクリック



(1) 「都道府県」を選択
(2) 「市区町村」を選択
(3) 「検索」ボタンをクリック
(4) 境界データの欄に表示されている
    「世界測地系緯度経度・Shape形式」のデータをクリック

ダウンロードしたzipファイルを展開すると、

　xxx.shp
　xxx.dbf
　xxx.shx
　xxx.prj

という、拡張子だけが違うファイルが4つできるんですが、これがshapeデータの1セットになります。

それぞれのファイルの中身は、シェープファイルとは？ | 用語集とGISの使い方 | 株式会社パスコによると、


　shp：    図形の座標が保存
　dbf：    属性の情報が保存
　shx：    shpの図形とdbfの属性の対応関係が保存

だそうです。あと、「.prj」は測地系の情報が入っています。（テキストとして開くとそれっぽい文言が入っている）

では単純にshapeファイルをRで表示させてみましょう。

境界情報だけをプロット
# GIS関係のパッケージ
library(maptools)
 
# 座標系の情報を持つオブジェクトを作る
# longlatは緯度経度指定、WGS84は世界測地系であることを表す
pj <- CRS("+proj=longlat +datum=WGS84")
 
# シェープファイルの読み込み
x  <- readShapePoly("h22ka14137.shp", proj4string=pj)
 
# 表示
plot(x)


あら簡単。実質4行で表示できました。

じゃあ、いろいろやってみましょう。

coordinates()という関数は、ポリゴンの重心の位置を返してくれます。

> coordinates(x)
        [,1]     [,2]
0   139.4639 35.62169
1   139.5178 35.62177
2   139.5138 35.62332
3   139.5057 35.62163
4   139.4617 35.61308
         ...      ...

そのままプロットすると↓こんな感じですね。

重心の位置をプロット
plot(x)
points(coordinates(x))

このshapeを読み込んだオブジェクト（SpatialPolygonsDataFrame）には、境界のポリゴンだけでなく、属性情報も入っています。

> head(x@data)
        AREA PERIMETER ... KEN_NAME ... MOJI           ...
0  206286.40  2421.068 ... 神奈川県 ... はるひ野５丁目 ...
1   67760.92  1765.906 ... 神奈川県 ... 細山６丁目     ...
2  119563.80  1457.849 ... 神奈川県 ... 細山７丁目     ...
3  225729.80  2633.794 ... 神奈川県 ... 細山           ...
4 2017504.00 11135.110 ... 神奈川県 ... 黒川           ...
5   95170.74  1569.407 ... 神奈川県 ... 細山８丁目     ...
         ...       ... ... ...      ... ...            ...


どうやら、この「MOJI」というカラムにそのポリゴンに対応する地名が入っているようです。これを地図上にプロットしてみましょう。プロット位置はさっきのやり方で求めた重心の位置にしてみます。あと、MOJIは因子型で入っているので、as.character()で文字列に変換してやる必要がありました。

地名を表示させたが、文字サイズが大きすぎた
plot(x)
text(coordinates(x), as.character(x@data$MOJI))

うわ、ぐちゃぐちゃ。

字を小さくしてみよう↓

文字サイズが小さいと見にくい・・・
plot(x)
text(coordinates(x), as.character(x@data$MOJI), cex=0.5)

今度は、字が小さくて見づらいですね。

画面に出すだけのときは、描画用のウィンドウをググッと広げると、ポリゴンは大きく、文字はそのままで表示される、見やすくなります。

直接、画像ファイルとして出力し、その際に解像度を高めに指定するというやり方もあります。


高解像度のpngに出力（拡大してみればいい感じです）
library(png) # png出力のためのパッケージ
 
# 今から出力先はpngファイルだぞ、という指示
png("output.png", width=2000, height=2000)
 
plot(x)
text(coordinates(x), as.character(x@data$MOJI))
 
# pngファイルを忘れずに閉じる
dev.off()


解像度やら文字の大きさやらは、お好みで調整してみてください。

他にもいろいろと遊べるのですが、長くなったので、今回はここまで。■

[2016年3月30日追記]
広いエリアの行政界シェープファイルのようなものを画像ファイルに出力するのは、けっこう無理があります。そのようなときはPDFファイルに出力すると、いい感じになります↓
Rでshapeを出力するときはPDFファイルにすると便利 - Rプログラミングの小ネタ
■

コメント

コメントを投稿

このブログの人気の投稿

Rのグラフで軸の目盛りの刻み幅を変更する方法

-
イメージ
Rでplotなどを使ってグラフを描くとき、x軸やy軸の目盛りは勝手に調整してくれて、大抵の場合はそれで問題ないのですが、たまにちょっと変えたい時があります。そのたびに必死で検索して調べているような気がするので、ここに書き留めておきます。 （最初の方の例はcurve関数を使ってますが、plot関数でも同様の方法でいけます） もちろん、特に指定しなくても、目盛りと目盛りのラベル（例えば、x軸の -4, -2, 0, 2, 4）は書いてくれます↓ 指定しなくても目盛りは表示される # 目盛りの幅に関しては指定しない curve(dnorm, xlim=c(-4,4)) 上記では、2刻みになっていますが、例えば、1刻みにしたいなあという場合↓ 目盛りが1刻みになるように指定 # -4から4までを8区間に分ける（1刻み） curve(dnorm, xlim=c(-4, 4), xaxp=c(-4, 4, 8)) さらに刻みを細かくして0.5にしてみましたが、目盛りのラベルは適宜間引かれてしまうようです↓（作図領域を横に伸ばせばラベルが表示されるようになります） 目盛りが0.5刻みになるように指定 # 16区間に分けると0.5刻みになるが # すべての目盛りに数値ラベルが書かれるわけではない curve(dnorm, xlim=c(-4, 4), xaxp=c(-4, 4, 16)) いくつに分割するかちゃんと考えないと、↓こんなことになることも・・・ 割り切れないところに目盛りがきてしまう例 # 8の幅を9で割るとか・・・ curve(dnorm, xlim=c(-4, 4), xaxp=c(-4, 4, 9)) 次は違うやり方。 最初は軸の目盛りを書かないようにしておいて、あとから、axis関数をつかって目盛りを書き足すという方法です↓ axis関数で目盛りを後から追記したもの（見た目おんなじだけど・・・） curve(dnorm, xlim=c(-4, 4), xaxt="n") # 最初は目盛りを書かない axis(side=1, at=-4:4) # atに目盛りをベクトルで指定する side=1というのがx軸の意味で、y軸に追記したい場合はsi...
続きを読む

Rで繰り返しを含む数列の生成（rep関数、seq関数）

-
Rで、「123123123...」とか「111222333...」とか「13579...」（奇数）とか「02468...」（偶数）とか、そんな感じの繰り返しを含むような数字の列を作りたい場合、1行でさらっと書けてしまいます。 【繰り返し】 > rep(7, 10)  [1] 7 7 7 7 7 7 7 7 7 7 「○を○個だけrepeat」って感じで。 【1ずつ増加】 > 1:7  [1]  1  2  3  4  5  6  7 コロンでつなげば「○から○までインクリメント」 【○ずつ増加】 > seq(1, 7, 2) [1]  1  3  5  7 sequence関数で、「○から○まで○ずつ増加」 【1ずつ増加する数列を繰り返す】   > rep(1:7, 3)  [1] 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 repeat関数の第1引数に、「1ずつ増加する数列」を指定すればOK。 【○ずつ増加する数列を繰り返す】 > rep(seq(1,7,2), 3)  [1] 1 3 5 7 1 3 5 7 1 3 5 7 repeat関数の第1引数に、「○ずつ増加する数列」を指定すればOK。 【それぞれ○個】   > rep(1:7, each=3)  [1] 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 【○ずつ増加、それぞれ○個】 > rep(seq(1,7,2), each=3)  [1] 1 1 1 3 3 3 5 5 5 7 7 7 コロンで指定していたところを、repにするだけなので簡単ですね。 【○ずつ増加、○回繰り返す、それぞれ○個】 > rep(seq(1,7,2), 2, each=3)  [1] 1 1 1 3 3 3 5 5 5 7 7 7 1 1 1 3 3 3 5 5 5 7 7 7 rep関数の第2引数は繰り...
続きを読む

Rで複数のヒストグラムを比較する方法あれこれ

-
イメージ
2つ以上のヒストグラムを1つのグラフ上に描画して、比較したいようなときってありますよね。 例えば、集団aのテストの点数と、集団bのテストの点数が、それぞれ1000件ずつあるとしましょう。 # サンプルデータを作っておく a <- rnorm ( 1000 , 40 , 10 ) b <- rnorm ( 1000 , 60 , 10 ) borderオプションを使って、ヒストグラムの線に色をつけるってのは1つのやり方ですね↓ # 線に色をつける hist ( a , breaks= seq ( 0 , 100 , 5 ) , border= "red" , xlim= c ( 0 , 100 ) , main= "" , xlab= "" ) hist ( b , breaks= seq ( 0 , 100 , 5 ) , border= "blue" , add=T ) borderオプションでヒストグラムの線に色をつける なんとなく見えてはいるのですが、重なった縦線は後から描画した方で上書きされてしまい、そのへんがちょっと見にくいですね。 colオプションで面に色をつけることもできますが・・・ # 色で塗りつぶす hist ( a , breaks= seq ( 0 , 100 , 5 ) , col = "red" , xlim= c ( 0 , 100 ) , main= "" , xlab= "" ) hist ( b , breaks= seq ( 0 , 100 , 5 ) , col = "blue" , add=T ) colオプションでヒストグラムを塗りつぶす ↑普通の色で塗っちゃうと重なったところが見えなくなってしまいます。 そういうときは、透過色を使えばいいです↓ # 透過色を使う hist ( a , breaks= seq ( 0 , 100 , 5 ) , col = "#FF00007F" , x...
続きを読む