Bonjour à tous,
je travaille depuis quelques jours sur un algorithme basique, très répandu dans la littérature d'agronomie appliquée, qui vise à calculer le nombre de degré-jour au-dessus d'un certain seuil de température pour déterminer le moment d'émergence d'une plante donnée.
Comme illustré sur ce lien: http://www.oardc.ohio-state.edu/gdd/glossary.htm , j'ai crée ma fonction dans R comme suit:
Or si je crée des données aléatoires:
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30 DD <- function(Tmin, Tmax, Tseuil, meanT, method = "DDsin") ### function that calculates the degree-days based on ### minimum and maximum recorded temperatures and the ### minimal threshold temperature (lower growth temperature) { ### method arcsin if(method == "DDsin"){ cond1 <- (Tmax <= Tseuil) cond2 <- (Tmin >= Tseuil) amp <- ((Tmax - Tmin) / 2) print(amp) print((Tseuil-meanT)/amp) alpha <- asin((Tseuil - meanT) / amp) print(alpha) DD_ifelse3 <- ((1 / pi) * ((meanT - Tseuil) * ((pi/2) - alpha)) + amp*cos(alpha)) DD <- ifelse(cond1, 0, ifelse(cond2, (meanT - Tseuil), DD_ifelse3)) } ### method (Tmin + Tmax) / 2 else if(method == "DDt2"){ cond1 <- (meanT > Tseuil) DD <- ifelse(cond1,(meanT - Tseuil),0) } else{ stop("\nMethod name is invalid.\nMethods available = DDsin (sinus) or DDt2 (mean)\n") } return(DD) }
J'obtiens un warning sur le calcul d'alpha (production de NaN); en effet, les valeurs que je donne en argument à asin() sont en-dehors du range [-1:1]. Je comprends pas comment régler ce problème, car le même algorithme sur Excel (visual basic) est fonctionnel.library(reshape2)
library(plyr)
station <- rep(c("station1","station2","station3"), 20)
values <- sample(-5:20, size = 60, replace = T)
values2<- sample(20:40, size = 60, replace = T)
meanT <- ((values+values2)/2)
d <- data.frame(station,values,values2,meanT)
names(d) <- c("station", "values_min","values_max","meanT")
x<-ddply(d, .(station), transform, t1 = cumsum(DD(values_min,values_max,0,meanT)))
Je comprends vraiment pas comment régler ce contre-temps, il semblerait pourtant que l'algorithme est très courant, mais je n'ai vu aucune signalisation de cette erreur sur internet.
Toute aide est bienvenue merci
Partager