2008年11月21日

時系列データのプロット[R] (「おんどとり」を例にした解説)

Rには日付・時間データをオブジェクトとして扱うための"POSIXlt"、"POSIXct"というクラスがある。例えば、Sys.time()関数を使うと現在の時間が"POSIXct"オブジェクトとして返ってくる。

> Sys.time()
[1] "2008-11-21 18:39:48 JST"

日時の情報は例えばデータロガーによって温度、湿度、二酸化炭素濃度などの環境要素の経時変化を観察するような場合に重要である。そこで、データロガーの出力する日時の情報をPOSIXctオブジェクトにでも変換して扱えればいいのだが、そのためには多少の作業が必要となる。

「おんどとり」というデータロガーがある。比較的安価で(といっても1台3万程度だが)温度・湿度情報を採取できるために研究・一般を問わず広く使われているデータロガーだ。「おんどとり」にはいくつか種類があるが、その中のある機種が吐き出すデータ(を付属ソフトで「カンマ区切りテキスト」に変換したもの)は次のような形式になっている。

"< Measurement Data Text Output >"
""
"Date/Time","@date()","ch.1","ch.2","ch.3","ch.4","ch.5","ch.6","ch.7","ch.8"
"","","゚C","゚C","","","","","",""
"2008/10/30 16:00'00","39751.6666666667","  20.8000","  20.4000","","","","","","",
"2008/10/30 16:20'00","39751.6805555556","  20.7000","  20.0000","","","","","","",
"2008/10/30 16:40'00","39751.6944444444","  20.4000","  19.4000","","","","","","",
"2008/10/30 17:00'00","39751.7083333333","  20.2000","  19.0000","","","","","","",
"2008/10/30 17:20'00","39751.7222222222","  19.9000","  18.6000","","","","","","",
"2008/10/30 17:40'00","39751.7361111111","  19.6000","  18.3000","","","","","","",
.
.
.

見づらいが、Excelで開くと見やすくなる。つまりExcel用である。

Rには直接関係無いが説明のためにちょっと詳しく見よう。

まず行を見ると、最初に余分な2行がついていて、3行目にデータのラベルが記述されている。そして4行目に採取したデータの単位が記述されており、その後に続くのがデータだ。

次に列を見る。左端の1列目は日時情報を可読性重視で出力している。そして2列目にある大きな数字だが、これはシリアル値と呼ばれるWindows特有の日時表現形式で、1900年1月1日0時を1としたときの経過日数となっている。ここで面倒なのがUNIX界(POSIX規約も同様)では同様に1970年1月1日を0とし、統計ソフトSASでは1960年1月1日を0とするような値によって日時を表現するということだ。要するにこのシリアル値というのはWindows専用、もっといえばExcel専用というわけ。実際、Excelで表示したときに書式を適当なものへ変更してやればこの値が正確に記録時刻を表す。が、とりあえずRでは使えない。使えないことも無いが、面倒なので使わない。で、3列目移行がデータとなっている。

さて、まずはデータを読み取ろう。データはRの作業ディレクトリに「data.csv」という名前で保存されているものとする。

mydata <- read.table("data.csv", skip=4, sep=",", header=F)
names(mydata) <- c("Date.Time", "d", "ch.1", "ch.2", "ch.3", "ch.4", "ch.5", "ch.6", "ch.7", "ch.8")

ここでは頭4行を飛ばして読んだあとに改めてラベルをつけなおしていることに注意。特に頭2行を読み込むと列数が合わなくてややこしいことになるので必ず飛ばさないといけない。

次に1列目のデータを利用してPOSIXltオブジェクトを作成する(読み込んだ時点では文字データになっている)。

attach(mydata)
x <- strptime(Date.Time, "%Y/%m/%d %H:%M'%S", tz="")

とりあえずattach()してデータを扱いやすく。

次にstrptime()関数が文字列をPOSIXltオブジェクトへ変換してくれる。一つ目の引数に文字列ベクトルを、2つ目の引数に文字列の書式を、3つ目にはタイムゾーン(普通は空で問題ない)を入れている。2つ目の書式設定が重要。「おんどとり」では"2008/10/30 17:40'00"といった具合なのでそれにあわせて"%Y/%m/%d %H:%M'%S"となる。もしこれが"2008年10月30日17時40分00秒"だったら、"%Y年%m月%d日%H時%M分%S秒"とやるわけだ。ちなみに、もし年の表記が2桁なら%yと小文字にする。

