日本全国の地震の推移をggplot2で可視化する

RubyやRでスクリプトを書いて可視化してみた.

tenki.jpの地震情報を取得するRubyスクリプト (fetch_quake_info.rb)

以下は2011/3/11 0:00から現在時刻までのtenki.jp地震情報を取得するRubyスクリプト.Nokogiriを使ってスクレイピングする.

#!/usr/bin/env ruby
# -*- coding: utf-8 -*-

require "open-uri"
require "nokogiri"

puts "datetime\tplace\tmagnitude"

i = 1
cont = true
stop = Time.parse(ARGV[0] || "2011/03/11 0:00")

while cont
  doc = Nokogiri::HTML(URI("http://tenki.jp/earthquake/entries?p=#{i}").read)
  first = true
  doc.search("table#seismicInfoEntries/tr").each do |row|
    # skip the first row
    if first
      first = false
      next
    end

    entry = row.search("td").map{|e| e.text}
    entry[0].gsub!("", "/"); entry[0].gsub!("", "/"); entry[0].gsub!(/.+/, "")
    entry[1].gsub!("", ":"); entry[1].delete!("")
    entry[2].gsub!("", ":"); entry[2].delete!("")

    # datetime
    datetime = Time.parse("#{entry[0]} #{entry[2]}")
    if datetime < stop
      cont = false
      break
    end

    # place
    place = entry[3].gsub("---", "不明")

    # magnitude
    magnitude = entry[4][0] == ?M ? entry[4][1..-1].to_f : nil

    puts [datetime, place, magnitude].join("\t")
  end
  i += 1
end

以下のように実行すれば,2011/3/11 0:00から現在時刻までの地震情報を取得することができる.Ruby 1.9で動作確認済み.

$ ruby fetch_quake_info.rb > quake.tsv
$ head quake.tsv
datetime	place	magnitude
2011-03-19 19:41:00 +0900	茨城県南部	3.0
2011-03-19 18:57:00 +0900	不明	
2011-03-19 18:56:00 +0900	茨城県北部	6.1
2011-03-19 16:24:00 +0900	長野県中部	2.9
2011-03-19 16:14:00 +0900	長野県北部	2.3
2011-03-19 13:07:00 +0900	新潟県上越地方	3.3
2011-03-19 12:58:00 +0900	長野県北部	2.3
2011-03-19 12:08:00 +0900	福島県会津	2.9
2011-03-19 10:59:00 +0900	新潟県中越地方	2.0

なお,スクリプトを実行するとtenki.jpへのアクセスが発生するため,何度も実行して負荷をかけ過ぎたりしないように注意する必要がある.

Rのggplot2ライブラリによる可視化

Rを起動してデータの読み込みと整形を行う.

library(ggplot2)

# タブ区切りデータなのでread.delim関数を使う
d <- read.delim("quake.tsv")

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

# 日時の文字列をPOSIXctオブジェクトに変換する
d$datetime <- as.POSIXct(d$datetime)
地震の発生日時とマグニチュードの可視化

qplot関数をgeom = "bar", stat = "identity", position = "identity"で呼び出し,ヒストグラムライクな垂線プロットで可視化する.地震の発生時刻,マグニチュードの大きさ,発生頻度がひと目で分かって便利.

qplot(datetime, magnitude, data = d, color = I("red"),
  main = "Earthquakes in Japan (3/11 0:00 ~ 3/19 20:00)",
  geom = "bar", stat = "identity", position = "identity") +
  scale_x_datetime(format = "%m/%d")

震源地の情報を無視しているため,余震以外の地震もプロットされているので注意.

各日の地震のマグニチュードと回数の可視化

facetsパラメータを使えば,d$dateの値ごとにグラフを作ることができる(条件付きプロット).

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

# 各日のヒストグラムをプロットする
qplot(magnitude, data = d, color = I("red"),
  main = "Earthquakes in Japan (3/11 0:00 ~ 3/19 20:00)",
  facets = ~ date, binwidth = 0.5)

ここ数日の地震はそれほど減ってないように見える.

以上

Rubyによるスクレイピングやggplot2による可視化の例として参考になれば幸いです.