ホーム | 統計 Top | 害虫発生予察のための各日有効温量をRで計算する Estimating daily effective temperature for a pest-occurrence forecast with R.

（まだ執筆中です）

古典的な有効温量モデルに加え、SSIモデルに関する日本語情報の集約を行う予定。時期は......今投稿している論文がアクセプトされたら早いうちに。

昆虫やダニ等の変温動物において、卵や幼虫といった各発育ステージの所要時間は気温の関数で（おおよそ）表すことができる（総説：桐谷 2012）。

【自作関数】三角法による有効積算温度（日度）を算出する関数【しょぼい】

以下の行をRコンソールにコピペしてエンターキー、で関数が定義される。すでに sankaku という名称のオブジェクトをR内で使用している場合は、
sankaku <- function(......の部分を適宜変更する。

##三角法による有効積算温度（日度）を算出する自作関数sankaku
## 坂神・是永（1981）日本応用動物昆虫学会誌 25 (1): 52–54. を基にシコシコ計算
## Copyright by Sudo Masaaki (NIAES) 
# (2011.12.16初版)
# (2014.04.05改訂版) 初版において最低気温がK1以上K2未満のときの係数が間違っていたので修正。
## n1, m=n1, n2=n1, K1=13.6, K2=30, K3=35を引数の初期値として定義
## 各引数は以下の通り
## n1 = (当日最低気温),
## m  = (当日最高気温、指定がなければ当日最低気温を代入),
## n2 = (翌日最低気温、指定がなければ当日最低気温を代入),
## K1 = (発育零点、指定がなければ摂氏13.6度),
## K2 = (発育上限温度、指定がなければ摂氏30.0度),
## K3 = (発育停止温度、指定がなければ摂氏35.0度)

sankaku <- function(n1,m=n1,n2=n1,K1=13.6,K2=30,K3=35) {
    if (is.na(prod(n1,m,n2, K1,K2,K3))==TRUE) return(NA) else
    if (n1 > m) return(NA) else
    if (n2 > m) return(NA) else
#######内部データを定義########
i1a <- 8*(K1-n1)/(m-n1) ##午前にK1となる時刻
i2a <- 8*(K2-n1)/(m-n1) ##午前にK2となる時刻
i3a <- 8*(K3-n1)/(m-n1) ##午前にK3となる時刻
i1b <- 16*(K1-n2)/(m-n2) ##午後にK1となる時刻
i2b <- 16*(K2-n2)/(m-n2) ##午後にK2となる時刻
i3b <- 16*(K3-n2)/(m-n2) ##午後にK3となる時刻
#######上昇期######## n1からmまで、8時間かけて上昇
morning <- if (m < K1) { # 1) m < K1
    0    ##最高気温が発育零点以下のとき、もちろん有効積算温量はゼロ。 
} else {
    if (m < K2) { # 2) K1 < m < K2
        if (n1 < K1) {
            (8-i1a)*(m-K1)/2
        } else {
            4*(m+n1-2*K1) # 最初の時点ですでに有効温度(n1-K1)に達している
        } 
    } else { # 3) K2 < m < K3
        if (m < K3) { # 最高気温がK2以上 K3未満
            if (n1 < K1) { # 最低気温がK1以下
                (K2-K1)*(i2a-i1a)/2+(K2-K1)*(8-i2a) 
                # 時刻i1aからi2aまでかけて0からK2-K1まで上昇。
            } else { 
                if (n1 < K2) { # 最低気温がK1以上K2未満
                    (K2+n1-2*K1)*i2a/2 + (K2-K1)*(8-i2a) 
                    # 時刻0からi2aまでかけてn1-K1からK2-K1まで上昇。前半部分の係数を修正
                } else { # 最低気温がK2以上K3未満、よって全区間定温で積分
                    (K2-K1)*8
                }
            }
        } else { # 4) K3 < m
            if (n1 < K1) { # 最低気温がK1以下
                (K2-K1)*(i2a-i1a)/2 + (K2-K1)*(i3a-i2a)
            } else {
                if (n1 < K2) { # 最低気温がK1以上K2未満
                    (K2+n1-2*K1)*i2a/2 + (K2-K1)*(i3a-i2a)
                    # 時刻0からi2aまでかけてn1-K1からK2-K1まで上昇。前半部分の係数を修正
                } else {
                    if (n1 < K3) { # 最低気温がK2以上 K3未満
                        (K2-K1)*i3a 
                    } else { # 最低気温がK3以上
                        0
                    }
                }
            }
        }
    }
}
#######下降期######## n2からmまで、16時間かけて上昇（実世界では「mからn2まで、16時間かけて下降」なのだが計算の便宜上時間を逆転させている。詳しくは元文献を参照）
afternoon <- if (m < K1) { # 1) m < K1
    0    ##最高気温が発育零点以下のとき、もちろん有効積算温量はゼロ。 
} else {
    if (m < K2) { # 2) K1 < m < K2
        if (n2 < K1) {
            (16-i1b)*(m-K1)/2
        } else {
            8*(m+n2-2*K1) # 最初の時点ですでに有効温度(n2-K1)に達している
        } 
    } else { # 3) K2 < m < K3
        if (m < K3) { # 最高気温がK2以上 K3未満
            if (n2 < K1) { # 最低気温がK1以下
                (K2-K1)*(i2b-i1b)/2+(K2-K1)*(16-i2b) 
                # 時刻i1bからi2bまでかけて0からK2-K1まで上昇。
            } else { 
                if (n2 < K2) { # 最低気温がK1以上K2未満
                    (K2+n2-2*K1)*i2b/2 + (K2-K1)*(16-i2b) 
                    # 時刻0からi2bまでかけてn1-K1からK2-K1まで上昇。前半部分の係数を修正
                } else { # 最低気温がK2以上K3未満、よって全区間定温で積分
                    (K2-K1)*16
                }
            }
        } else { # 4) K3 < m
            if (n2 < K1) { # 最低気温がK1以下
                (K2-K1)*(i2b-i1b)/2 + (K2-K1)*(i3b-i2b)
            } else {
                if (n2 < K2) { # 最低気温がK1以上K2未満
                    (K2+n2-2*K1)*i2b/2 + (K2-K1)*(i3b-i2b)
                    # 時刻0からi2bまでかけてn2-K1からK2-K1まで上昇。前半部分の係数を修正
                } else {
                    if (n2 < K3) { # 最低気温がK2以上 K3未満
                        (K2-K1)*i3b 
                    } else { # 最低気温がK3以上
                        0
                    }
                }
            }
        }
    }
}
#######戻り値を合成########
return((morning+afternoon)/24 )
}
#######ここまで########

今のところ公式なパッケージ等への登録を目指すつもりはないが、バグ等にお気付きの際は教えていただければ幸いである。

Reference

• Kiritani K (2012) The low development threshold temperature and the thermal constant in insects and mites in Japan (2nd edition). Bulletin of National Institute for Agro-Environmental Sciences 31: 1–74 (in Japanese)
• Sakagami Y, Korenaga R (1981) "Triangle Method", a simple method for the estimation of total effective temperature. Jpn J Appl Entomol Zool 25: 52–54 (in Japanese)