Rで複数のシェープファイルを結合する

-
地図で可視化する際の行政界のシェープは、統計局のものをよく利用させてもらっています。↓以下参考。

e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ

統計局で公開されているシェープは比較的小さい単位に分かれています。例えば、福岡市だと、区単位（東区、博多区、中央区、南区、西区、城南区、早良区）の7つのシェープに分かれています。

こういうのをまとめて１つにしたい時ってありますよね。

QGISのようなGUIのあるGISソフトでやる方法もありますが、数が多くなるとしんどい。ということで、Rのスクリプトでシェープの結合をする方法を紹介します。

ざっくり言うと、readShapePolyで読み込んで、spRbindで結合させるだけなんですけどね。

でも、そう簡単にはいかない時があるんですよね。というのも、ポリゴンのIDがかぶっているとspRbindでエラーが出てしまう。例えば、東区のシェープの中のあるポリゴンに「1」というIDが付いていて、博多区のシェープの中のあるポリゴンにも「1」というIDが付いていると、この２つを結合させようとする際にエラーになってしまいます。

なので、東区のポリゴンはのIDは、"1.1", "1.2", "1.3", ...、博多区のポリゴンのIDは"2.1", "2.2", "2.3", ... みたいにリネームしたあとに、結合させましょうというのが、↓このスクリプトのポイントでございます。

library(maptools)
 
# シェープのファイル名を取得
shape_files <- list.files(path="shape_folder", pattern=".*\\.shp", full.names=T)
 
# 緯度経度、世界測地系
pj <- CRS("+proj=longlat +datum=WGS84")
 
# ポリゴンのIDが重複しないように前処理したあと、
# リストに格納しておく
shape_list <- list() # 空のリストを用意
 
for ( i in 1:length(shape_files) ) {
  # シェープの読み込み
  shp <-readShapePoly(shape_files[i], proj4string=pj)
 
  # id文字列のベクトルを作る
  ids <- paste(i, 1:length(shp@polygons), sep=".")
 
  # polygonsのIDを変更する
  for ( j in 1:length(shp@polygons) ) {
    shp@polygons[[j]]@ID <- ids[j]
  }
 
  # dataの行名もそれに対応して変更する
  row.names(shp@data) <- ids
 
  # リストに追加する
  shape_list <- c(shape_list, shp)
}
 
# リストに格納してある、シェープを結合する
merged <- shape_list[[1]]
for ( i in 2:length(shape_list) ) {
  merged <- spRbind(merged, shape_list[[i]])
}
 
# シェープの書き出し
writePolyShape(merged, "merge.shp")


↓あとは、おまけみたいなもんです。

ちゃんとくっついているかplotしてみます。ifelseでごちゃごちゃ書いているのは、区ごとに色分けするためです。

ややこしいのが海の部分で、ここを同じ色で塗りたいのですが、区によってMOJIの属性（通常は町丁目名が入っている）に「博多湾」と入っていたり、「博多港」と入っていたり、「水面」と入っていたり、NULLが入っていたりする。博多湾とか博多港は分かるけど、「水面」ってなんだよ。

あと、地図に文字をプロットするようなときには、文字が大きいと重なっちゃうし、小さいとつぶれてしまうことがよくあるので、PDFに出力するとベクタとして書き出せるので吉です。

library(RColorBrewer)
cols <- brewer.pal(8, "Pastel2")
 
merged$col <- ifelse ( merged$MOJI=="博多湾" |
                       merged$MOJI=="博多港" |
                       merged$MOJI=="水面"   |
                       is.na(merged$MOJI)       , cols[1],
              ifelse ( merged$CSS_NAME=="東区"  , cols[2],
              ifelse ( merged$CSS_NAME=="博多区", cols[3],
              ifelse ( merged$CSS_NAME=="中央区", cols[4],
              ifelse ( merged$CSS_NAME=="南区"  , cols[5],
              ifelse ( merged$CSS_NAME=="西区"  , cols[6],
              ifelse ( merged$CSS_NAME=="城南区", cols[7],
              ifelse ( merged$CSS_NAME=="早良区", cols[8], "NG" ))))))))
 
pdf("out.pdf", family="Japan1GothicBBB", width=11.69, height=16.54)
plot(merged, col=merged$col, lwd=0.1)
text(merged$X_CODE, merged$Y_CODE, merged$MOJI, cex=0.1)
dev.off()


↓こんな感じで表示されます。

福岡市の行政界地図

↓PDFファイル自体も置いておきますね。

福岡市の行政界地図（PDFファイル）

コメント

コメントを投稿

このブログの人気の投稿

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...
続きを読む