Download: Windows application, OS X Console application, C++ code
# ========= NEW VERSION: ReGen 2.0 ================================
# Parameters derived from stations in Israel and western Jordan
# Time units of parameters can be set to days or months
# Simulation of hourly values based on analysis of data from four field stations
# More flexible Gauss function (variable 'shape' parameter)
#
# Martin Köchy, 2008-09-02
# Gauss peak curve function for calculating daily precipitation
# probability and mean rain volume
# NUMBERING OF MONTHS AND DAYS STARTS WITH August 1 = 0
# The Gauss curve for time measured in DAYS
gauss.d <- function (day, amplitude, location, width, shape=2)
{ day = (day + 182 - location) %% 365 + location - 182
G = amplitude*exp(-abs(day-location)^shape/(2*width^shape))
return(G)
}
# The Gauss curve for time measured in MONTHS
gauss.m <- function (month, amplitude, location, width, shape=2)
{ month = (month + 6 - location) %% 12 + location - 6
G = amplitude*exp(-abs(month-location)^shape/(2*width^shape))
return(G)
}
# Generic Gauss curve, picks the appropriate function
gauss <- function (time.units, amplitude, location, width, shape=2)
{
if (max(time.units)>300) return(gauss.d(time.units, amplitude, location, width, shape) else
return(gauss.m(time.units, amplitude, location, width, shape)
}
# Calculation of PARAMETERS based on mean annual precipitation (MAP)
# and change of daily mean rain (relCh)
# Parameters derived from stations in ISRAEL and eastern JORDAN.
# THIS IS THE VERSION FOR DAYS
rmp.ILWJO.d<-function(MAP, relCh=0) {
if(min(MAP) < 80 | max(MAP) > 960 | min(relCh) < -0.3 | max(relCh) > 0.3) {
warning("This function has been validated only for MAP between 90 and 800 mm
and a factor between -0.3 and +0.3 for the rain volume") }
MAP = (MAP+6)/1.03 # fine-tuning of output
A = 0.1082772 + 0.0004189*MAP # R2 = 0.78, n=64
A2 = A*1*(1-relCh + (0.3615545 - 0.0579928*log(MAP))*relCh + relCh^2 -0.4*relCh^3 )
L = 182.35245 - 26.135966*A # R2 = 0.44, n=64
W = 31.24488 + 4.7929879*log(MAP) # R2 = 0.43, n=63
S = -0.360935 -2.55397*A + 0.05713082*W + 0.000465807*MAP # R2 = 0.82, n=64
a = -5.3822417 + 0.01660176* MAP + 4.17701974*log(MAP) -40.243339*A -0.0666435*W # R2 = 0.94, n=64
a2 = a*(1+relCh)
l = 191.510119 - 48.62739*A # R2 = 0.43, n=62
w = -3.501397 + 0.0231263*MAP - 2.1062518*a + 0.67337071*l # R2 = 0.34, n=63
s = 72.8194938 -0.3949018*W + 6.20975069*S - 1.4508933*w + 0.00861436*w^2 # R2 = 0.90, n=54
return(data.frame(P.amplitude=A2, P.location=L, P.width=W, P.shape=S, V.amplitude=a2, V.location=l, V.width=w, V.shape=s)) }



# =========== PUBLISHED VERSION: ReGen 1.0 Parameters ==============================================
# Calculation of parameters based on mean annual precipitation (MAP)
# and change of daily mean rain (relCh)
# Parameters derived from stations in ISRAEL only. Corresponds to the published version.
# THIS IS THE VERSION FOR DAYS
rmp.iIL.d<-function (MAP, relCh=0)
{ if(min(MAP) < 80 | max(MAP) > 960 | min(relCh) < -0.3 | max(relCh) > 0.3) {
warning("This function has been validated only for MAP between 90 and 800 mm
and a factor between -0.3 and +0.3 for the rain volume") }
A = 0.130034 + 0.000408 * MAP # rain occurance
A = A*(1-relCh +1.33*relCh^2-(0.61+1.57*A)*relCh^3)
L = 177
W = 52 + 0.007 * MAP
a = -14 + 4 * log(MAP) # rain volume
a = a * (1+relCh)
l = 170
w = 102.4109 # differs from published version to reflect that 'shape' in the Gauss function
# is now a variable but produces identical results
return(data.frame(P.amplitude=A, P.location=L, P.width=W, P.shape=2, V.amplitude=a, V.location=l, V.width=w, V.shape=4))
}
# Calculation of parameters based on mean annual precipitation (MAP)
# and change of daily mean rain (relCh)
# Parameters derived from stations in ISRAEL only.
# THIS IS THE ORIGINAL VERSION FOR MONTHS.
rmp.iIL.m<-function (MAP, relCh=0) {
if(min(MAP) < 80 | max(MAP) > 960 | min(relCh) < -0.3 | max(relCh) > 0.3) {
warning("This function has been validated only for MAP between 90 and 800 mm
and a factor between -0.3 and +0.3 for the rain volume") }
A = 0.130034 + 0.000408 * MAP # rain occurance
A = A*(1-relCh +1.33*relCh^2-(0.61+1.57*A)*relCh^3)
L = 5.735
W = 1.704746 + 0.000242 * MAP
a = -13.7178 + 4.0508301 * log(MAP) # rain volume
a = a * (1+relCh)
l = 5.586
w = 11.336 # differs from published version to reflect that 'shape' in the Gauss function is now
# a variable but produces identical results
return(data.frame(P.amplitude=A, P.location=L, P.width=W, V.amplitude=a, V.location=l, V.width=w))}
# ======================End ReGen 1.0 Parameters ====================================================
# Wrapper function to produce the parameters for the stochastic time series.
# The user can choose time units (default is 'days')
# and the region (interior Israel or Jordan River valley) used for the parameterization;
# 'interior Israel' corresponds to the publsihed version.
rmp<-function(MAP, relCh=0, region='jordanvalley', time.units="d") {
switch(region,
'jordanvalley' = {
switch(time.units,
'd' = rmp.ILWJO.d(MAP, relCh),
warning("No function for this time units for this region"))},
'interiorIsrael' = {
switch(time.units,
'd' = rmp.iIL.d(MAP, relCh),
'm' = rmp.iIL.m(MAP, relCh),
warning("No function for this time unit for this region"))},
warning("No function for this region."))
}
############# Production of stochastic time series ############################
# Parameters can be entered directly or as 'scenario' calculated by function #
# 'rmp'. #
# ReGen returns a matrix with 365 colummns and one row for each year. #
ReGen<-function (years, P.amplitude=0, P.location=0, P.width=0, P.shape=2,
amplitude=0, location=0, width=0, shape=4, scenario=0)
{ DM = matrix(rep(1:365, years),nrow=years, byrow=T)
if(is.data.frame(scenario))
{ P.amplitude=as.numeric(scenario$P.amplitude)
P.location=as.numeric(scenario$P.location)
P.width=as.numeric(scenario$P.width)
P.shape=as.numeric(scenario$P.shape)
V.amplitude=as.numeric(scenario$V.amplitude)
V.location=as.numeric(scenario$V.location)
V.width=as.numeric(scenario$V.width)
V.shape=as.numeric(scenario$V.shape)
}
# probability of a rainy day
G<-gauss.d(DM, P.amplitude, P.location, P.width, P.shape)
RT<-ifelse(runif(DM)0.05, 1, 0)
# daily mean rain volume:
RV<- gauss.d(DM, V.amplitude, V.location, V.width, V.shape)
# actual rain volume:
R<- RT * ifelse(RV>0.01,rexp(DM, 1/pmax(RV,0.01)),0) # pmax: prevent division by zero
return(R)
}
# ========= Break down the daily rain amount into hourly amounts ===============================
ReGen.Tag <- function(daily.amount) {
# derived from measurements at 4 Israeli field sites
a = 0.07*daily.amount # a = p/(1-p), p = rainy hours/d
h = round(a/(1+a)*24)
d1 = seq(0, floor((24-h)/2))
d2 = 24-d1
d = c(d1, seq(daily.amount/h, h), seq(0, d2))
}