地震データの経時分析
東日本大震災（東北地方太平洋沖地震）が起こったのがきっかけで、地震データのモニタリングができないかなと思いました。
サーバー負荷も考え、webから直接データを取得するのではなくページの表をエクセルにコピペしてデータに起こしています。
データソースは気象庁のtenki.jp。
まず取得したのは、2/27～3/22 19:20までの676件。
tenki.jpでは東北地方太平洋沖地震のマグニチュードらしきデータがM7になっていたので、M9に修正しています。
マグニチュードが不明な観測値は削除しています。
震源地をいくつかの地域に分類して、マグニチュード毎に分割表にすると次のようになります。
当然ですが関東と東北が桁違いに多いです。
このデータを日付とマグニチュードで図にしてみます。
これは地域毎に色分けしており、右に行くほど最近の地震を示しています。
M9の地震のところに点線を引いています。
ピンク色が東北地方で黄緑色が関東地方ですが、東北地方太平洋沖地震が起こった後は似たように地震が起こっています。
これは同じプレート上で余震が起こっている事を示しているのでしょう。
またもう1つ気になるのは、3/11の2日前くらいから急に地震が増えていること。
しかもM6クラスが連発しています（M5とM6のところに点線で補助線を引いています）。
別のデータで過去3年間分チェックしましたが、このようにM6クラスが2日間連続している状況はありませんでした。
これを巨大地震の前兆と考えるかどうかは、もっと多くのデータを統計的に見ないとダメでしょうが、1つの指標になる可能性はあります。
このプロットに平滑化曲線（スプライン）をあてはめるとこのようになります。
スプラインは極端な値に影響を受けやすいので解釈にはとても注意が必要です。
しかしこれも1つの目安にはなるかと思います。
このように予測は難しいにせよ、常時モニタリングができるだけでも、被害の大きさが少しは違うのかもしれませんね。。。
↓以下Rのプログラムです。
Jishin <- read.csv("/地震.csv", as.is = T) Jishin$Mag <- as.numeric(substr(Jishin$M, 2, 5)) Jishin$Time <- as.numeric(strptime(paste(paste(Jishin$月, gsub("日", "", Jishin$日), substr(Jishin$年, 3, 4), sep = "/"), Jishin$発生時刻), "%m/%d/%y %H:%M:%S")) Jishin$Eng <- 10 ^ (4.5 + 1.5 * Jishin$Mag) table(Jishin$震源地) Jishin$地域 <- ifelse(regexpr("愛知", Jishin$震源地) != -1, "関東", ifelse(regexpr("伊豆", Jishin$震源地) != -1, "関東", ifelse(regexpr("茨城", Jishin$震源地) != -1, "関東", ifelse(regexpr("奄美", Jishin$震源地) != -1, "関東", ifelse(regexpr("岩手", Jishin$震源地) != -1, "東北", ifelse(regexpr("岐阜", Jishin$震源地) != -1, "中部", ifelse(regexpr("宮城", Jishin$震源地) != -1, "東北", ifelse(regexpr("京都", Jishin$震源地) != -1, "関西", ifelse(regexpr("群馬", Jishin$震源地) != -1, "関東", ifelse(regexpr("三重", Jishin$震源地) != -1, "関西", ifelse(regexpr("三陸", Jishin$震源地) != -1, "東北", ifelse(regexpr("山形", Jishin$震源地) != -1, "東北", ifelse(regexpr("山梨", Jishin$震源地) != -1, "関東", ifelse(regexpr("滋賀", Jishin$震源地) != -1, "関西", ifelse(regexpr("秋田", Jishin$震源地) != -1, "東北", ifelse(regexpr("新潟", Jishin$震源地) != -1, "東北", ifelse(regexpr("神奈川", Jishin$震源地) != -1, "関東", ifelse(regexpr("青森", Jishin$震源地) != -1, "東北", ifelse(regexpr("静岡", Jishin$震源地) != -1, "関東", ifelse(regexpr("長野", Jishin$震源地) != -1, "中部", ifelse(regexpr("東京", Jishin$震源地) != -1, "関東", ifelse(regexpr("栃木", Jishin$震源地) != -1, "関東", ifelse(regexpr("福島", Jishin$震源地) != -1, "東北", ifelse(regexpr("兵庫", Jishin$震源地) != -1, "関西", ifelse(regexpr("千葉", Jishin$震源地) != -1, "関東", ifelse(regexpr("熊本", Jishin$震源地) != -1, "九州", ifelse(regexpr("宮崎", Jishin$震源地) != -1, "九州", ifelse(regexpr("広島", Jishin$震源地) != -1, "中国", ifelse(regexpr("島根", Jishin$震源地) != -1, "中国", ifelse(regexpr("トカラ", Jishin$震源地) != -1, "九州", ifelse(regexpr("千葉", Jishin$震源地) != -1, "関東", ifelse(regexpr("新島", Jishin$震源地) != -1, "関東", ifelse(regexpr("西表", Jishin$震源地) != -1, "九州", ifelse(regexpr("宮古", Jishin$震源地) != -1, "九州", ifelse(regexpr("和歌山", Jishin$震源地) != -1, "関西", "その他"))))))))))))))))))))))))))))))))))) table(Jishin$地域) Jishin[Jishin$地域 == "その他",] Chiiki <- names(table(Jishin$地域)) #---histogram hist(Jishin$Mag) hist(Jishin$Eng) library(gridExtra) plot(1, axes = F, xlab = "", ylab = "") grid.draw(tableGrob(table(Jishin$地域, trunc(Jishin$Mag)))) XLow <- 1298800000 XHigh <- 1300800000 plot(Jishin$Time, Jishin$Mag, type = "n", ylab = "マグニチュード", xlab = "日付", xlim = c(XLow, XHigh), ylim = c(0, 10), main = "観測値") for(i in 1:length(Chiiki)){ lines(Jishin$Time[Jishin$地域 == Chiiki[i, Jishin$Mag[Jishin$地域 == Chiiki[i, col = rainbow(length(Chiiki))[i]) points(Jishin$Time[Jishin$地域 == Chiiki[i, Jishin$Mag[Jishin$地域 == Chiiki[i, col = rainbow(length(Chiiki))[i], pch = i, cex = 0.7) } abline(v = as.numeric(Jishin$Time[Jishin$Mag == max(Jishin$Mag)]), lty = 2) abline(h = c(0, 5, 6), lty = 3) text(as.numeric(Jishin$Time[Jishin$Mag == max(Jishin$Mag)]) + 30000, 9.5, Jishin$Time[Jishin$Mag == max(Jishin$Mag)], adj = 0) text(as.numeric(Jishin$Time[Jishin$Mag == max(Jishin$Mag)]) + 30000, 9, "M9.0", adj = 0) legend("topleft", Chiiki, col = rainbow(length(Chiiki)), lty = 1, pch = 1:7) Spline <- smooth.spline(Jishin$Time, Jishin$Mag) plot(predict(Spline, seq(XLow, XHigh, len = 100)), type = "n", ylim = c(0, 10), xlab = "時間", ylab = "マグニチュード", main = "スプライン", xaxt = "n") for(i in c(1:4, 6, 7)){ Spline <- smooth.spline(Jishin$Time[Jishin$地域 == Chiiki[i, Jishin$Mag[Jishin$地域 == Chiiki[i) lines(predict(Spline, seq(XLow, XHigh, len = 1000)), col = rainbow(length(Chiiki))[i]) points(Jishin$Time[Jishin$地域 == Chiiki[i, Jishin$Mag[Jishin$地域 == Chiiki[i, col = rainbow(length(Chiiki))[i], pch = i, cex = 0.7) } abline(v = as.numeric(Jishin$Time[Jishin$Mag == max(Jishin$Mag)]), lty = 2) abline(h = c(0, 5, 6), lty = 3) text(as.numeric(Jishin$Time[Jishin$Mag == max(Jishin$Mag)]) + 30000, 9.5, Jishin$Time[Jishin$Mag == max(Jishin$Mag)], adj = 0) text(as.numeric(Jishin$Time[Jishin$Mag == max(Jishin$Mag)]) + 30000, 9, "M9.0", adj = 0) legend("topleft", Chiiki, col = rainbow(length(Chiiki)), lty = 1, pch = 1:7)