そうしたらこれを使ってプロットをしてみよう。今日の目標は時系列データをプロットすることだ。

普通にplot(x, ch.1)とかやってもいいけど、きっとx軸のラベルが気に入らないと思うので、最初プロットするときは空にしておいて後でPOSIXオブジェクト専用の軸を描く関数、axis.POSIXct()関数で追加するのがオススメ。とりあえず軸なしでグラフを書き上げる。

plot(x, ch.1, type="l", ylim=c(5,25),
	 xlab="Date", ylab="Air Temperature (℃)", xaxt="n")
lines(x, ch.2, col=2)
legend(locator(1), c("ch.1", "ch.2"), col=1:2, lty=c(1,1), box.lty=0)

X座標の値と位置の関係が分かりにくいので凡例の位置は上のようにlocator()で対話的に決めたほうが良いような気がする。他のデバイスへ出力したいときはあらかじめlocator()で調べておくしかない?一応今回使ったデータで適当にlocator()を実行した結果はこんな感じ。

> a <- locator()
> a
$x
[1] 1225740207 1226309024 1226763384

$y
[1]  9.146104  8.647322 24.699051

x軸のこの大きな値は何だろう?秒数?

で、次に軸を描く。ここが重要。といってもまあaxis.POSIXct()関数のヘルプをそのまんま和訳したようなもんですが。

r <- as.POSIXct(round(range(x), "days"))
axis.POSIXct(1, at=seq(r[1],r[2], by="1 week"), format="%d日 0時")

最初にrというオブジェクトに入れているのはデータの「端」を丸めたもの。ここで"days"に丸めるとちょうど0時になる。"hours"なら0分、"secs"なら0秒といった具合になる。これはプロットするときに最初の目盛りが入る位置になるので、データの長さや目的に応じて換えるといい。

そしてaxis.POSIXct()だ。状況に応じて変更すべきは真ん中のby="1 week"とformat="%d日 0時"の部分。by=""には目盛りを刻む間隔を入れる。sec, min, day, weekなどの単位をつけて時間の長さで指定する。先に"days"で丸めたからといって"1 day"より長くする必要は無い。format=""には目盛りのラベルをどのような書式で書くかを指定する。strptime()のときに使った%Y, %m, %dなどが使える。ここでは0時と直接書いたが、%Hでもいい。ただし「00時」という表記になる。

では完成品を。

datetimeplot.jpg

なんだかずいぶん長くなってしまったけれど作業量自体はたいしたことが無いので、あるいはExcelより楽かもしれない。それにExcelでは扱えないような量のデータも扱えるはず。それに各種の時系列解析も使えるハズ(やったことないけど)。

まあ見た目の問題はあるっちゃあるけど。プレゼンには向かない?

posted by Rion778 at 18:35 | Comment(1) | TrackBack(0) | PC関連。HTMLとか,Linuxとか,Rとか | このブログの読者になる | 更新情報をチェックする
この記事へのコメント
ここで紹介されているコマンドを参考にさせていただきました。
ずっと懸案だったこと(RでX軸が日付けのグラフ作成)ができたので、ひとりで感動してます。ありがとうございます!!!

試行錯誤したところ、attachだと、detachしても、別のファイルを読み込むときに同じラベルを使うとエラーが出てしまいました・・・withで囲ってやるほうがいい感じみたいです。
Posted by Nic at 2011年06月02日 19:51
コメントを書く
お名前: [必須入力]

メールアドレス: [必須入力]

ホームページアドレス: [必須入力]

コメント: [必須入力]

認証コード: [必須入力]


※画像の中の文字を半角で入力してください。
この記事へのトラックバックURL(言及がない、関連がないと判断した場合削除することがあります)
http://blog.seesaa.jp/tb/110002044
※言及リンクのないトラックバックは受信されません。

この記事へのトラックバック
×

この広告は1年以上新しい記事の投稿がないブログに表示されております。