福島原発のデータ fukushima*.csv ( http://oku.edu.mie-u.ac.jp/~okumura/stat/data/ ) をRでグラフにするスクリプトです．使い方についてはコメント欄を参照してください．

demo316open.R

# -*- coding: utf-8 -*-

## 1)http://oku.edu.mie-u.ac.jp/~okumura/stat/data/
## で配布されている、fukushima1.csvをダウンロードします
## https://gist.github.com/869400
## で配布されているスクリプトでデータを整形します

## テキストエディタ等で次のファイルのようにデータを更に整形します
## https://gist.github.com/871960
## 1列目から 日付, 時刻, 観測箇所, 観測値 に対応しています
## 観測箇所に関しては、MP-1からMP-8を1から8、正門 11, 管理 12, 体育館 13, 西門 14 と数字で対応させます


## 初回実行時には以下のコメントアウトを外してください。
## ----------------- setting ----------------------
## パッケージのダウンロード先を変更します
## options(repos="http://cran.cnr.berkeley.edu/")
## zooライブラリをインストールします
## install.packages("zoo")
## -------------------------------------------------


## ライブラリ呼びこみ
library(zoo)

## データの読み込み
dat <- read.csv("fukushima316-mod.csv")

## 計測時間を24時間単位で算出します
day <- min(dat[,1]):max(dat[,1])
dayidx <- cbind(c(1,which(diff(dat[,1])==1)+1),
                c(which(diff(dat[,1])==1),nrow(dat)))
time1 <- ((dat[,2]%/%100)*60 +dat[,2]%%100)/(60*24) + dat[,1]

## 各計測点毎にデータを整頓
place <- c(1:8,11:14)
pidx <- vector("list",length(place))
ptime <- vector("list",length(place))
dat1 <- matrix(0,length(time1),length(place))
for(i in 1:length(place)){
  pidx[[i]] <- which(dat[,3]==place[i])
  ptime[[i]] <- ((dat[pidx[[i]],1]-11)*24*60 + dat[pidx[[i]],3])/(60*24) + 11
  dat1[pidx[[i]],i] <- dat[pidx[[i]],4]
}

## 全計測時刻毎にデータを整頓
level <- length(unique(time1)) 
ltime <- numeric(length(time1))
epoint <- numeric(level)
for(i in 1:level){
  flag <- TRUE
  if(i==1){
    spoint <- 1
  }else if(i==level){
    epoint[i] <- spoint + 1
    next
  }else{
    spoint <- epoint[i-1]+1
  }
  tmp <- time1[1]
  while(flag){
    if(time1[spoint] == time1[spoint + 1]){
      spoint <- spoint + 1
    }else{
      epoint[i] <- spoint
      flag <- FALSE
    }
  }
}
ltime <- time1[epoint]

dat2 <- matrix(0,level,length(place))
for(i in 1:level){
  if(i == level){
    J <- length(epoint[i]:level)
    Jidx <- epoint[i]:level
  }else{
    J <- length(epoint[i]:(epoint[i+1]-1))
    Jidx <- epoint[i]:(epoint[i+1]-1)
  }
  for(j in 1:J){
    dat2[i,which(dat[Jidx[j],3] == place)] <- dat[Jidx[j],4]
  }
}


## zoo形式のデータを作ります
## x: 計測値, order.by: 計測時刻
dat3 <- cbind(dat2[,9:12],dat2[,1:8])
tmp <- zoo(x=dat3,order.by=ltime)
## グラフィックのパラメータを調整します
par("bg" = "grey80")
leg <- c("Front","Kanri","Gym","West",
         "MP-1","MP-2","MP-3","MP-4","MP-5","MP-6","MP-7","MP-8")

## プロットします
## type="h"をtype="o"に変更すると、計測点間を直線で繋ぎます
## ylimでy軸の範囲を設定しています
ymax <- 13000
matplot(x=ltime,y=dat3,type="h",lty=1,lwd=1.5,col=rainbow(length(leg)),
        ylim=c(0,ymax),
        xlab="Time[day], from '3/11/2011 17:30' to '3/15/2011 23:35'",
        ylab="Gamma ray [uSv/h]",main="Fukushima nuclear power plant, No.1")

## 補助線を描きます
abline(v=day,lwd=2,col="grey95",lty=1)
abline(v=day+0.5,lwd=1,col="grey90",lty=2)
abline(h=seq(1000,ymax,by=2000),lwd=2,col="grey95",lty=1)
abline(h=seq(1000,ymax,by=1000),lwd=1,col="grey90",lty=2)

## 日にちを図上に記入します
for(i in 1:length(day)){
  text(x=day[i]+0.2,y=ymax,labels=paste("3/",day[i],sep=""),cex=2)
}
## 凡例を描きます
legend(x=12,y=12000,legend=leg,col=rainbow(ncol(dat3)),lty=1,lwd=4)

## 画像を保存します
##dev.copy2eps(file="0044h.eps")
##dev.copy2pdf(file="004h.pdf")