| Top Page | プログラミング | R 自動化 目次 | 8章へ |

#  渡された温度データのベクトルを，最初からひとつづつ積算して積算温度のベクトル
#  を作製する関数．積算温度を計算するときの基準温度 threshold はディフォルトでは
#  5 度となる．

calc.deg.day <- function(t.data, threshold = 5) 
{
    x <- t.data[1] - threshold  #  1 日めの積算温度
    if (x < 0) {
        x <- 0
    }

    deg.day = c(x)  #  要素がひとつのベクトル

    for (i in 2:length(t.data)) {     #  2日めからの温度
        if (t.data[i] > threshold) {  #  基準温度より高ければ…
            x <-  (t.data[i] - threshold) + deg.day[i - 1]  # 前日までの積算温度に加える．
        }
        else {
            x <- deg.day[i - 1]        #  そうでなければ，積算温度は前日の値のまま．
        }
        deg.day <- c(deg.day, x)     #  積算温度のベクトルの後尾に x を付け加える．
    }

    return (deg.day)     #  t.data の最後まで行ったあと，できたベクトルを返す．
}


########################################################
# calc.deg.day を使って積算温度を計算し，グラフを描く

t <- read.table ('temperature.txt', header = TRUE)  # データの読み込み

n.sites <- ncol(t) - 1    # 測定地点数

for (i in 1:n.sites) {              # 各地点について…
    calc.deg.day(t[[i + 1]]) -> dd  #  積算温度を計算する（温度データは２列目から）
    t <- cbind(t, dd)               #  データフレーム t に列を加える．
    colnames(t)[ncol(t)] <- sprintf("DegDay%02d", i)  #  列に名前をつける．
}

# グラフの座標系だけつくる．
plot(0, 0, type = 'n',              #  type = 'n' なので，渡すデータはなんでもよい．
     xlim = range(t[1]),                        #  横軸の範囲
     ylim = range(t[(n.sites + 2):ncol(t)]),    #  縦軸の範囲
     xlab = 'Day', ylab = 'Degree Day')

for (i in (n.sites + 2):ncol(t)) {     #  ２列めから最後の列まで
    lines(t[[1]], t[[i]], col = i - n.sites - 1)  #  日付を横軸に，温度を縦軸に線グラフ
}

#  ↑ 上で， (n.sites + 2):ncol(t) としていることに注意．
#  最初のカッコをなしに n.sites + 2:ncol(t) とすると，
#  2:ncol(t) で作られるすべての要素に n.sites が足されたベクトルになる．

par(usr = c(0, 1, 0, 1))    # 凡例が書きやすいように座標系を設定する．

#　色の設定は lines での設定と対応させることに注意．
legend(0.1, 0.9, colnames(t[(n.sites + 2):ncol(t)]), col = 1:n.sites, lty = 1)