Skip to content
GitLab
Select Git revision
  • master default
  • BRIEFCOM
  • ctfs
  • test.dplyr
  • package
Find file BlameHistoryPermalink
FUN.climate.R 4.61 KiB
1 
#####################################
2 
#####################################
3 
#### FUNCTION TO COMPUTE CLIMATIC DATA FOR FRANCE
4
5 
## function to compute sum of degree days above 5.56
6 
##' .. Compute teh sum of degree days above 5.56C from monthly data with a spline ..
7 
##'
8 
##' .. A spline is used to smooth the variation of T over each day from mothly data. Require package 'season' ..
9 
##' @title
10 
##' @param temp monthly temperature in C
11 
##' @param threshold.sgdd threshold of temperature to compte sum of degree day default 5.56
12 
##' @return
13 
##' @author Kunstler
14 
fun.sgdd <-  function(temp,threshold.sgdd =  5.56){
15 
  require(season)
16 
  temp <-  unlist(temp)
17 
  ndays.month <- flagleap(data.frame(year=2013, month=1:12), F)$ndaysmonth
18 
  ndays.year <- sum(ndays.month)
19 
  x <- c(-ndays.month[12]/2,(ndays.month/2) + c(0, cumsum(ndays.month[1:11])),ndays.year+ndays.month[1]/2)
20 
   ## plot(x,c(temp[12],temp,temp[1]),xlim=c(0,365),type='b')
21 
  myfit <- loess(y~x, data=data.frame(y=c(temp[12],temp,temp[1]), x=x),span=0.4,degree=2)
22 
  mypred <- predict(myfit, 1:ndays.year)
23 
   ## lines( 1:ndays.year,mypred,col="red")
24 
  sgdd <- sum(mypred[mypred>=threshold.sgdd]-threshold.sgdd)
25 
  return(sgdd)
26 
}
27
28
29
30 
#################################################################################################
31 
### function to compute max soil water content based on Piedallu 2010 Geoderma
32 
fun.swhc <-  function(affroc,cailloux,text2,text1,prof2,prof1,codeprof,codecaillou,data.texture){
33 
## browser()
34
35 
### transform code of prof= soil depth in cm and code of percentage of rock in %
36 
prof2.t <- prof2
37 
prof2.t[is.na(prof2)] <- 0
38 
prof2.cm.temp <- codeprof[prof2+1]
39 
prof2.cm.temp[is.na(prof2)] <-0
40 
prof1.t <-  prof1
41 
prof1.t[is.na(prof1)] <- 0
42 
text1[is.na(text1)] <- 0
43 
prof1.cm <- codeprof[prof1.t+1]
44 
prof1.cm[is.na(prof1)] <- 0
45 
cailloux.t <-  cailloux
46 
cailloux.t[is.na(cailloux)] <-  0
47 
cailloux.perc <- codecaillou[cailloux.t+1]
48 
cailloux.perc[is.na(cailloux)] <- 0
49
50 
## compute depth of second horizon
51 
prof2.cm <- prof2.cm.temp - prof1.cm
52 
perc.top1 <- apply(cbind(10/prof1.cm,rep(1,length(prof1.cm))),MARGIN=1,FUN=min)
53 
perc.top2 <- apply(cbind(apply(cbind(rep(0,length(prof1.cm)),(10-prof1.cm)/prof2.cm.temp),MARGIN=1,FUN=max),
54 
                         rep(1,length(prof1.cm))),MARGIN=1,FUN=min)
55 
caill <- (1-codecaillou[affroc+1]/100)*(1-(sqrt(cailloux.perc/100))^3)
56 
## print(c(prof2.cm,perc.top2))
57 
top.hor <- (prof1.cm*perc.top1*data.texture[text1+1,"U_almajou_top"]+prof2.cm*perc.top2*data.texture[text2+1,"U_almajou_top"])
58
59 
sub.hor <- (prof1.cm*(1-perc.top1)*data.texture[text1+1,"U_almajou_sub"]+prof2.cm*(1-perc.top2)*data.texture[text2+1,"U_almajou_sub"])
60
61 
swhc <-  caill*(top.hor+sub.hor)
62 
return(swhc)
63 
}
64
65 
#################################################################
66 
## function to compute PET
67 
fun.PET <-  function(i,rad,temp){
68 
## conversion data
69 
  require(season)
70 
  ndays.month <- flagleap(data.frame(year=2013, month=1:12), F)$ndaysmonth
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121
Rs <- unlist(rad[i,])/100/30*1000 Ta <- unlist(temp[i,]) # compute TURC FORMULA vec.pet <- 0.4*((0.0239001*Rs)+50)*(Ta/(Ta+15))/30*ndays.month ## with Rs daily global (total) solar radiation (kJ/m2/day) ## and Ta monthly mean temp in C and ndays.month number of day in each month return(vec.pet) } #################################################################################### #################################################################################### #################################################################################### ### FUNCTION TO COMPUTE WATER BUDGET based on Bugmann & Cramer 1998 fun.WaterBudget <-function(i,prcp.m,PET.m,Ta.m,SWHC.v,n=2){ WATER <- rep(NA,12*n) WATER[1] <- unlist(SWHC.v)[i] SWHC <- unlist(SWHC.v)[i] AET <- rep(NA,12*n) D <- rep(NA,12*n) prcp <- rep(unlist(prcp.m[i,]),n) PET <- rep(unlist(PET.m[i,]),n) Ta <- rep(unlist(Ta.m[i,]),n) for (i in 2:(12*n)){ Pi <- min(0.3*prcp[i] , PET[i]) Ps <- prcp[i] - Pi S <- 120*WATER[i-1]/SWHC D[i] <- PET[i]-Pi AET[i] <- max(0,min(D[i],S)) WATER[i] <- max(0,min(WATER[i-1]+Ps-AET[i],SWHC)) } ## plot(1:(12*n),WATER,type="b") WATER <- WATER[13:24] Ta <- Ta[13:24] AET <- AET[13:24] D <- D[13:24] ## NEED to compute mean water availability of the growing season WB.s <- mean(WATER[Ta>5.56],na.rm=T) ## NEED to compute mean water availability of the year WB.y <- mean(WATER,na.rm=T) ## NEED to compute water stress index of the growing season WS.s <- sum(na.omit(AET[Ta>5.56 & D>0]))/sum(na.omit(D[Ta>5.56 & D>0])) ## NEED to compute water stress index of the year WS.y <- sum(na.omit(AET[D>0]))/sum(na.omit(D[D>0])) return(c("WB.s"=WB.s,"WB.y"=WB.y, "WS.s"=WS.s,"WS.y"=WS.y)) }