Commit 9812d7f2 authored by Manuela's avatar Manuela

function prsim.weather

validation script
adapt documentation and example files
parent 752181b9
......@@ -7,3 +7,6 @@ src/*.o
src/*.so
*.*~
# function different margins
R/fun_stoch_sim_wave_diff_margins.R
\ No newline at end of file
......@@ -6,13 +6,13 @@ importFrom('splus2R', 'ifelse1')
importFrom("methods", "as")
importFrom("stats", "deltat", "time")
importFrom('wavScalogram', 'cwt_wst')
importFrom('splus2R', 'ifelse1')
# importFrom("ismev","gev.fit")
# importFrom("evd","rgev","pgev")
# importFrom("GB2","mlfit.gb2","rgb2","pgb2")
export("prsim")
export("prsim.wave")
export("prsim.weather")
useDynLib(PRSim, .registration = TRUE)
This diff is collapsed.
# short demo for prsim.weather
### load data for four stations
data(weather_sim_multi_sites)
sim <- weather_sim_multi_sites
### define plotting colors
col_sim <- adjustcolor("#fd8d3c",alpha=0.8)
col_sim_tran <- adjustcolor("#fd8d3c",alpha=0.2)
col_obs <- adjustcolor( "black", alpha.f = 0.2)
### greys
col_vect_obs <- c('#cccccc','#969696','#636363','#252525')
### oranges
col_vect_sim <- c('#fdbe85','#fd8d3c','#e6550d','#a63603')
### plot time series for multiple sites
### Temperature (first list entry)
par(mfrow=c(2,1),mar=c(3,3,2,1))
### determine ylim
ylim_max <- max(sim[[1]][[1]]$Temp)*1.5
### observed
plot(sim[[1]][[1]]$Temp[1:1000],ylab=expression(bold(paste("Temperature [degrees]"))),xlab="Time [d]",type="l",col=col_vect_obs[1],ylim=c(0,ylim_max),main='Observations')
for(l in 2:4){
lines(sim[[l]][[1]]$Temp[1:1000],col=col_vect_obs[l])
}
# legend('topleft',legend=c('Station 1','Station 2','Station 3','Station 4'),lty=1,col=col_vect_obs[1:4])
### simulated (one run)
plot(sim[[1]][[1]]$r1[1:1000],ylab=expression(bold(paste("Temperature [degrees]"))),xlab="Time [d]",type="l",col=col_vect_sim[1],ylim=c(0,ylim_max),main='Stochastic simulations')
for(l in 2:4){
lines(sim[[l]][[1]]$r1[1:1000],col=col_vect_sim[l])
}
### precipitation (second list entry)
ylim_max <- max(sim[[1]][[2]]$Prec)*1
### observed
plot(sim[[1]][[2]]$Prec[1:1000],ylab=expression(bold(paste("Precipitation [mm/d]"))),xlab="Time [d]",type="l",col=col_vect_obs[1],ylim=c(0,ylim_max),main='Observations')
for(l in 2:4){
lines(sim[[l]][[2]]$Prec[1:1000],col=col_vect_obs[l])
}
# legend('topleft',legend=c('Station 1','Station 2','Station 3','Station 4'),lty=1,col=col_vect_obs[1:4])
### simulated (one run)
plot(sim[[1]][[2]]$r1[1:1000],ylab=expression(bold(paste("Precipitation [mm/d]"))),xlab="Time [d]",type="l",col=col_vect_sim[1],ylim=c(0,ylim_max),main='Stochastic simulations')
for(l in 2:4){
lines(sim[[l]][[2]]$r1[1:1000],col=col_vect_sim[l])
}
###===============================###===============================###
### Application of PRSim.weather for temperature and precipitation
###===============================###===============================###
### (1) apply function with default distributions
### temperature: SEP, precipitation: E-GPD
### does not allow for extrapolation to yet unobserved values
###===============================###===============================###
out <- prsim.weather(data_p=data_p, data_t=data_t, number_sim=5, p_margin='egpd',t_margin='sep')
### save example simulation data
# setwd("~/PRSim-devel/data")
# save(file='weather_sim_multi_sites.rda',out)
### (2) example with alternative distributions
### rCDF and CDF_fit need to be defined
### temperature: normal, precipitation: GEV
###===============================###===============================###
### define normal distribution
library(fitdistrplus)
rNORM <- function(n, theta) rnorm(n, theta[1], theta[2])
pNORM <- function(x, theta) pnorm(x, theta[1], theta[2])
NORM_fit <- function( xdat, ...) fitdistr(xdat, 'normal', show=FALSE, ...)$estimate
### define GEV distribution
require("evd")
require("ismev")
rGEV <- function(n, theta) rgev(n, theta[1], theta[2], theta[3])
pGEV <- function(x, theta) pgev(x, theta[1], theta[2], theta[3])
GEV_fit <- function( xdat, ...) gev.fit(xdat, show=FALSE, ...)$mle
### apply function using alternative distributions
out <- prsim.weather(data_p=data_p, data_t=data_t, number_sim=1,p_margin='GEV',t_margin='NORM')
out <- prsim.weather(data_p=data_p, data_t=data_t, number_sim=1,p_margin='NORM',t_margin='NORM')
......@@ -18,7 +18,9 @@ Contains two functions for the stochastic simulation of continuous discharge tim
prsim: Simulation in the frequency domain is based on the randomization of the phases of the Fourier transform. We here combine phase randomization simulation with the flexible, four-parameter kappa distribution, which allows for the extrapolation to yet unobserved low and high flows. Alternative distributions or the empirical distribution can be used instead. The
simulation approach consists of eight steps: (1) fitting of theoretical Kappa distribution, (2) normalization and deseasonalization, (3) Fourier transformation, (4) Fourier phases computation, (5) random phase generation, (6) inverse Fourier transformation, (7) back transformation, and (8) simulation.
prsim.wave: Simulation in the frequency domain based on the randomization of the phases of the continuous wavelet transform. We combine phase randomization with the flexible, four-parameter kappa distribution. Alternative theoretical distributions or the empirical distribution can be used instead. The simulation procedure consists of five steps: (1) Derivation of random phases from a white noise time series, (2) Fitting of kappa distribution, (3) Wavelet transform, (4) Inverse wavelet transform, and (5) Transformation to the kappa distribution (or the distribution of choice).
prsim.wave: Simulation for multiple sites in the frequency domain based on the randomization of the phases of the continuous wavelet transform. We combine phase randomization with the flexible, four-parameter kappa distribution. Alternative theoretical distributions or the empirical distribution can be used instead. The simulation procedure consists of five steps: (1) Derivation of random phases from a white noise time series, (2) Fitting of kappa distribution, (3) Wavelet transform, (4) Inverse wavelet transform, and (5) Transformation to the kappa distribution (or the distribution of choice).
prsim.weather: Simulation of two variables (temperature and precipitation) for multiple sites in the frequency domain based on the randomization of the phases of the continuous wavelet transform. We combine phase randomization with the flexible, skewed exponential power (sep) and extended generalized pareto distributions (egpd). Alternative theoretical distributions can be used instead. The simulation procedure consists of five steps: (1) Derivation of random phases from a randomly sampled time series, (2) Fitting of temperature and precipitation disstributions, (3) Wavelet transform, (4) Inverse wavelet transform, and (5) Transformation to the desired distributions.
}
\author{
\packageAuthor{PRSim}
......@@ -29,6 +31,8 @@ Maintainer: \packageMaintainer{PRSim}
Brunner, M. I., A. Bárdossy, and R. Furrer (2019). Technical note: Stochastic simulation of streamflow time series using phase randomization. Hydrology and Earth System Sciences, 23, 3175-3187, https://doi.org/10.5194/hess-23-3175-2019.
Brunner, M. I., and E. Gilleland (2020). Stochastic simulation of streamflow and spatial extremes: a continuous, wavelet-based approach, Hydrology and Earth System Sciences, https://doi.org/10.5194/hess-24-3967-2020.
Brunner, M. I., and E. Gilleland (2021). Spatial compound hot-dry events in the United States: assessment using a multi-site multi-variable weather generator, in preparation.
}
\keyword{ package }
\examples{
......@@ -37,5 +41,7 @@ demo("PRSim")
demo("PRSim-validate")
demo("PRSim_wave")
demo("PRSim_wave-validate")
demo("PRSim_weather")
demo("PRSim_weather-validate")
}
}
\name{pRsim.weather}
\alias{PRsim.weather}
\alias{prsim.weather}
\alias{prsim_weather}
\title{Weather simulation (temperature and precipitation) for multiple stations}
\description{
Applies the wavelet-based weather simulation algorithm to multiple sites (single site possible as well)
}
\usage{
prsim.weather(data_p, data_t, station_id_p="Precip",station_id_t="Temp", number_sim=1, win_h_length=15, n_wave=100,verbose=TRUE,t_margin='sep',p_margin='egpd',...)
}
\arguments{
\item{data_p}{list of precipitation data frames. One list entry, i.e. data frame, corresponds to one station/grid cell. Each data frame contains the time indications and precipitation of one station. See \sQuote{Details}.}
\item{data_t}{list of temperature data frames. One list entry, i.e. data frame, corresponds to one station/grid cell. Each data frame contains the time indications and temperature of one station. See \sQuote{Details}.}
\item{station_id_p}{identifies the precipitation variable name in case several time series are present in \code{data_p}. See \sQuote{Details}.}
\item{station_id_t}{identifies the temperature variable name in case several time series are present in \code{data_t}. See \sQuote{Details}.}
\item{number_sim}{number of simulations to be carried out.}
\item{win_h_length}{(half-)length of moving window size.}
\item{t_margin}{marginal distribution to be used for the backtransformation of temperature. Can be either \code{"sep"} or any type of CDF (see \sQuote{Details}). \code{"sep"} uses the four-parameter skewed-exponential power for backtransformation. CDF allows for specifying any distribution \sQuote{Examples}.}
\item{p_margin}{marginal distribution to be used for the backtransformation of precipitation. Can be either \code{"egpd"} or any type of CDF (see \sQuote{Details}). \code{"egpd"} uses the extended GPD for backtransformation. CDF allows for specifying any distribution \sQuote{Examples}.}
\item{verbose}{logical. Should progress be reported?}
\item{n_wave}{number of scales to be considered in the continuous wavelet transform.}
\item{...}{any other argument passed to the sub-function specifying the cdf for fitting. See \sQuote{Details} and \sQuote{Examples}. }
}
\details{
Time can be given with three columns named \code{"YYYY"}, \code{"MM"}, \code{"DD"}, or as in POSIXct format \code{YYYY-MM-DD}.
All leap days (Feb 29th) will be omitted from the analysis, but no missing observations are allowed.
Stations are identified by list index.
Alternative distributions can be specified by providing three functions: (1) a function fitting the parameters of a distributions and providing a vector of these parameters as output (CDF_fit), (2) a function simulating random numbers from this distribution (rCDF), and (3) a function specifying the distribution (pCDF). See \sQuote{Examples} for the generalized beta for the second kind and for the Generalized Extreme Values (GEV) distribution.
}
\value{A list with elements temperature and precipitation of
\item{simulation}{data frames with time information, observations, and
\code{number_sim} columns containing the simulated data.}
}
\references{
Brunner, M. I., and E. Gilleland (2021). Spatial compound hot-dry events in the United States: assessment using a multi-site multi-variable weather generator, in preparation.}
\author{
Manuela Brunner
}
\examples{
data(weather_multi_sites)
\dontrun{ # The following call requires half minute or so to execute.
prsim.weather(data_p=data_p, data_t=data_t, number_sim=1, p_margin='egpd',t_margin='sep')
}
\dontrun{ # The following call requires 5 seconds to execute
### define normal distribution
library(fitdistrplus)
rNORM <- function(n, theta) rnorm(n, theta[1], theta[2])
pNORM <- function(x, theta) pnorm(x, theta[1], theta[2])
NORM_fit <- function( xdat, ...) fitdistr( xdat, 'normal', show=FALSE, ...)$estimate
### define GEV distribution
require("evd")
require("ismev")
rGEV <- function(n, theta) rgev(n, theta[1], theta[2], theta[3])
pGEV <- function(x, theta) pgev(x, theta[1], theta[2], theta[3])
GEV_fit <- function( xdat, ...) gev.fit(xdat, show=FALSE, ...)$mle
### apply function using alternative distributions
out <- prsim.weather(data_p=data_p, data_t=data_t, number_sim=1,p_margin='GEV',t_margin='NORM')
}
}
\keyword{ts}
\name{weather_multi_sites}
\alias{weather_multi_sites}
\alias{weather multi sites}
\docType{data}
\title{
Sample temperature and precipitation of four catchments derived from the ERA5-Land gridded dataset
}
\description{
Reanalysis data of four grid cells from ERA5-Land.
}
\usage{data("weather_multi_sites")}
\format{
Contains two lists data_p and data_t containing precipitation and temperature
data, respectively.
Each list consists of four data frames (one list per station/grid cell) of the following 4 variables.
\describe{
\item{\code{YYYY}}{a numeric vector, year}
\item{\code{MM}}{a numeric vector, month}
\item{\code{DD}}{a numeric vector, day}
\item{\code{Precip/Temp}}{a numeric vector, observed precipitation/temperature}
}
}
\details{
The data contains data for four grid cells in the Pacific Northwest.
}
\source{
The preciptation data were downloaded from ERA5-Land
\url{https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land?tab=overview}.
}
\references{
Brunner, M. I., and E. Gilleland (2021). Spatial compound hot-dry events in the United States: assessment using a multi-site multi-variable weather generator, in preparation.
}
\examples{
data(weater_multi_sites)
str(weather_multi_sites)
weather_multi_sites[[1]]$timestamp <- paste(weather_multi_sites[[1]]$YYYY,
weather_multi_sites[[1]]$MM, weather_multi_sites[[1]]$DD, sep=" ")
weather_multi_sites[[1]]$timestamp <-
as.POSIXct(strptime(weather_multi_sites[[1]]$timestamp,format="\%Y \%m \%d", tz="GMT"))
plot(weather_multi_sites[[1]]$timestamp[1:1000], weather_multi_sites[[1]]$Qobs[1:1000], type="l",
xlab="Time [d]", ylab=expression(paste("Temperature [degrees]")))
}
\keyword{datasets}
\name{weather_sim_multi_sites}
\alias{weather.sim.multi.sites}
\alias{weather_sim_multi_sites}
\docType{data}
\title{
Simulated temperature and precipitation for four grid cells
}
\description{
The dataset is generated with the package own routines and represent 5 series of 38 years of meteorological data for four grid cells
}
\usage{data("weather_sim_multi_sites")}
\format{
Two lists (one per variable) of four elements (one per catchment), containing a data frame each holding information about the observed time series and the stochastic simulations
\describe{
\item{\code{YYYY}}{a numeric vector, year}
\item{\code{MM}}{a numeric vector, month}
\item{\code{DD}}{a numeric vector, day}
\item{\code{timestamp}}{\code{POSIXct} vector of the daily runoff}
\item{\code{Prec/Temp}}{observed precipitation/temperature}
\item{\code{r1},\dots,\code{r5}}{5 simulated data series}
}
}
\details{
The data is included to illustrate the validation and visualization routines in \code{demo("PRSim_weather-validate")}.
}
\source{
The data has been generated with
\code{prsim.weather(data_p=data_p, data_t=data_t, number_sim=5, p_margin='egpd',t_margin='sep')}
(default values for all other arguments).
}
\references{
Brunner, M. I., and E. Gilleland (2021). Spatial compound hot-dry events in the United States: assessment using a multi-site multi-variable weather generator, in preparation.
}
\examples{
### define plotting colors
col_sim <- adjustcolor("#fd8d3c",alpha=0.8)
col_sim_tran <- adjustcolor("#fd8d3c",alpha=0.2)
col_obs <- adjustcolor( "black", alpha.f = 0.2)
### greys
col_vect_obs <- c('#cccccc','#969696','#636363','#252525')
### oranges
col_vect_sim <- c('#fdbe85','#fd8d3c','#e6550d','#a63603')
### plot time series for multiple sites
### Temperature (first list entry)
par(mfrow=c(2,1),mar=c(3,3,2,1))
### determine ylim
ylim_max <- max(sim[[1]][[1]]$Temp)*1.5
### observed
plot(sim[[1]][[1]]$Temp[1:1000],ylab=expression(bold(paste("Temperature [degrees]"))),xlab="Time [d]",type="l",col=col_vect_obs[1],ylim=c(0,ylim_max),main='Observations')
for(l in 2:4){
lines(sim[[l]][[1]]$Temp[1:1000],col=col_vect_obs[l])
}
# legend('topleft',legend=c('Station 1','Station 2','Station 3','Station 4'),lty=1,col=col_vect_obs[1:4])
### simulated (one run)
plot(sim[[1]][[1]]$r1[1:1000],ylab=expression(bold(paste("Temperature [degrees]"))),xlab="Time [d]",type="l",col=col_vect_sim[1],ylim=c(0,ylim_max),main='Stochastic simulations')
for(l in 2:4){
lines(sim[[l]][[1]]$r1[1:1000],col=col_vect_sim[l])
}
### precipitation (second list entry)
ylim_max <- max(sim[[1]][[2]]$Prec)*1
### observed
plot(sim[[1]][[2]]$Prec[1:1000],ylab=expression(bold(paste("Precipitation [mm/d]"))),xlab="Time [d]",type="l",col=col_vect_obs[1],ylim=c(0,ylim_max),main='Observations')
for(l in 2:4){
lines(sim[[l]][[2]]$Prec[1:1000],col=col_vect_obs[l])
}
# legend('topleft',legend=c('Station 1','Station 2','Station 3','Station 4'),lty=1,col=col_vect_obs[1:4])
### simulated (one run)
plot(sim[[1]][[2]]$r1[1:1000],ylab=expression(bold(paste("Precipitation [mm/d]"))),xlab="Time [d]",type="l",col=col_vect_sim[1],ylim=c(0,ylim_max),main='Stochastic simulations')
for(l in 2:4){
lines(sim[[l]][[2]]$r1[1:1000],col=col_vect_sim[l])
}
}
\keyword{datasets}
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