福島原発のγ線量測定データをggplot2で可視化してみた

東京電力が福島原発周辺におけるγ線量などの測定データをPDFで随時公開しているが，これを奥村先生がCSV形式に加工して配布なさっている．

http://oku.edu.mie-u.ac.jp/~okumura/stat/data/

このデータをRのggplot2ライブラリで可視化してみた．

CSVデータの特徴

各CSVファイルは以下のような形式になっている．

$ head fukushima1.csv
【別紙】福島第一原子力発電所モニタリングカーによる計測状況,,,,,,,,,
計測日,計測時間,計測場所,γ線,中性子線,風向,風速(m/s),,,
3月11日,午後5時30分,体育館付近,49nGy/h,-,-,-,,,
,午後5時40分,正門付近,56nGy/h,-,-,-,,,
,午後5時50分,管理棟,64nGy/h,-,-,-,,,
,午後6時45分,MP-6,56nGy/h,-,-,-,,,
,午後7時00分,MP-7,57nGy/h,-,-,-,,,
,午後7時10分,MP-5,55nGy/h,-,-,-,,,
,午後7時15分,MP-4,59nGy/h,-,-,-,,,
,午後7時20分,MP-3,59nGy/h,-,-,-,,,

線量のカラムには測定値の単位も含まれており，このままでは扱いづらい．また，データが欠損している時間帯もあるため，可視化の際にはそのような時間帯が一目で分かるようにすることが望ましい．

文字コード，改行コードの修正

fukushima*.csvの中身をlessやheadでちゃんと見られない場合は，ファイルの文字コードや改行コードがおかしい可能性がある．この場合はnkfを使って文字コードや改行コードを変換すれば大丈夫．

以下を実行すれば，文字コードを直して元のファイルを上書きすることができる．

$ nkf --overwrite fukushima1.csv

CSVファイル加工スクリプト (convert_fukushima.rb)

fukushima*.csvをRなどで扱いやすい形式に変換するためのRubyスクリプト(convert_fukushima.rb)を用意した．Gistはこちら．

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

require "time"

puts %w[datetime mon day hour min place gamma neutron wind_dir wind_speed].join("\t")

lines = ARGF.each_line.to_a[2..-1]
date = ""
lines.each do |line|
  list = line.split(",")

  # date
  date = list[0].empty? ? date : list[0]
  date =~ /(\d+)\D+(\d+)/
  mon = $1
  day = $2

  # time
  time = list[1].gsub(/午後(\d+)/){$1.to_i + 12}
  time =~ /(\d+)\D+(\d+)/
  hour = $1
  min = $2

  # datetime
  datetime = "2011/#{mon}/#{day} #{hour}:#{min}"

  # place
  place = list[2].gsub("付近", "")

  # gamma
  gamma = list[3].to_f * (list[3].include?("μ") ? 1 : 1.0/1000)

  # neutron
  neutron = list[4].to_f * (list[4].include?("μ") ? 1 : 1.0/1000)

  # wind_dir
  wind_dir = list[5]

  # wind_speed
  wind_speed = list[6].chomp

  puts [datetime, mon, day, hour, min, place, gamma, neutron, wind_dir, wind_speed].join("\t")
end

使い方は以下の通り．

$ ruby convert_fukushima.rb fukushima1.csv > fukushima1.tsv

Ruby 1.8の場合は-Kオプションで文字コードを指定する必要があるかもしれない．

$ ruby -Ku convert_fukushima.rb fukushima1.csv > fukushima1.tsv

生成されたファイルはTSV形式(タブ区切り)となる．各列はそれぞれ日時，月，日，時，分，地点名，γ線量，中性子線量，風向，風量となっている．各線量の単位はμSv/h．

