東北沖地震について

f:id:joker1110:20110502025439j:image:w360:left
3/11,東北沖地震が発生して以来,自分の無力さを痛感した.何かをしなければ,と思っても,せいぜい募金ぐらい.

もちろんお金は復興のためにも重要だけれど,こういった被害が二度とないようにしたいという思いもある.そこで,今回の地震について,統計的な視点から何かできないかと考えた.

まずは,未だに続く地震そのものについて.

地震発生のデータは日本気象協会が地震情報として公開しているものを利用する.

手順は,

  • データ取得には Java を利用.
  • データ図示には R を利用.

で行う.

ソースは長文なので別の機会に.

  • 地震データ取得用ソース GWT.java
  • csv ファイルの読み込み用ソース ReadCSV.java
  • csv ファイルの書き込み用ソース WriteCSV.java
  • (# ついでに,Gui パネル Panel.java Panel_abstract.java)

データは magnitude.csv というファイル名で保存.
これを R によって図示する.
図示には ggplot2 を利用.

#パッケージggplot2を読み込み
library(ggplot2)
library(splines)

#日本語対応フォント指定
quartzFonts(HiraMaru=quartzFont(rep("HiraMaruProN-W4", 4)))
par(family="HiraKakuProN-W3") # Mac専用

x = read.csv("magnitude.csv",header=TRUE, fileEncoding="SJIS")

# マグニチュードが欠損している行を取り除く
x <- x[!is.na(x$magnitude),]
x <- x[!is.na(x$power),]
x$depth <- -x$depth

# 日時の文字列をPOSIXctオブジェクトに変換する
x$datetime <- as.POSIXct(x$datetime)

# 日付のカラムを追加する
x$date <- as.Date(x$date, tz = "Asia/Tokyo")

# マグニチュードで時間軸にプロット
pdf("magnitude-to-date.pdf",width=13,height=8,family="Japan1GothicBBB")
h1 <- ggplot(x, aes(x=datetime, y=magnitude))
h1 <- h1 + geom_point(aes(colour=depth,size=power))
h1 <- h1 + geom_line(alpha=0.3)
h1 <- h1 + scale_x_datetime(format = "%m/%d")
h1 <- h1 + theme_grey(base_size=24)
h1 <- h1 + geom_smooth()
print(h1)
dev.off()

f:id:joker1110:20110502024043j:image:w360
余震は減ってきている.ただ,4/11あたりが気になる.

日毎に分けてヒストグラムにしてみる.

# 各日のマグニチュードヒストグラム
pdf("magnitude-by-date.pdf",width=13,height=8,family="Japan1GothicBBB")
h <- ggplot(x, aes(magnitude))
h <- h + geom_histogram(aes(fill=..count..))
h <- h + facet_wrap(~date)
h <- h + xlab("Magnitude")
h <- h + theme_grey(base_size=20,base_family="HiraMaru")
print(h)
dev.off()

f:id:joker1110:20110502024217j:image:w360
4/11-12が増加しているが,それ以降は落ち着いている.
全体でのヒストグラムを見てみる.

# 全日のマグニチュードヒストグラム
png("magnitude-all-date.png",width=1280,height=800)
h_all <- ggplot(x, aes(magnitude))
h_all <- h_all + geom_histogram(aes(fill=..count..))
h_all <- h_all + xlab("Magnitude")
h_all <- h_all + theme_grey(base_size=20,base_family="HiraMaru")
print(h_all)
dev.off()

f:id:joker1110:20110502024249j:image:w360
マグニチュード3と4.5あたりに二つの山ができている.
これだと日本全国をまとめて見ているので,場所ごとに分けてみる.

# 震源地別マグニチュードヒストグラム
png("magnitude-by-place.png",width=1280,height=800)
hi <- ggplot(x, aes(magnitude))
hi <- hi + geom_histogram(aes(fill=..count..))
hi <- hi + facet_wrap(~place)
hi <- hi + xlab("Magnitude")
hi <- hi + theme_grey(base_size=20,base_family="HiraMaru")
print(hi)
dev.off()

f:id:joker1110:20110502024314j:image:w360
うーん,わかりづらい.
ちなみに震度別.

# 震源地別震度ヒストグラム
png("SeismicIntensity-by-place.png",width=1280,height=800)
si <- ggplot(x, aes(power))
si <- si + geom_histogram(aes(fill=..count..))
si <- si + facet_wrap(~place)
si <- si + xlab("Seismic Intensity")
si <- si + theme_grey(base_size=20,base_family="HiraMaru")
print(si)
dev.off()

f:id:joker1110:20110502024349j:image:w360
これも,わかりづらい.
やはり日毎に場所ごとに分けないといけないか.
そこで,緯度経度上にマグニチュードデータをプロットしてみる.

# 緯度経度でプロット
png("jpnMap-to-magnitude.png",width=1280,height=800)
p <- ggplot(jpn2,aes(longitude, latitude))
p <- p + geom_path(alpha=0.1)
map <-  p + geom_point(aes(x=x$longitude, y=x$latitude,size=x$magnitude, colour=x$depth, alpha=0.1))
map <- map + scale_shape(solid = FALSE)
map <- map + scale_size(to = c(1, 20))
map <- map + xlab("longitude") + ylab("latitude")
map <- map + theme_grey(20)
print(map)
dev.off()

f:id:joker1110:20110502024412j:image:w360
時間毎に表示する動画にしないと,感覚的にはつかめない.