Generation of Precipitation in One Station
Occurrence
A precipitation occurrence generator model is created in order to reproduce a time series consistent with observed time series in:
it1 <- "T0083" #station[2]
it1
#> [1] "T0083"Let assume that precipitation occurrence also depend on the diurnal temperature variation which is taken into account as an exogenous variable. Therefore, a model based on a logistic regression with a generalized linear model is created and saved as model1 variable;
exogen <- Tx_mes[,it1]-Tn_mes[,it1]
months <- factor(prec_mes$month)
model1 <- PrecipitationOccurrenceModel(x=prec_mes[,it1],exogen=exogen,monthly.factor=months)The model is an S3 opject:
class(model1)
#> [1] "PrecipitationOccurrenceModel"constituting of a list of 5 elements:
class(unclass(model1))
#> [1] "list"
names(model1)
#> [1] "predictor" "glm" "p" "valmin" "id.name"The first element is a data frame containing the predictor values: the exogenous variable, the month of the year and the precipiation occurrence at the previous days:
str(model1$predictor)
#> 'data.frame': 10957 obs. of 3 variables:
#> $ exogen: num 7.59 3.59 5.59 5.59 14.59 ...
#> $ month : Factor w/ 12 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
#> $ x_l01 : logi NA FALSE TRUE FALSE TRUE FALSE ...Then, a genaralized linear model :
summary(model1$glm)
#>
#> Call:
#> glm(formula = df, family = "binomial")
#>
#> Deviance Residuals:
#> Min 1Q Median 3Q Max
#> -2.6037 -0.6806 -0.3892 0.5783 3.0965
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 1.440154 0.122650 11.742 < 2e-16 ***
#> exogen -0.330272 0.008389 -39.368 < 2e-16 ***
#> month2 0.164982 0.141512 1.166 0.2437
#> month3 0.904704 0.136761 6.615 3.71e-11 ***
#> month4 1.222853 0.134677 9.080 < 2e-16 ***
#> month5 1.736783 0.130607 13.298 < 2e-16 ***
#> month6 2.034278 0.130296 15.613 < 2e-16 ***
#> month7 1.883444 0.130523 14.430 < 2e-16 ***
#> month8 1.819358 0.130535 13.938 < 2e-16 ***
#> month9 1.250177 0.136727 9.144 < 2e-16 ***
#> month10 0.810081 0.138679 5.841 5.18e-09 ***
#> month11 0.272739 0.134311 2.031 0.0423 *
#> month12 -0.291167 0.139652 -2.085 0.0371 *
#> x_l01TRUE 0.791204 0.053903 14.678 < 2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> Null deviance: 12790.1 on 10955 degrees of freedom
#> Residual deviance: 9572.5 on 10942 degrees of freedom
#> AIC: 9600.5
#>
#> Number of Fisher Scoring iterations: 5The third element is the order of the autoregression, i.e. the number of previous days whose data are used as predictors in the genaralized linear model:
model1$p
#> [1] 1The minimum daily precipitation value above which the day ids considered as wet, e.g. 1 mm:
model1$valmin
#> [1] 0.5Then, a field name is considered:
model$id.nameOnce built the model, the time series of daily probabilities of precipitation occurrence can be predicted:
probs <- predict(model1$glm,type="response")where probs is the predicted probabability of precipitation occurrence releted to the time series of the observed values: each element of probs corresponds to the day of the corresponding row of the input dataset prec_mes[,it]. In case that the model is applied to a daily time series of precipitation occurrence and exogenous variale, i.e. predictors, they must be included as a new data frame equal to newdata argument:
row_test <- 2000:2007
newdata <- model1$predictor[row_test,]For instance, subsetting the observiations, a new dataset is created:
head(newdata)
#> exogen month x_l01
#> 2000 12.76 6 FALSE
#> 2001 13.76 6 FALSE
#> 2002 15.76 6 FALSE
#> 2003 16.76 6 TRUE
#> 2004 12.76 6 FALSE
#> 2005 14.76 6 FALSE
newdata$probs2 <- predict(model1,newdata=newdata)
head(newdata)
#> exogen month x_l01 probs2
#> 2000 12.76 6 FALSE 0.3230402
#> 2001 13.76 6 FALSE 0.2553828
#> 2002 15.76 6 FALSE 0.1505044
#> 2003 16.76 6 TRUE 0.2193055
#> 2004 12.76 6 FALSE 0.3230402
#> 2005 14.76 6 FALSE 0.1977561Once obtained the model, a random stochastic time series of daily precipitation occurrence is generated through generate method. First of all, the following predictors are needed:
- a time series of exogenous variables is needeed.
- a time series of month factor indicating the corresponding month of the year for each day of the new stochastically generated time series;
These two time series must be of the same length which corresponds to the length of the time series that is about to generate. For instance, if the time series is a replication of the observed one, it is:
exogen <- exogen
months <- factor(prec_mes$month)
set.seed(1235)
prec_gen1_occ <- generate(model1,exogen=exogen,monthly.factor=months,n=length(months))Amount
Once generated precipitation occurrence values, precipitation amount values are genereted through a stochastic generator model. The precipitation amount stochastic generator is built as follows:
model1_amount <- PrecipitationAmountModel(prec_mes,station=it1,origin=origin)where model_amount is an S3 object:
class(model1_amount)
#> [1] "PrecipitationAmountModel"containing the following items:
names(model1_amount)
#> [1] "T0083" "station" "sample" "x" "origin" "valmin"The first item is a generalized linear model, in which Gaussianized value of precipitation depends on month:
summary(model1_amount[[it1]])
#>
#> Call:
#> lm(formula = df)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.74134 -0.68036 -0.00999 0.66302 3.07671
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -0.01865 0.08008 -0.233 0.816
#> month2 0.08147 0.11382 0.716 0.474
#> month3 0.03716 0.10579 0.351 0.725
#> month4 0.02625 0.10285 0.255 0.799
#> month5 0.04111 0.09830 0.418 0.676
#> month6 0.04275 0.09792 0.437 0.662
#> month7 0.01087 0.09960 0.109 0.913
#> month8 0.04105 0.10012 0.410 0.682
#> month9 0.03853 0.10659 0.362 0.718
#> month10 0.04228 0.10849 0.390 0.697
#> month11 0.04162 0.10483 0.397 0.691
#> month12 0.05511 0.11699 0.471 0.638
#>
#> Residual standard error: 0.9873 on 2608 degrees of freedom
#> Multiple R-squared: 0.0003119, Adjusted R-squared: -0.003905
#> F-statistic: 0.07398 on 11 and 2608 DF, p-value: 1Other items are the following variables: names of the considered station (e.g. only “T0083”}), the origin date of the time series used to build the model and the minimum value of precipipition to consider as precipitation occurrence, in this case 1 mm.
model1_amount$station
#> [1] "T0083"
model1_amount$sample
#> [1] "monthly"
model1_amount$origin
#> [1] "1961-1-1"
model1_amount$valmin
#> [1] 1Finally, the x item is the data frame containing observed time series of precipitation:
str(model1_amount$x)
#> 'data.frame': 10957 obs. of 1 variable:
#> $ T0083: num 0 6.76 0 61.5 0 ...Finally, a time series, here called prec_gen1 , is generate by applying generate method to the model of precipitation account:
prec_gen1 <- generate(model1_amount,newdata=prec_gen1_occ)
str(prec_gen1)
#> 'data.frame': 10957 obs. of 1 variable:
#> $ T0083: num 0 14.87 31.25 2.66 0 ...Comparison
The stochastic coherence between observed and generated time series is verified through the Quantile-Quantile plot:
library(ggplot2)
library(lubridate)
library(reshape2)
df <- data.frame(obs=prec_mes[,it1],gen=prec_gen1[,it1])
df$date <- as.Date(origin)+days(1:nrow(df))-1
df$month <- factor(month(df$date))
df$season <- "none"
df$season[df$month %in% c(12,2,1)] <- "`1.DJF"
df$season[df$month %in% c(3,4,5)] <- "2.MAM"
df$season[df$month %in% c(6,7,8)] <- "3.JJA"
df$season[df$month %in% c(9,10,11)] <- "4.SON"
qqplot_ <- function(df) {
df <- as.list(df)
o <- qqplot(df[[1]],df[[2]],plot.it=FALSE)
names(o) <- names(df)[1:2]
o <- as.data.frame(o)
return(o)
}
qqdf <- split(df,f=df$season) %>% lapply(FUN=qqplot_) %>% melt(id=names(df)[1:2])
names(qqdf)[names(qqdf)=="L1"] <- "season"
g <- ggplot(data=qqdf)+geom_point(aes(x=obs,y=gen))+theme_bw()+geom_abline()+facet_grid(. ~ season)+xlab("observed")+ylab("generated")
show(g)
A comparison of the dry/wet spells, considered as periods of consecutive dry/rainy days, between observed and generated time series are plotted below respectively. The following lines of codes makes use of dw.spell function that calculates the length expressed in days for the dry or wet spells:
dw <- list()
dw$obs <- dw.spell(prec_mes[,it1],origin=origin)[[1]]
dw$obs$variable <- "obs"
dw$gen <- dw.spell(prec_gen1[,it1],origin=origin)[[1]]
dw$gen$variable <- "gen"
dw <- do.call(what=rbind,args=dw)
dw$season <- "none"
dw$season[dw$month %in% c(12,2,1)] <- "`1.DJF"
dw$season[dw$month %in% c(3,4,5)] <- "2.MAM"
dw$season[dw$month %in% c(6,7,8)] <- "3.JJA"
dw$season[dw$month %in% c(9,10,11)] <- "4.SON"
qqdw <- split(dw,f=paste(dw$spell_state,dw$season,sep="_")) %>% lapply(FUN=function(x){list(obs=x$spell_length[x$variable=="obs"],gen=x$spell_length[x$variable=="gen"])}) %>% lapply(FUN=qqplot_) %>% melt(id=c("obs","gen"))
qqdw_add <- str_split(qqdw$L1,pattern="_") %>% do.call(what=rbind) %>% as.data.frame()
names(qqdw_add) <- c("state","season")
qqdw <- cbind(qqdw,qqdw_add)Therefore, the quantile-quantile plots for the dry spells are:
qqdry <- qqdw[qqdw$state=="dry",]
## ggplot plot
ggdry <- ggplot(data=qqdry)+geom_point(aes(x=obs,y=gen))+theme_bw()+geom_abline()+facet_grid(. ~ season)+xlab("observed")+ylab("generated")
show(ggdry)
Analogously for the wet spells it is:
qqwet <- qqdw[qqdw$state=="wet",]
## ggplot plot
ggwet <- ggplot(data=qqwet)+geom_point(aes(x=obs,y=gen))+theme_bw()+geom_abline()+facet_grid(. ~ season)+xlab("observed")+ylab("generated")
show(ggwet)
It is a comparison of sample probability distributtions of daily precipitation, dry and wet spells between the one obtained by the observed time series and the generated time series.
and the quantile-quantile plots for the wet spells are:
