Commit fc6577a2 authored by Daniel Falster's avatar Daniel Falster
Browse files

Reorgansie files

No related merge requests found
Showing with 2 additions and 144 deletions
+2 -144
output
output_cluster
data/*
ignore
data/raw
figs
*.xls
*.xlsx
......
##################################### FUNCTION TO COMPUTE CLIMATIC DATA FOR FRANCE
## function to compute sum of degree days above 5.56
##' .. Compute teh sum of degree days above 5.56C from monthly data with a spline ..
##'
##' .. A spline is used to smooth the variation of T over each day from mothly data. Require package 'season' ..
##' @title
##' @param temp monthly temperature in C
##' @param threshold.sgdd threshold of temperature to compte sum of degree day default 5.56
##' @return
##' @author Kunstler
fun.sgdd <- function(temp, threshold.sgdd = 5.56) {
require(season)
temp <- unlist(temp)
ndays.month <- flagleap(data.frame(year = 2013, month = 1:12), F)$ndaysmonth
ndays.year <- sum(ndays.month)
x <- c(-ndays.month[12]/2, (ndays.month/2) + c(0, cumsum(ndays.month[1:11])),
ndays.year + ndays.month[1]/2)
## plot(x,c(temp[12],temp,temp[1]),xlim=c(0,365),type='b')
myfit <- loess(y ~ x, data = data.frame(y = c(temp[12], temp, temp[1]), x = x),
span = 0.4, degree = 2)
mypred <- predict(myfit, 1:ndays.year)
## lines( 1:ndays.year,mypred,col='red')
sgdd <- sum(mypred[mypred >= threshold.sgdd] - threshold.sgdd)
return(sgdd)
}
################################################################################################# function to compute max soil water content based on Piedallu 2010 Geoderma
fun.swhc <- function(affroc, cailloux, text2, text1, prof2, prof1, codeprof, codecaillou,
data.texture) {
## browser()
### transform code of prof= soil depth in cm and code of percentage of rock in %
prof2.t <- prof2
prof2.t[is.na(prof2)] <- 0
prof2.cm.temp <- codeprof[prof2 + 1]
prof2.cm.temp[is.na(prof2)] <- 0
prof1.t <- prof1
prof1.t[is.na(prof1)] <- 0
text1[is.na(text1)] <- 0
prof1.cm <- codeprof[prof1.t + 1]
prof1.cm[is.na(prof1)] <- 0
cailloux.t <- cailloux
cailloux.t[is.na(cailloux)] <- 0
cailloux.perc <- codecaillou[cailloux.t + 1]
cailloux.perc[is.na(cailloux)] <- 0
## compute depth of second horizon
prof2.cm <- prof2.cm.temp - prof1.cm
perc.top1 <- apply(cbind(10/prof1.cm, rep(1, length(prof1.cm))), MARGIN = 1,
FUN = min)
perc.top2 <- apply(cbind(apply(cbind(rep(0, length(prof1.cm)), (10 - prof1.cm)/prof2.cm.temp),
MARGIN = 1, FUN = max), rep(1, length(prof1.cm))), MARGIN = 1, FUN = min)
caill <- (1 - codecaillou[affroc + 1]/100) * (1 - (sqrt(cailloux.perc/100))^3)
## print(c(prof2.cm,perc.top2))
top.hor <- (prof1.cm * perc.top1 * data.texture[text1 + 1, "U_almajou_top"] +
prof2.cm * perc.top2 * data.texture[text2 + 1, "U_almajou_top"])
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"])
swhc <- caill * (top.hor + sub.hor)
return(swhc)
}
################################################################# function to compute PET
fun.PET <- function(i, rad, temp) {
## conversion data
require(season)
ndays.month <- flagleap(data.frame(year = 2013, month = 1:12), F)$ndaysmonth
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))
}
File moved
File moved
File moved
File moved
File moved
File moved
File moved
File moved
File moved
File moved
File moved
File moved
File moved
# use formatR package to tidy code
library(formatR)
for (f in dir(".", pattern = "*.R")) {
cat("Cleaning ", f, "\n")
tidy.source(f, file = f)
}
for (f in c(dir(".", pattern = "*.R"), dir("R", pattern = "*.R", full.names = TRUE))) {
cat("Cleaning ", f, "\n")
tidy.source(f, file = f)
}
File moved
File moved
File moved
File moved
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment