Rで放射線モニタリング情報を時系列プロットする

データを公開するときはCSV形式みたいに、RAWデータかそれに近い形で提供してくれると、いろいろと利用できて嬉しいですよね。

原子力規制委員会のサイトでは、モニタリングポストで測定された空間線量のデータをCSV形式で公開しているので、これを使って時系列データの可視化を試してみましょう。

下記のページで、都道府県やエリア、期間の日付を選択してダウンロードできます。

ダウンロード | 東日本大震災関連情報 放射線モニタリング測定結果等 | 原子力規制委員会

試しに東京都の都健康安全研究センターのモニタリングポストのデータを2014年1月1日～2014年12月31日の1年分ダウンロードしてみました。10分ごとのデータなのでかなりの量ですね。約5.5MBくらいありました。

まずは、読み込んで、冒頭部分を表示。

> d <- read.csv("新宿2014.csv", header=F, as.is=T)
> head(d)
          V1                             V2       V3       V4
1 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987
2 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987
3 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987
4 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987
5 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987
6 東京都全域 新宿区　都健康安全研究センター 35.70608 139.6987
                V5    V6     V7    V8     V9 V10 V11
1 2014/12/31 23:50 0.034 μSv/h 0.061 μSv/h  ‐  NA
2 2014/12/31 23:40 0.033 μSv/h 0.059 μSv/h  ‐  NA
3 2014/12/31 23:30 0.033 μSv/h 0.059 μSv/h  ‐  NA
4 2014/12/31 23:20 0.033 μSv/h 0.059 μSv/h  ‐  NA
5 2014/12/31 23:10 0.033 μSv/h 0.059 μSv/h  ‐  NA
6 2014/12/31 23:00 0.033 μSv/h 0.059 μSv/h  ‐  NA

V5が測定日時時刻、V6が線量の測定値です。（これだけじゃ分かりませんが、V8は高さ1mで換算した推計値らしいです）

V5を横軸、V6を縦軸にすれば、線量の時系列変化が確認できますね。

ただそのままではV5は文字列データなので、日付時刻型に変換する必要がありますね。詳細は以前の記事↓に書きましたが、

Rで文字列を日付時刻型に変換する - Rプログラミングの小ネタ


このデータの場合は as.POSIXlt関数１つで簡単に変換できます。

plot(as.POSIXct(d$V5), d$V8, type="l")

放射線量測定データを時系列プロット（1年分）

さすがに10分ごとに1年分だとデータが多いので、ひと月ぶんだけ取り出してプロットしてみます。

n <- 6 * 24 * 31 # 1ヶ月ぶんの10分の個数
plot(as.POSIXct(d$V5[1:n]), d$V8[1:n], type="l")

放射線量測定データを時系列プロット（ひと月分）

上記は12月ぶんのプロットになっていますね。

細かい変動には意味がないと思われますので、移動平均を取ってみましょう。やり方は過去記事↓にも書きました。

Rで移動平均を求める - Rプログラミングの小ネタ

# 何度か使うので関数化しておく
moving_average <- function(x, n){
  filter(x, rep(1,n)) / n
}
 
時刻   <- as.POSIXct(d_part$V5[1:n])
測定値 <- d_part$V8[1:n]
plot(時刻, moving_average(測定値,3), type="l")

10分前、現時刻、10分後で移動平均を取ったもの

多少見やすくなりましたね。移動平均を取っていないものと、移動平均を取る期間の幅を変えたものを1つのグラフに重ねてプロットし、比較してみます↓

plot( 時刻,                測定値,      type="l", col="gray")
lines(時刻, moving_average(測定値,  6), type="l", col="green")
lines(時刻, moving_average(測定値,144), type="l", col="blue")
legend( "topright",
        lty=c(1,1,1),
        col=c("gray","green","blue"),
        c("測定データそのまま", "1時間の幅で移動平均", "1日の幅で移動平均") )

移動平均なし（グレー）、1時間幅（緑）、1日幅（青）

1日幅の移動平均（青の線）だと、だいぶ傾向が見やすくなったんじゃないでしょうか。山の現れ方には周期があるように見えなくもないです。もし本当に周期性があるなら、それが現れる原因は何か、などを探ってみると面白いかもしれません。