Rで、時系列データを扱うときは、POSIX形式の時間データとして扱わないといけないので、ほかの処理系に比べて、少し、癖があります。
で、V4368 Sgrの光度曲線を描きました。
本当は、横軸をJDにすれば、こんなに、苦労しないのですが、必要にかられて、日付にしました。
data <- read.table('V4368Sgr.std') # stdデータの読み込み date_std <- data$V2 mag <- as.numeric(data$V3) band <- data$V4 date <- strptime(date_std, "%Y%m%d%H%M") # 年月日時間として変換 # グラフ table <- data.frame(date, mag , band) B_table <- subset(table, band=="B", c(date, mag)) # バンドごとに抽出 V_table <- subset(table, band=="V", c(date, mag)) Rc_table <- subset(table, band=="Rc", c(date, mag)) Ic_table <- subset(table, band=="Ic", c(date, mag)) y_table <- subset(table, band=="y", c(date, mag)) r <- as.POSIXct(round(range(date), "days")) # 横軸の範囲を決める endtime <- as.POSIXct("2015-03-22 00:00:00") #plot(V_table, xlim=c(r[1],r[2]), ylim=c(13, 9), xaxt="n", col="green", pch=1) # バンドごとにプロット plot(V_table, xlim=c(r[1], endtime), ylim=c(13, 9), xaxt="n", col="green", pch=19) # バンドごとにプロット points(B_table, col="blue", pch=2) points(Rc_table, col="red", pch=5) points(Ic_table, col="black", pch=6) points(y_table, col="yellow", pch=3) axis.POSIXct(1, at=seq(r[1],endtime, by="1825 days"), format="%y/%m/%d") # 横軸を書く
もとのコードは、Iakさんに書かれたものです。
最初から、UTのstdファイルにしたり、測光バンドはあらかじめ置換して、スペースを入れて分離しておいて、少し、手間を減らしてします。
#あれ、Icのばらつきが大きいですね。元データを確認しないと。
(追記)
attach(data)
の行は、わたしが、いつもの流儀でdata$V2の様にframe内の要素を参照するように書き換えた段階で要らなくなったと、Katさんに教えていただいたので、削除しました。