$ head fukushima1.tsv
datetime        mon     day     hour    min     place   gamma   neutron wind_dir        wind_speed
2011/3/11 17:30 3       11      17      30      体育館  0.049   0.0     -       -
2011/3/11 17:40 3       11      17      40      正門    0.056   0.0     -       -
2011/3/11 17:50 3       11      17      50      管理棟  0.064   0.0     -       -
2011/3/11 18:45 3       11      18      45      MP-6    0.056   0.0     -       -
2011/3/11 19:00 3       11      19      00      MP-7    0.057   0.0     -       -
2011/3/11 19:10 3       11      19      10      MP-5    0.055   0.0     -       -
2011/3/11 19:15 3       11      19      15      MP-4    0.059   0.0     -       -
2011/3/11 19:20 3       11      19      20      MP-3    0.059   0.0     -       -
2011/3/11 19:52 3       11      19      52      MP-6    0.057   0.0     -       -

ggplot2でγ線量を可視化

加工したデータを使って各地点のγ線量をggplot2で可視化してみる．

> library(ggplot2)

TSV形式のファイルはRのread.delim関数で読み込むことができる．

> d <- read.delim("fukushima1.tsv")
> head(d)
         datetime mon day hour min  place gamma neutron wind_dir wind_speed
1 2011/3/11 17:30   3  11   17  30 体育館 0.049       0        -          -
2 2011/3/11 17:40   3  11   17  40   正門 0.056       0        -          -
3 2011/3/11 17:50   3  11   17  50 管理棟 0.064       0        -          -
4 2011/3/11 18:45   3  11   18  45   MP-6 0.056       0        -          -
5 2011/3/11 19:00   3  11   19   0   MP-7 0.057       0        -          -
6 2011/3/11 19:10   3  11   19  10   MP-5 0.055       0        -          -

あらかじめ日時の文字列をPOSIXctオブジェクトに変換しておく．これをやっておかないとqplot関数で日時データを使ったプロットがうまくいかないので注意．

d$datetime <- as.POSIXct(d$datetime)

また，ggplot2の日本語文字化け対策を行っておく(参考：Using CJK Fonts in R and ggplot2 – Hi!!)．

> quartzFonts(HiraMaru=quartzFont(rep("HiraMaruProN-W4", 4)))

正門付近のγ線量

まず，正門付近のγ線量のみをプロットする．ここではqplot関数を使ってみる．

x軸は経過時間(emin / 60で単位はhours)，y軸はγ線量(単位はμSv/h)
geom，stat，positionを設定し，ヒストグラムライクの垂線プロットを行う(plot関数のtype = "h"と同じ)
scaleレイヤーで日時のフォーマットを指定する
文字化け対策としてthemeレイヤーを追加する

> qplot(datetime, gamma, data = d[d$place == "正門",], color = I("blue"),
    geom = "bar", stat = "identity", position = "identity",
    xlab = "Date", ylab = "Gamma ray [μSv/h]") +
    scale_x_datetime(format = "%m/%d") +
    theme_grey(base_family="HiraMaru")

γ線量のピークが数ヶ所あり，特に最後のピークは飛び抜けて大きいことが分かる．また，何もプロットされていない時間帯は，その時間帯のデータがないことを示している．垂線プロットのおかげでデータが欠損している時間帯が一目で分かる．

測定地点ごとに色分け

次に，各測定地点のデータを色分けしてひとつのグラフで描画してみる．color(colour)パラメータを使えば地点ごとに色分けしたグラフを描画することができる．

> qplot(datetime, gamma, data = d, color = place,
    geom = "bar", stat = "identity", position = "identity",
    xlab = "Date", ylab = "Gamma ray [μSv/h]") +
    scale_x_datetime(format = "%m/%d") +
    theme_grey(base_family="HiraMaru")

測定地点ごとにグラフ分割

最後に，各測定地点を別々のグラフに分けて描画してみる．facetsパラメータを使えば地点ごとにグラフを分割して描画することができる．

> qplot(datetime, gamma, data = d, color = place, facets = ~ place,
    geom = "bar", stat = "identity", position = "identity",
    xlab = "Date", ylab = "Gamma ray [μSv/h]") +
    scale_x_datetime(format = "%d") +
    theme_grey(base_family="HiraMaru")