2014年11月22日土曜日

Rで度数分布表を作る

まずは今回のサンプル用として、身長っぽいデータを作ってみます。

a <- round(rnorm(30, mean=170, sd=5), 1)
a
 
 [1] 173.0 168.6 168.5 164.2 167.3 170.6 162.7 168.0 170.8
[10] 175.7 166.0 166.5 162.4 172.6 170.1 170.2 164.2 167.1
[19] 163.8 163.4 168.1 168.7 171.1 164.6 166.3 177.0 170.0
[28] 173.3 170.7 169.0


このaに対して、度数分布表を作りたい。

質的データ(カテゴリカルデータ)なら、table関数を使って度数を集計できるのですが、量的データに対してtableを使うと、↓こんな感じになってしまいます。

table(a)
 
a
162.4 162.7 163.4 163.8 164.2 164.6   166 166.3 166.5 167.1 
    1     1     1     1     2     1     1     1     1     1 
167.3   168 168.1 168.5 168.6 168.7   169   170 170.1 170.2 
    1     1     1     1     1     1     1     1     1     1 
170.6 170.7 170.8 171.1 172.6   173 173.3 175.7   177 
    1     1     1     1     1     1     1     1     1


量的データだとうまくいかないですね。

Rには、お馴染みhist関数があるんで、ヒストグラムは簡単に描けます。

hist(a)




実は、このhist関数の戻り値のオブジェクトに各階級の度数が入っているんです。知ってました?

h <- hist(a)
h
 
$breaks
[1] 162 164 166 168 170 172 174 176 178
 
$counts
[1] 4 4 5 6 6 3 1 1
 
$density
[1] 0.06666667 0.06666667 0.08333333 0.10000000 0.10000000
[6] 0.05000000 0.01666667 0.01666667
 
$mids
[1] 163 165 167 169 171 173 175 177
 
$xname
[1] "a"
 
$equidist
[1] TRUE
 
attr(,"class")
[1] "histogram"


breaksが階級を区切る値で、countsが度数ですね。

つまり、155~160の範囲の度数が1、160~165の範囲の度数が7、・・・のような意味です。

細かく言うと、境目のところは

  155 < 階級に含まれる値 ≦ 160

のように扱われています。

これらの値を使えば、度数分布表が作れそうです。

前置きが長くなりましたが、↓こんな感じでどうでしょう?

h <- hist(a)
n <- length(h$counts) # 階級の数
class_names <- NULL # 階級の名前格納用
for(i in 1:n) {
  class_names[i] <- paste(h$breaks[i], "~", h$breaks[i+1])
}
frequency_table <- data.frame(class=class_names, frequency=h$counts)
frequency_table
 
       class frequency
1 162164         4
2 164166         4
3 166168         5
4 168170         6
5 170172         6
6 172174         3
7 174176         1
8 176178         1


データフレームの形で度数分布表ができました。

階級の取り方を指定したいときには、hist関数を呼ぶときにbreaksオプションで指定すればOKです。↓こんな感じ。

h <- hist(a, breaks=seq(155, 185, 10))


データフレームの形になってしまえば、あとは好きなように使えますよね。




2014年11月3日月曜日

Rを使ってシェープファイルをKMLファイルに変換する(土砂災害危険箇所マップ)

↓こちらで公開したKMLファイルは、シェープファイルをRで変換処理したものなのですが、その手順です。

神奈川県の土砂災害危険箇所ハザードマップ(Google Earth閲覧用KML形式): 主張

元にしたシェープは↓こちらで入手しました。

国土数値情報ダウンロードサービス

<災害・防災>のカテゴリに、「土砂災害危険箇所」というのがあります。シェープファイルは都道府県別に分かれて公開されています。

データ形式のタブ(?)で、「JPGIS2.1」のところを選択すると、JPGIS2.1形式とともになぜかShapeファイルも置いてあるという、微妙なトラップがあるので要注意です。

このShapeファイルをRでKMLに変換するんですが、スクリプトはごく簡単です↓

# 必要なパッケージの読み込み
library(maptools)
library(plotKML)
 
# 座標系の情報を持つオブジェクトを作る
# longlatは緯度経度指定、WGS84は世界測地系であることを表す
pj <- CRS("+proj=longlat +datum=WGS84")
 
# シェープファイルの読み込み
hazard.area <-readShapePoly("xxx.shp", proj4string=pj)
 
# もし、表示して確認したいときはplotすればOK
# plot(hazard.area)
 
# KMLファイルとして出力する
kml(hazard.area, file = "xxx.kml")


あとは、GoogleEarthで開けば閲覧できます。

↓こんな感じですね。
国交省が公開している「土砂災害危険箇所」シェープをRでKMLに変換後、Google Earthで表示したもの



閲覧時の色の変更などは、冒頭のリンク先を参照してみてください。■

[2014年11月23日追記]東京都、埼玉県、千葉県のKMLファイルも追加しました。あと、Google Earthの設定で色を変更しなくて済むよう、KMLファイル内のタグで色を赤に変えときました↓
土砂災害危険箇所のKMLファイル(東京、神奈川、埼玉、千葉): 主張