東北沖地震について

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

余震は減ってきている．ただ，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()

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

マグニチュード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()

うーん，わかりづらい．
ちなみに震度別．

# 震源地別震度ヒストグラム
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()

これも，わかりづらい．
やはり日毎に場所ごとに分けないといけないか．
そこで，緯度経度上にマグニチュードデータをプロットしてみる．

# 緯度経度でプロット
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()

時間毎に表示する動画にしないと，感覚的にはつかめない．