Skip to content

GitLab

Commit cd73ff5b authored by unknown's avatar unknown
Browse files

v1.0.5.16 function DataAltiExtrapolation_Valery() cleaned

parent 27f26347
Hide whitespace changes
Inline Side-by-side
DESCRIPTION
View file @ cd73ff5b
Package: airGR
Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
Version: 1.0.5.15
Version: 1.0.5.16
Date: 2017-01-23
Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl")),
......
R/DataAltiExtrapolation_Valery.R
View file @ cd73ff5b
DataAltiExtrapolation_Valery <- function(DatesR, Precip, PrecipScale = TRUE, TempMean, TempMin = NULL, TempMax = NULL, ZInputs, HypsoData, NLayers, verbose = TRUE){
DataAltiExtrapolation_Valery <- function(DatesR,
Precip, PrecipScale = TRUE,
TempMean, TempMin = NULL, TempMax = NULL,
ZInputs, HypsoData, NLayers,
verbose = TRUE) {
##Altitudinal_gradient_functions_______________________________________________________________
##unique_gradient_for_precipitation
GradP_Valery2010 <- function(){
return(0.00041); ### value from Val? PhD thesis page 126
}
##daily_gradients_for_mean_min_and_max_air_temperature
GradT_Valery2010 <- function(){
RESULT <- matrix(c(
##Altitudinal_gradient_functions_______________________________________________________________
##unique_gradient_for_precipitation
GradP_Valery2010 <- function() {
return(0.00041) ### value from Valerie PhD thesis page 126
}
##daily_gradients_for_mean_min_and_max_air_temperature
GradT_Valery2010 <- function() {
RESULT <- matrix(c(
01, 01, 0.434, 0.366, 0.498,
02, 01, 0.434, 0.366, 0.500,
03, 01, 0.435, 0.367, 0.501,
......@@ -374,130 +377,170 @@ DataAltiExtrapolation_Valery <- function(DatesR, Precip, PrecipScale = TRUE, Tem
28, 12, 0.431, 0.366, 0.495,
29, 12, 0.431, 0.366, 0.495,
30, 12, 0.432, 0.366, 0.496,
31, 12, 0.433, 0.366, 0.497),ncol=5,byrow=TRUE);
dimnames(RESULT) <- list(1:366,c("day","month","grad_Tmean","grad_Tmin","grad_Tmax"));
return(RESULT);
}
31, 12, 0.433, 0.366, 0.497), ncol = 5, byrow = TRUE)
dimnames(RESULT) <- list(1:366, c("day", "month", "grad_Tmean", "grad_Tmin", "grad_Tmax"))
return(RESULT)
}
##Format_______________________________________________________________________________________
HypsoData <- as.double(HypsoData);
ZInputs <- as.double(ZInputs);
HypsoData <- as.double(HypsoData)
ZInputs <- as.double(ZInputs)
##ElevationLayers_Creation_____________________________________________________________________
ZLayers <- as.double(rep(NA,NLayers));
if(!identical(HypsoData,as.double(rep(NA,101)))){
nmoy <- 100 %/% NLayers;
nreste <- 100 %% NLayers;
ncont <- 0;
for(iLayer in 1:NLayers){
if(nreste > 0){ nn <- nmoy+1; nreste <- nreste-1; } else { nn <- nmoy; }
if(nn==1){ ZLayers[iLayer] <- HypsoData[ncont+1]; }
if(nn==2){ ZLayers[iLayer] <- 0.5 * (HypsoData[ncont+1] + HypsoData[ncont+2]); }
if(nn>2 ){ ZLayers[iLayer] <- HypsoData[ncont+nn/2]; }
ncont <- ncont+nn;
ZLayers <- as.double(rep(NA, NLayers))
if (!identical(HypsoData, as.double(rep(NA, 101)))) {
nmoy <- 100 %/% NLayers
nreste <- 100 %% NLayers
ncont <- 0
for (iLayer in 1:NLayers) {
if (nreste > 0) {
nn <- nmoy + 1
nreste <- nreste - 1
} else {
nn <- nmoy
}
if (nn == 1) {
ZLayers[iLayer] <- HypsoData[ncont + 1]
}
if (nn == 2) {
ZLayers[iLayer] <- 0.5 * (HypsoData[ncont + 1] + HypsoData[ncont + 2])
}
if (nn > 2) {
ZLayers[iLayer] <- HypsoData[ncont + nn / 2]
}
ncont <- ncont + nn
}
}
}
##Precipitation_extrapolation__________________________________________________________________
##Initialisation
if (identical(ZInputs, HypsoData[51]) & NLayers == 1) {
LayerPrecip <- list(as.double(Precip))
} else {
##Elevation_gradients_for_daily_mean_precipitation
GradP <- GradP_Valery2010() ### single value
TabGradP <- rep(GradP, length(Precip))
##Extrapolation
##Thresold_of_inputs_median_elevation
Zthreshold <- 4000
LayerPrecip_df <- sapply(1:NLayers, function(iLayer) {
##If_layer_elevation_smaller_than_Zthreshold
if (ZLayers[iLayer] <= Zthreshold) {
prcp <- as.double(Precip * exp(TabGradP * (ZLayers[iLayer] - ZInputs)))
##If_layer_elevation_greater_than_Zthreshold
} else {
##If_inputs_median_elevation_smaller_than_Zthreshold
if (ZInputs <= Zthreshold) {
prcp <- as.double(Precip * exp(TabGradP*(Zthreshold - ZInputs)))
##If_inputs_median_elevation_greater_then_Zthreshold
if (identical(ZInputs, HypsoData[51]) & NLayers == 1) {
LayerPrecip <- list(as.double(Precip))
} else {
##Elevation_gradients_for_daily_mean_precipitation
GradP <- GradP_Valery2010() ### single value
TabGradP <- rep(GradP, length(Precip))
##Extrapolation
##Thresold_of_inputs_median_elevation
Zthreshold <- 4000
LayerPrecip_df <- sapply(1:NLayers, function(iLayer) {
##If_layer_elevation_smaller_than_Zthreshold
if (ZLayers[iLayer] <= Zthreshold) {
prcp <- as.double(Precip * exp(TabGradP * (ZLayers[iLayer] - ZInputs)))
##If_layer_elevation_greater_than_Zthreshold
} else {
prcp <- as.double(Precip)
##If_inputs_median_elevation_smaller_than_Zthreshold
if (ZInputs <= Zthreshold) {
prcp <- as.double(Precip * exp(TabGradP * (Zthreshold - ZInputs)))
##If_inputs_median_elevation_greater_then_Zthreshold
} else {
prcp <- as.double(Precip)
}
}
return(prcp)
})
if (PrecipScale) {
LayerPrecip_df <- LayerPrecip_df / rowMeans(LayerPrecip_df) * Precip
LayerPrecip_df[is.nan(LayerPrecip_df)] <- 0
}
return(prcp)
}
)
if (PrecipScale) {
LayerPrecip_df <- LayerPrecip_df / rowMeans(LayerPrecip_df) * Precip
LayerPrecip_df[is.nan(LayerPrecip_df)] <- 0
LayerPrecip <-
lapply(seq_len(ncol(LayerPrecip_df)), function(x) LayerPrecip_df[, x])
}
LayerPrecip <- lapply(seq_len(ncol(LayerPrecip_df)), function(x) LayerPrecip_df[, x])
}
##Temperature_extrapolation____________________________________________________________________
##Initialisation
LayerTempMean <- list(); LayerTempMin <- list(); LayerTempMax <- list();
if(identical(ZInputs,HypsoData[51]) & NLayers==1) {
LayerTempMean[[1]] <- as.double(TempMean);
if(!is.null(TempMin) & !is.null(TempMax)){ LayerTempMin[[1]] <- as.double(TempMin); LayerTempMax[[1]] <- as.double(TempMax); }
} else {
##Elevation_gradients_for_daily_mean_min_and_max_temperature
GradT <- as.data.frame(GradT_Valery2010())
iday <- match(format(DatesR,format="%d%m"), sprintf("%02i%02i", GradT[, "day"], GradT[, "month"]))
TabGradT <- GradT[iday, c("grad_Tmean", "grad_Tmin", "grad_Tmax")]
##Extrapolation
##On_each_elevation_layer...
for(iLayer in 1:NLayers){
LayerTempMean[[iLayer]] <- as.double(TempMean + (ZInputs-ZLayers[iLayer])*abs(TabGradT[, "grad_Tmean"])/100);
if(!is.null(TempMin) & !is.null(TempMax)){
LayerTempMin[[iLayer]] <- as.double(TempMin + (ZInputs-ZLayers[iLayer])*abs(TabGradT[, "grad_Tmin"])/100);
LayerTempMax[[iLayer]] <- as.double(TempMax + (ZInputs-ZLayers[iLayer])*abs(TabGradT[, "grad_Tmax"])/100);
LayerTempMean <-
list()
LayerTempMin <- list()
LayerTempMax <- list()
if (identical(ZInputs, HypsoData[51]) & NLayers == 1) {
LayerTempMean[[1]] <- as.double(TempMean)
if (!is.null(TempMin) & !is.null(TempMax)) {
LayerTempMin[[1]] <- as.double(TempMin)
LayerTempMax[[1]] <- as.double(TempMax)
}
} else {
##Elevation_gradients_for_daily_mean_min_and_max_temperature
GradT <- as.data.frame(GradT_Valery2010())
iday <- match(format(DatesR, format = "%d%m"),
sprintf("%02i%02i", GradT[, "day"], GradT[, "month"]))
TabGradT <-
GradT[iday, c("grad_Tmean", "grad_Tmin", "grad_Tmax")]
##Extrapolation
##On_each_elevation_layer...
for (iLayer in 1:NLayers) {
LayerTempMean[[iLayer]] <- as.double(TempMean + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmean"]) / 100)
if (!is.null(TempMin) & !is.null(TempMax)) {
LayerTempMin[[iLayer]] <- as.double(TempMin + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmin"]) / 100)
LayerTempMax[[iLayer]] <- as.double(TempMax + (ZInputs - ZLayers[iLayer]) * abs(TabGradT[, "grad_Tmax"]) / 100)
}
}
}
}
##Solid_Fraction_for_each_elevation_layer______________________________________________________
LayerFracSolidPrecip <- list();
##Thresold_of_inputs_median_elevation
Zthreshold <- 1500;
##On_each_elevation_layer...
for(iLayer in 1:NLayers){
Option <- "USACE";
if(!is.na(ZInputs)){ if(ZInputs < Zthreshold & !is.null(TempMin) & !is.null(TempMax)){ Option <- "Hydrotel"; } }
##Turcotte_formula_from_Hydrotel
if(Option=="Hydrotel"){
TempMin <- LayerTempMin[[iLayer]];
TempMax <- LayerTempMax[[iLayer]];
SolidFraction <- 1 - TempMax/(TempMax - TempMin);
SolidFraction[TempMin >= 0] <- 0;
SolidFraction[TempMax <= 0] <- 1;
}
##USACE_formula
if(Option=="USACE"){
USACE_Tmin <- -1.0;
USACE_Tmax <- 3.0;
TempMean <- LayerTempMean[[iLayer]];
SolidFraction <- 1- (TempMean - USACE_Tmin)/(USACE_Tmax - USACE_Tmin);
SolidFraction[TempMean > USACE_Tmax] <- 0;
SolidFraction[TempMean < USACE_Tmin] <- 1;
LayerFracSolidPrecip <- list()
##Thresold_of_inputs_median_elevation
Zthreshold <- 1500
##On_each_elevation_layer...
for (iLayer in 1:NLayers) {
Option <- "USACE"
if (!is.na(ZInputs)) {
if (ZInputs < Zthreshold & !is.null(TempMin) & !is.null(TempMax)) {
Option <- "Hydrotel"
}
}
##Turcotte_formula_from_Hydrotel
if (Option == "Hydrotel") {
TempMin <- LayerTempMin[[iLayer]]
TempMax <- LayerTempMax[[iLayer]]
SolidFraction <- 1 - TempMax / (TempMax - TempMin)
SolidFraction[TempMin >= 0] <- 0
SolidFraction[TempMax <= 0] <- 1
}
##USACE_formula
if (Option == "USACE") {
USACE_Tmin <- -1.0
USACE_Tmax <- 3.0
TempMean <- LayerTempMean[[iLayer]]
SolidFraction <- 1 - (TempMean - USACE_Tmin) / (USACE_Tmax - USACE_Tmin)
SolidFraction[TempMean > USACE_Tmax] <- 0
SolidFraction[TempMean < USACE_Tmin] <- 1
}
LayerFracSolidPrecip[[iLayer]] <- as.double(SolidFraction)
}
LayerFracSolidPrecip[[iLayer]] <- as.double(SolidFraction);
}
##END__________________________________________________________________________________________
return(list(LayerPrecip=LayerPrecip,LayerTempMean=LayerTempMean,LayerTempMin=LayerTempMin,LayerTempMax=LayerTempMax,
LayerFracSolidPrecip=LayerFracSolidPrecip,ZLayers=ZLayers));
return(list(LayerPrecip = LayerPrecip,
LayerTempMean = LayerTempMean,
LayerTempMin = LayerTempMin,
LayerTempMax = LayerTempMax,
LayerFracSolidPrecip = LayerFracSolidPrecip,
ZLayers = ZLayers))
}
Markdown is supported
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