ホーム | 統計 Top | 害虫発生予察のための各日有効温量をRで計算する Estimating daily effective temperature for a pest-occurrence forecast with R.
2015/04/04追記:個人的に使うため再チェックしてみたところ、以前に公開していたバージョンにバグ(最低気温がK1以上K2未満のときの係数の間違い)があったため修正します。
(まだ執筆中です)
古典的な有効温量モデルに加え、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)