In HydrosystemsGroup/Weather-Generator: Daily weather generator for performing a climate stress test.
Load Libraries
library(forecast)
library(dplyr)
library(tidyr)
library(stringr)
library(lubridate)
library(ggplot2)
theme_set(theme_bw())
library(weathergen)
data(maurer)
# set RNG seed for reproducibility
set.seed(2345)
Load Monthly Climate Data
This analysis will use the gridded climate dataset by Maurer et al, 2002, which is provided in the maurer dataset in this package. This dataset was created using the makefile and R scripts in the /scripts/maurer_mon directory of the package.
Here is a map of the climate data grid points.
df.grid <- maurer$grid
library(ggmap)
map <- get_map(location=c(lon=mean(range(df.grid$LON)),
lat=mean(range(df.grid$LAT))),
zoom=6, maptype="satellite", color='bw')
ggmap(map, darken=c(0.25, "white"), extent="device") +
geom_point(aes(x=LON, y=LAT), data=df.grid, color='red', size=1)
Zoomed into new england. The green points will be used to compute a regional average.
ggmap(map, darken=c(0.25, "white"), extent="device") +
geom_point(aes(x=LON, y=LAT, color=SELECT),
data=mutate(df.grid, SELECT=(LAT<=43)), size=1) +
geom_hline(yintercept=43, color='green', linetype=2) +
scale_color_manual(values=c("TRUE"="green", "FALSE"="red"), guide=FALSE)
Monthly Climate Dataset
Using the selected grid points, retrieve the monthly timeseries from the database. Also, compute the variable TAVG as the arithmetic mean of TMIN and TMAX.
clim.mon <- maurer$data %>%
filter(LAT<=43) %>%
collect %>%
select(-WIND) %>%
mutate(TAVG=(TMIN+TMAX)/2) %>%
gather(VAR, VALUE, PRCP, TMIN, TMAX, TAVG)
summary(clim.mon)
Compute annual timeseries of each climate variable and each location. Note that precipitation is simply the annual sum, while the other three variables are computed as weighted means based on the number of days in each month to get a correct annual mean.
clim.yr <- spread(clim.mon, VAR, VALUE) %>%
mutate(MONTHYEAR=ymd(paste(YEAR,MONTH,1,sep='-')),
N_DAY=days_in_month(MONTHYEAR)) %>%
group_by(LAT, LON, YEAR) %>%
summarise(PRCP=sum(PRCP),
TMAX=sum(TMAX*N_DAY)/sum(N_DAY),
TMIN=sum(TMIN*N_DAY)/sum(N_DAY),
TAVG=sum(TAVG*N_DAY)/sum(N_DAY)) %>%
gather(VAR, VALUE, PRCP:TAVG)
summary(clim.yr)
Regional Average
Compute the regional average monthly value of each climate variable.
clim.reg.mon <- clim.mon %>%
group_by(YEAR, MONTH, DATE, VAR) %>%
summarise(MEAN=mean(VALUE),
SD=sd(VALUE))
summary(clim.reg.mon)
This figure shows the regional average monthly climate timeseries. The shaded region shows mean +/- 1 standard deviation.
clim.reg.mon %>%
ggplot(aes(DATE, MEAN)) +
geom_ribbon(aes(ymin=MEAN-SD, ymax=MEAN+SD), fill='grey50') +
geom_line(color='blue') +
facet_wrap(~VAR, ncol=1, scales='free_y') +
labs(x="Month/Year", y="Mean +/- StDev")
Compute the regional average annual value of each climate variable.
clim.reg.yr <- clim.yr %>%
group_by(VAR, YEAR) %>%
summarise(MEAN=mean(VALUE),
SD=sd(VALUE))
summary(clim.reg.yr)
This figure shows the regional-average annual climate timeseries. The shaded region shows the mean +/- 1 standard deviation.
ggplot(clim.reg.yr, aes(YEAR, MEAN)) +
geom_ribbon(aes(ymin=MEAN-SD, ymax=MEAN+SD), fill='grey80') +
geom_line(color='blue') +
facet_wrap(~VAR, ncol=1, scales='free_y') +
labs(x="Year", y="Mean +/- StDev")
Convert the regional average datasets wide format for use in the weather generator.
clim.reg.mon <- select(clim.reg.mon, -SD) %>%
ungroup %>%
spread(VAR, MEAN)
head(clim.reg.mon)
clim.reg.yr <- select(clim.reg.yr, -SD) %>%
ungroup %>%
spread(VAR, MEAN)
head(clim.reg.yr)
Wavelet Analysis
Perform wavelet analysis on regional average annual precipitation timeseries.
wt.reg <- wavelet_analysis(clim.reg.yr$PRCP, years=clim.reg.yr$YEAR, sig.level=0.90, noise.type='white')
par(mfrow=c(1,2))
plot(wt.reg, plot.cb=TRUE, plot.phase=FALSE)
plot(wt.reg$gws, wt.reg$period, type="b",
xlab="Global Wavelet Spectrum", ylab="",
log="y", ylim=rev(range(wt.reg$period)),
xlim=range(c(wt.reg$gws, wt.reg$gws.sig$signif)))
lines(wt.reg$gws.sig$signif, wt.reg$period, lty=2, col='red')
Although there are some significant wavelet periods, we'll use a simple arima model of the regional annual precipitation instead of arima models of the wavelet components.
models <- list(PRCP=auto.arima(clim.reg.yr$PRCP,
max.p=2,max.q=2,max.P=0,max.Q=0,
stationary=TRUE))
models[["PRCP"]]
This figure shows a 20-year forecast of the arima model reflecting no low-frequency oscillations (since we are not using the wavelet decomposition).
par(mfrow=c(1,1))
plot(forecast(models[['PRCP']], h=20), main="ARIMA Forecast of Regional Annual Precip",
xlab="Timestep (yr)", ylab='Annual Precip (mm/yr)')
Generate 50 simulations of annual precipitation using the AR1 model of the regional annual precipitation.
sim.prcp <- weathergen::arimas_monte(models=models, n=nrow(clim.reg.yr), n.iter=50)
str(sim.prcp)
Compare the GWS power and the mean, standard deviation, and skewness statistics across these simulations (red point is the overall mean, boxplots are the distribution across simulations) to regional annual precipitation timeseries (blue point).
par(mfrow=c(2,2))
x1 <- 1
x2 <- min(length(wt.reg$period), nrow(sim.prcp$gws))
ymin <- min(wt.reg$period[x1:x2], wt.reg$gws[x1:x2],
sim.prcp$gws.stat[x1:x2, 'MEAN'],
sim.prcp$gws.stat[x1:x2, 'Q025'])
ymax <- max(wt.reg$period[x1:x2], wt.reg$gws[x1:x2],
sim.prcp$gws.stat[x1:x2, 'MEAN'],
sim.prcp$gws.stat[x1:x2, 'Q975'])
plot(wt.reg$period[x1:x2], wt.reg$gws[x1:x2],
type="n", ylim=c(ymin,ymax), xlab="Period (years)", ylab="",
main="", log="y")
mtext(side=2, expression(paste("Power (",mm^2,")")), line=2.5)
polygon(c(wt.reg$period[x1:x2],
rev(wt.reg$period[x1:x2])),
c(sim.prcp$gws.stat[x1:x2, 'Q025'],
rev(sim.prcp$gws.stat[x1:x2, 'Q975'])),
col="grey")
lines(wt.reg$period[x1:x2], wt.reg$gws[x1:x2])
lines(wt.reg$period[x1:x2], sim.prcp$gws.stat[x1:x2, 'MEAN'], lty=2, col="blue")
lines(wt.reg$period[x1:x2], wt.reg$gws.sig$signif[x1:x2], col="red", lty=3)
boxplot(colMeans(sim.prcp$x), main="Mean")
points(mean(clim.reg.yr$PRCP),col="red",pch=19)
points(mean(colMeans(sim.prcp$x)), col="blue", pch=19)
boxplot(apply(sim.prcp$x, 2, sd), main="Standard Deviation")
points(sd(clim.reg.yr$PRCP), col="red", pch=19)
points(mean(apply(sim.prcp$x, 2, sd)), col="blue", pch=19)
boxplot(apply(sim.prcp$x, 2, skewness), main="Skew")
points(skewness(clim.reg.yr$PRCP), col="red", pch=19)
points(mean(apply(sim.prcp$x, 2, skewness)),col="blue",pch=19)
This figure shows the simulated and observed regional-average annual precipitation. The simulations were generated using the AR1 model described above.
as.data.frame(sim.prcp$x) %>%
mutate(TIMESTEP=row_number()) %>%
gather(TRIAL, VALUE, -TIMESTEP) %>%
ggplot(aes(TIMESTEP, VALUE)) +
geom_line(aes(group=TRIAL), alpha=0.3) +
stat_summary(fun.y=mean, geom='line', color='red', size=1) +
geom_line(aes(x=TIMESTEP, y=PRCP), data=mutate(clim.reg.yr, TIMESTEP=row_number()), color='blue', size=1) +
labs(x="Year", y="Annual Precipitation (mm/yr)",
title="WARM Annual Precipitation Simulations\nRed Line=Mean of Simulations, Blue Line=Observd Regional Mean")
Monthly Weather Generator
Select Location
First, select a random location from the complete climate dataset.
clim.locs <- select(clim.yr, LAT, LON) %>%
unique()
loc <- clim.locs[sample(nrow(clim.locs), size=1),] %>% as.numeric
names(loc) <- c('LAT', 'LON')
loc
Now extract the monthly and annual climate datasets for this location.
clim.loc.yr <- filter(clim.yr, LAT==loc[['LAT']], LON==loc[['LON']]) %>%
spread(VAR, VALUE) %>%
arrange(YEAR)
clim.loc.mon <- filter(clim.mon, LAT==loc[['LAT']], LON==loc[['LON']]) %>%
spread(VAR, VALUE) %>%
arrange(YEAR, MONTH)
This figure shows the original climate timeseries for the regional-average dataset and the selected location.
clim.loc.mon %>%
select(-LAT, -LON) %>%
mutate(SOURCE="Location") %>%
rbind(clim.reg.mon %>% mutate(SOURCE="Region")) %>%
mutate(SOURCE=factor(SOURCE)) %>%
gather(VAR, VALUE, PRCP:TAVG) %>%
mutate(DATE=ymd(paste(YEAR, MONTH, 1, sep='-'))) %>%
ggplot(aes(DATE, VALUE, color=SOURCE)) +
geom_line() +
labs(x="Month/Year", y="") +
scale_color_manual('', values=c('steelblue', 'orangered')) +
facet_wrap(~VAR, ncol=1, scales='free_y')
Nearest Neighbors
Compute the euclidian distance between a single simulated annual precipitation value and each observed year. The simulated annual precipitation is for the first year in the first simulation trial.
j <- 1 # year number
samp <- 1 # simulation trial number
# compute distances between simulated WARM timeseries and observed regional precip
distances <- cbind(clim.reg.yr, DISTANCE=sqrt((sim.prcp$x[j,samp] - clim.reg.yr$PRCP)^2))
ggplot(distances, aes(YEAR, PRCP)) +
geom_line() +
geom_point(aes(color=DISTANCE), size=3) +
geom_hline(yintercept=sim.prcp$x[j,samp], linetype=2, color='red') +
geom_text(x=max(distances$YEAR), y=sim.prcp$x[j,samp], label="Sim", hjust=1, vjust=-0.5) +
scale_color_gradient(low='green', high='red') +
labs(x="Year", y="Annual Precipitation (mm/yr)",
title="Regional and Current Simulation Annual Precipitation")
Select the eight years from the observed regional annual simulation that are the most similar to the first simulation trial and year.
distances <- mutate(distances,
IN_TOP_8=rank(DISTANCE) %in% seq(1,8))
distances %>%
ggplot(aes(YEAR, PRCP)) +
geom_line() +
geom_point(aes(color=DISTANCE), size=3) +
geom_hline(yintercept=sim.prcp$x[j,samp], linetype=2, color='red') +
geom_point(mapping=aes(alpha=IN_TOP_8), color='black', shape=1, size=4) +
geom_text(aes(x=YEAR, y=PRCP, label=YEAR), data=filter(distances, IN_TOP_8), vjust=-1) +
geom_text(x=max(distances$YEAR), y=sim.prcp$x[j,samp], label="Sim", hjust=1, vjust=-0.5) +
scale_alpha_manual(values=c("TRUE"=1, "FALSE"=0), guide=FALSE) +
scale_color_gradient(low='green', high='red') +
labs(x="Year", y="Annual Precipitation (mm)", title="Regional Annual Precipitation with 8 Nearest Neighbors Selected")
Extract the observed regional precipitation for the candidate years and compute the sampling weights as inverse distance-squared weighted.
distances.top <- filter(distances, IN_TOP_8) %>%
mutate(WEIGHT=(1/DISTANCE)^2,
WEIGHT=WEIGHT/sum(WEIGHT))
distances.top %>%
ggplot(aes(PRCP, WEIGHT)) +
geom_point() +
geom_vline(xintercept=sim.prcp$x[j,samp], linetype=2, color='red') +
geom_text(aes(label=YEAR), hjust=-0.1) +
geom_text(x=sim.prcp$x[j,samp], y=0.4, label="Sim Precip", hjust=-0.1) +
labs(x="Annual Precipitation (mm)", y="Sample Weight")
Now we can randomly sample from the 8 nearest neighbors using the sampling weights based on the inverse distance squared. We'll do this 1000 times and compare the frequency of the samples to the sample weights.
sampled.years <- sample(distances.top$YEAR,
size=1000,
prob=distances.top$WEIGHT,
replace=TRUE)
sampled.years.tbl <- prop.table(table(sampled.years))
merge(distances.top,
data.frame(YEAR=names(sampled.years.tbl),
FREQ=as.vector(sampled.years.tbl)),
all.x=TRUE) %>%
select(YEAR, WEIGHT, FREQ) %>%
ggplot(aes(WEIGHT, FREQ)) +
geom_point(size=3) +
geom_abline(linetype=2) +
labs(x="Sample Weight", y="Sampled Frequency")
Now if we just sample for one year.
sampled.year <- sample(distances.top$YEAR, size=1, prob=distances.top$WEIGHT, replace=TRUE)
sampled.year
Then we can extract the observed monthly climate data for the sampled year from the observed dataset of the selected location.
filter(clim.loc.mon, YEAR==sampled.year)
Generate Complete Timeseries
Given the selected location, repeat the nearest neighbor process for each of the num_sim simulations and num_year years, and store the result in a 3-dimensional array.
num_sim <- ncol(sim.prcp$x)
num_year <- nrow(sim.prcp$x) # number of composited years
num_months <- num_year*12
climate_variables <- select(clim.reg.yr, -YEAR) %>%
names()
sampled.year <- array(NA, c(num_year, num_sim))
gen.loc.mon.arr <- array(NA, c(num_months, length(climate_variables), num_sim))
# loop through simulations (trials)
for (samp in 1:num_sim) {
# loop through years
for (j in 1:num_year) {
# compute distances between simulated WARM timeseries and observed regional precip
distances <- mutate(clim.reg.yr,
DISTANCE=sqrt((sim.prcp$x[j, samp] - clim.reg.yr$PRCP)^2),
IN_TOP_8=rank(DISTANCE) %in% seq(1,8))
distances.top <- filter(distances, IN_TOP_8) %>%
mutate(WEIGHT=(1/DISTANCE)^2,
WEIGHT=WEIGHT/sum(WEIGHT))
sampled.year[j,samp] <- sample(distances.top$YEAR,
size=1,
prob=distances.top$WEIGHT,
replace=FALSE)
row.idx <- ((1+12*(j-1)):(12+12*(j-1)))
gen.loc.mon.arr[row.idx, , samp] <- data.matrix(filter(clim.loc.mon,
YEAR==sampled.year[j, samp]) %>%
select(PRCP, TMAX, TMIN, TAVG))
}
}
Now convert the 3-dim array to a list of dataframes with each element in the list containing the simulated data for a single variable.
gen.loc.mon <- sapply(climate_variables, function(v) {
i.var <- which(climate_variables==v)
df <- as.data.frame(gen.loc.mon.arr[,i.var,]) %>%
mutate(TIMESTEP=row_number(),
SIM_YEAR=rep(seq(1, num_year), each=12),
SIM_MONTH=rep(seq(1, 12), times=num_year),
VAR=v) %>%
gather(TRIAL, VALUE, -TIMESTEP, -SIM_YEAR, -SIM_MONTH, -VAR) %>%
mutate(TRIAL=as.numeric(TRIAL))
return(list(df))
})
str(gen.loc.mon)
Aggregate the generated timeseries to compute annual values or each variable.
gen.loc.yr <- lapply(gen.loc.mon, function(df) {
if ('PRCP' %in% unique(df$VAR)) {
df.yr <- group_by(df, VAR, TRIAL, SIM_YEAR) %>%
summarise(VALUE=sum(VALUE))
} else {
df.yr <- group_by(df, VAR, TRIAL, SIM_YEAR) %>%
summarise(VALUE=mean(VALUE))
}
df.yr <- df.yr %>%
ungroup() %>%
as.data.frame()
return(df.yr)
})
str(gen.loc.yr)
This figure shows all the weather generator simulations compared to the observed climate dataset at the selected location.
do.call(rbind, gen.loc.yr) %>%
mutate(VAR=factor(VAR)) %>%
ggplot() +
geom_line(aes(SIM_YEAR, VALUE, group=TRIAL), alpha=0.5) +
geom_line(aes(SIM_YEAR, VALUE),
data=mutate(clim.reg.yr, SIM_YEAR=row_number()) %>%
gather(VAR, VALUE, PRCP:TAVG),
color='red') +
facet_wrap(~VAR, scales='free_y') +
labs(x="Simulation Year", y="Value", title="All Generated Datasets with Observed Location Dataset ")
This figure shows only the first weather generator simulation compared to the observed climate dataset at the selected location.
do.call(rbind, gen.loc.yr) %>%
mutate(VAR=factor(VAR)) %>%
filter(TRIAL==1) %>%
ggplot() +
geom_line(aes(SIM_YEAR, VALUE, group=TRIAL), alpha=0.5) +
geom_line(aes(SIM_YEAR, VALUE),
data=mutate(clim.reg.yr, SIM_YEAR=row_number()) %>%
gather(VAR, VALUE, PRCP:TAVG),
color='red') +
facet_wrap(~VAR, scales='free_y') +
labs(x="Simulation Year", y="Value", title="Single Generated Dataset with Observed Location Dataset ")
The following three figures compare the mean, standard deviation, and skewness of the annual precipitation timeseries between the weather generator simulations and the observed climate dataset at the selected location. The boxplots show the distribution of each statistic across the simulations, the red point shows the mean of the statistic for the simulations, and the blue point shows the statistic value for the observed annual timeseries.
do.call(rbind, gen.loc.yr) %>%
mutate(VAR=factor(VAR)) %>%
group_by(VAR, TRIAL) %>%
summarise(MEAN=mean(VALUE)) %>%
ggplot(aes(x=1, y=MEAN)) +
geom_boxplot() +
stat_summary(fun.y=mean, geom="point", color="red", size=3) +
geom_point(aes(x=1, y=MEAN),
data=gather(clim.loc.yr, VAR, VALUE, PRCP:TAVG) %>%
group_by(VAR) %>%
summarise(MEAN=mean(VALUE)),
color='blue', size=3) +
facet_wrap(~VAR, scales='free_y') +
labs(x="", y="Mean", title="Mean of Annual Precipitation for Generated (Red) and Observed (Blue)") +
theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())
do.call(rbind, gen.loc.yr) %>%
mutate(VAR=factor(VAR)) %>%
group_by(VAR, TRIAL) %>%
summarise(SD=sd(VALUE)) %>%
ggplot(aes(x=1, y=SD)) +
geom_boxplot() +
stat_summary(fun.y=mean, geom="point", color="red", size=3) +
geom_point(aes(x=1, y=SD),
data=gather(clim.loc.yr, VAR, VALUE, PRCP:TAVG) %>%
group_by(VAR) %>%
summarise(SD=sd(VALUE)),
color='blue', size=3) +
facet_wrap(~VAR, scales='free_y') +
labs(x="", y="Standard Deviation", title="StDev of Annual Precipitation for Generated (Red) and Observed (Blue)") +
theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())
do.call(rbind, gen.loc.yr) %>%
mutate(VAR=factor(VAR)) %>%
group_by(VAR, TRIAL) %>%
summarise(SKEW=skewness(VALUE)) %>%
ggplot(aes(x=1, y=SKEW)) +
geom_boxplot() +
stat_summary(fun.y=mean, geom="point", color="red", size=3) +
geom_point(aes(x=1, y=SKEW),
data=gather(clim.loc.yr, VAR, VALUE, PRCP:TAVG) %>%
group_by(VAR) %>%
summarise(SKEW=skewness(VALUE)),
color='blue', size=3) +
facet_wrap(~VAR, scales='free_y') +
labs(x="", y="Skewness", title="Skewness of Annual Precipitation for Generated (Red) and Observed (Blue)") +
theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())
weathergen Package
NOTE: The following code no longer works, the monthly simulation was removed from a previous version. I'll leave it here for reference.
knitr::opts_chunk$set(eval=FALSE)
The weathergen::gen_mon_arma() function provides a wrapper around these steps to generate multiple monthly timeseries based on the regional annual precipitation timeseries and the location monthly climate timeseries. This function returns a list object containing the following elements:
$models: list of arma models (only 1 if using ARMA model instead of WARM model)$sim.prcp.yr: output from weathergen::mc_arimas() containing annual precipitation simulations from the ARMA model$sim.mon: unnamed list of length n.iter with each element containing a data.frame with the generated monthly timeseries for the given location. The columns YEAR_OBS and PRCP_OBS_YEAR contain the year and corresponding annual precipitation sampled from the monthly timeseries of the location.
wgen_mon <- weathergen::gen_month_arma(x.reg.yr=clim.reg.yr,
x.loc.mon=clim.loc.mon,
n.iter=50, n.year=nrow(clim.reg.yr))
str(wgen_mon, max.level=1)
Here is an example of the simulated monthly timeseries for the first trial.
head(wgen_mon$sim.mon[[1]])
This figure compares the total annual precipitation for the generated monthly timeseries and the annual arma model simulation. Note that the differences between the two reflect the different annual precipitation in the sampled year for the monthly time series and the simulated year from the ARMA model.
wgen_mon$sim.mon[[1]] %>%
group_by(YEAR_SIM) %>%
summarise(PRCP=sum(PRCP),
YEAR_OBS=mean(YEAR_OBS),
PRCP_OBS=mean(PRCP_OBS_YEAR)) %>%
gather(VAR, VALUE, PRCP, PRCP_OBS) %>%
ggplot(aes(YEAR_SIM, VALUE, color=VAR)) +
geom_line() +
labs(x="Year", y="Annual Precipitation (mm)")
This figure shows the first simulation trial of the monthly climate timeseries for the selected location.
wgen_mon$sim.mon[[1]] %>%
select(YEAR_SIM, MONTH, PRCP:TAVG) %>%
gather(VAR, VALUE, PRCP:TAVG) %>%
mutate(DATE=ymd(paste(YEAR_SIM, MONTH, 1, sep='-'))) %>%
ggplot(aes(DATE, VALUE)) +
geom_line() +
facet_wrap(~VAR, scales='free_y', ncol=1) +
labs(x="Month/Year", y="")
Climate Change Trends
After generating the monthly climate timeseries, we can impose linear trends on the timeseries to perform a climate stress test.
Single Timeseries
The linear trends are incorporated using an additive temperature factor and a multiplicative precipitation factor. For this example, we'll consider an 150% increase in precipitation over the simulation period, and a 2 degC increase in temperature. Note that the same temperature factor is applied to all three temperature variables (TMIN, TMAX, TAVG). We'll impose these trends on only one of the simulated monthly timeseries (chosen at random).
trend.factors <- list(PRCP=1.5, # multiplicative trend factor
TEMP=2) # additive trend factor
# select one trial at random
sim.gen <- wgen_mon$sim.mon[[sample(length(wgen_mon$sim.mon), 1)]]
# create linear trend lines
trend.lines <- list(PRCP=seq(1, trend.factors[['PRCP']], length.out=nrow(sim.gen)),
TEMP=seq(0, trend.factors[['TEMP']], length.out=nrow(sim.gen)))
sim.trend <- mutate(sim.gen,
PRCP.TREND=PRCP * trend.lines[['PRCP']],
TMAX.TREND=TMAX + trend.lines[['TEMP']],
TMIN.TREND=TMIN + trend.lines[['TEMP']],
TAVG.TREND=TAVG + trend.lines[['TEMP']])
This figure compares the original generate timeseries to the adjusted timeseries for the given trend factors.
sim.trend %>%
mutate(PRCP.GEN=PRCP, TMAX.GEN=TMAX, TMIN.GEN=TMIN, TAVG.GEN=TAVG) %>%
gather(VAR.GROUP, VALUE, PRCP.GEN:TAVG.GEN, PRCP.TREND:TAVG.TREND) %>%
separate(VAR.GROUP, c("VAR", "GROUP")) %>%
mutate(DATE=ymd(paste(YEAR_SIM, MONTH, 1, sep='-'))) %>%
ggplot(aes(DATE, VALUE, color=GROUP)) +
geom_line() +
facet_wrap(~VAR, ncol=1, scales='free_y') +
labs(x="Month/Year", y="") +
scale_color_manual('', values=c(GEN='grey20', TREND='chartreuse3'),
labels=c(GEN='w/o trend', TREND='w/ trend'))
Multiple Trends for Stress Test
The weathergen::trends_mon() function applies a range of trend factors to the set of monthly timeseries generated by weathergen::gen_month_arma(). For each pair of temperature and precipitation factors, the function will create a single monthly timeseries by sampling a random iteration from the gen_month_arma() output.
trend.factors <- list(PRCP=seq(0.75, 1.25, by=0.05),
TEMP=seq(0, 5, by=0.5))
sim.trends <- trends_mon(sim.mon.list=wgen_mon$sim.mon,
temp.factors=trend.factors[['TEMP']],
prcp.factors=trend.factors[['PRCP']])
str(sim.trends[[1]])
Now we can merge these trends and add columns containing the trend factors for temperature (TEMP_FACTOR) and precipitaiton (PRCP_FACTOR).
df.trends <- lapply(sim.trends, function(x) {
df <- x$df
df$PRCP_FACTOR <- x$trend.factors[['PRCP']]
df$TEMP_FACTOR <- x$trend.factors[['TEMP']]
return(df)
}) %>%
do.call(rbind, .)
summary(df.trends)
This figure shows the monthly precipitation timeseries colored by trend magnitude. This figure only shows precipitation timeseries without any temperature trend. You can see some increase with the 'redder' lines that have higher trend factors.
df.trends %>%
filter(TEMP_FACTOR==0) %>%
group_by(PRCP_FACTOR) %>%
mutate(TIMESTEP=row_number()) %>%
ggplot(aes(TIMESTEP, PRCP.TREND, color=PRCP_FACTOR)) +
geom_line(alpha=0.2) +
scale_color_gradient('Precip Factor', low='green', high='red') +
labs(x="Simulation Month", y="Monthly Precipitation (mm)")
This figure shows the generated precipitation timeseries over the complete range of precipitation trend factors including only timeseries with no temperature trend (TEMP_FACTOR==0). This figure clearly shows the increasing trends with higher precipitation trend factors.
df.trends %>%
filter(TEMP_FACTOR==0) %>%
group_by(TEMP_FACTOR, PRCP_FACTOR, YEAR_SIM) %>%
summarise(PRCP=sum(PRCP.TREND)) %>%
ggplot(aes(YEAR_SIM, PRCP, color=PRCP_FACTOR, group=PRCP_FACTOR)) +
geom_line(alpha=0.7) +
scale_color_gradient('Precip Factor', low='green', high='red') +
labs(x="Simulation Year", y="Annual Precipitation (mm)")
This figure shows the generated average temperature timeseries over the range of trend factors including only timeseries with no precipitation trend (PRCP_FACTOR==0). This figure clearly shows the increasing trends with higher temperature trend factors.
df.trends %>%
filter(PRCP_FACTOR==1) %>%
mutate(DATE=ymd(paste(2001, MONTH, 1, sep='-'))) %>%
group_by(TEMP_FACTOR, PRCP_FACTOR, YEAR_SIM) %>%
summarise(TAVG=sum(TAVG.TREND*days_in_month(DATE))/sum(days_in_month(DATE))) %>%
ggplot(aes(YEAR_SIM, TAVG, color=TEMP_FACTOR, group=TEMP_FACTOR)) +
geom_line(alpha=0.7) +
scale_color_gradient('Precip Factor', low='green', high='red') +
labs(x="Simulation Year", y="Annual Mean Temperature (degC)")
Session Info
sessionInfo()
HydrosystemsGroup/Weather-Generator documentation built on May 8, 2019, 8:34 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Load Libraries
library(forecast) library(dplyr) library(tidyr) library(stringr) library(lubridate) library(ggplot2) theme_set(theme_bw()) library(weathergen) data(maurer) # set RNG seed for reproducibility set.seed(2345)
Load Monthly Climate Data
This analysis will use the gridded climate dataset by Maurer et al, 2002, which is provided in the maurer dataset in this package. This dataset was created using the makefile and R scripts in the /scripts/maurer_mon directory of the package.
Here is a map of the climate data grid points.
df.grid <- maurer$grid library(ggmap) map <- get_map(location=c(lon=mean(range(df.grid$LON)), lat=mean(range(df.grid$LAT))), zoom=6, maptype="satellite", color='bw') ggmap(map, darken=c(0.25, "white"), extent="device") + geom_point(aes(x=LON, y=LAT), data=df.grid, color='red', size=1)
Zoomed into new england. The green points will be used to compute a regional average.
ggmap(map, darken=c(0.25, "white"), extent="device") + geom_point(aes(x=LON, y=LAT, color=SELECT), data=mutate(df.grid, SELECT=(LAT<=43)), size=1) + geom_hline(yintercept=43, color='green', linetype=2) + scale_color_manual(values=c("TRUE"="green", "FALSE"="red"), guide=FALSE)
Monthly Climate Dataset
Using the selected grid points, retrieve the monthly timeseries from the database. Also, compute the variable TAVG as the arithmetic mean of TMIN and TMAX.
clim.mon <- maurer$data %>% filter(LAT<=43) %>% collect %>% select(-WIND) %>% mutate(TAVG=(TMIN+TMAX)/2) %>% gather(VAR, VALUE, PRCP, TMIN, TMAX, TAVG) summary(clim.mon)
Compute annual timeseries of each climate variable and each location. Note that precipitation is simply the annual sum, while the other three variables are computed as weighted means based on the number of days in each month to get a correct annual mean.
clim.yr <- spread(clim.mon, VAR, VALUE) %>% mutate(MONTHYEAR=ymd(paste(YEAR,MONTH,1,sep='-')), N_DAY=days_in_month(MONTHYEAR)) %>% group_by(LAT, LON, YEAR) %>% summarise(PRCP=sum(PRCP), TMAX=sum(TMAX*N_DAY)/sum(N_DAY), TMIN=sum(TMIN*N_DAY)/sum(N_DAY), TAVG=sum(TAVG*N_DAY)/sum(N_DAY)) %>% gather(VAR, VALUE, PRCP:TAVG) summary(clim.yr)
Regional Average
Compute the regional average monthly value of each climate variable.
clim.reg.mon <- clim.mon %>% group_by(YEAR, MONTH, DATE, VAR) %>% summarise(MEAN=mean(VALUE), SD=sd(VALUE)) summary(clim.reg.mon)
This figure shows the regional average monthly climate timeseries. The shaded region shows mean +/- 1 standard deviation.
clim.reg.mon %>% ggplot(aes(DATE, MEAN)) + geom_ribbon(aes(ymin=MEAN-SD, ymax=MEAN+SD), fill='grey50') + geom_line(color='blue') + facet_wrap(~VAR, ncol=1, scales='free_y') + labs(x="Month/Year", y="Mean +/- StDev")
Compute the regional average annual value of each climate variable.
clim.reg.yr <- clim.yr %>% group_by(VAR, YEAR) %>% summarise(MEAN=mean(VALUE), SD=sd(VALUE)) summary(clim.reg.yr)
This figure shows the regional-average annual climate timeseries. The shaded region shows the mean +/- 1 standard deviation.
ggplot(clim.reg.yr, aes(YEAR, MEAN)) + geom_ribbon(aes(ymin=MEAN-SD, ymax=MEAN+SD), fill='grey80') + geom_line(color='blue') + facet_wrap(~VAR, ncol=1, scales='free_y') + labs(x="Year", y="Mean +/- StDev")
Convert the regional average datasets wide format for use in the weather generator.
clim.reg.mon <- select(clim.reg.mon, -SD) %>% ungroup %>% spread(VAR, MEAN) head(clim.reg.mon) clim.reg.yr <- select(clim.reg.yr, -SD) %>% ungroup %>% spread(VAR, MEAN) head(clim.reg.yr)
Wavelet Analysis
Perform wavelet analysis on regional average annual precipitation timeseries.
wt.reg <- wavelet_analysis(clim.reg.yr$PRCP, years=clim.reg.yr$YEAR, sig.level=0.90, noise.type='white') par(mfrow=c(1,2)) plot(wt.reg, plot.cb=TRUE, plot.phase=FALSE) plot(wt.reg$gws, wt.reg$period, type="b", xlab="Global Wavelet Spectrum", ylab="", log="y", ylim=rev(range(wt.reg$period)), xlim=range(c(wt.reg$gws, wt.reg$gws.sig$signif))) lines(wt.reg$gws.sig$signif, wt.reg$period, lty=2, col='red')
Although there are some significant wavelet periods, we'll use a simple arima model of the regional annual precipitation instead of arima models of the wavelet components.
models <- list(PRCP=auto.arima(clim.reg.yr$PRCP, max.p=2,max.q=2,max.P=0,max.Q=0, stationary=TRUE)) models[["PRCP"]]
This figure shows a 20-year forecast of the arima model reflecting no low-frequency oscillations (since we are not using the wavelet decomposition).
par(mfrow=c(1,1)) plot(forecast(models[['PRCP']], h=20), main="ARIMA Forecast of Regional Annual Precip", xlab="Timestep (yr)", ylab='Annual Precip (mm/yr)')
Generate 50 simulations of annual precipitation using the AR1 model of the regional annual precipitation.
sim.prcp <- weathergen::arimas_monte(models=models, n=nrow(clim.reg.yr), n.iter=50) str(sim.prcp)
Compare the GWS power and the mean, standard deviation, and skewness statistics across these simulations (red point is the overall mean, boxplots are the distribution across simulations) to regional annual precipitation timeseries (blue point).
par(mfrow=c(2,2)) x1 <- 1 x2 <- min(length(wt.reg$period), nrow(sim.prcp$gws)) ymin <- min(wt.reg$period[x1:x2], wt.reg$gws[x1:x2], sim.prcp$gws.stat[x1:x2, 'MEAN'], sim.prcp$gws.stat[x1:x2, 'Q025']) ymax <- max(wt.reg$period[x1:x2], wt.reg$gws[x1:x2], sim.prcp$gws.stat[x1:x2, 'MEAN'], sim.prcp$gws.stat[x1:x2, 'Q975']) plot(wt.reg$period[x1:x2], wt.reg$gws[x1:x2], type="n", ylim=c(ymin,ymax), xlab="Period (years)", ylab="", main="", log="y") mtext(side=2, expression(paste("Power (",mm^2,")")), line=2.5) polygon(c(wt.reg$period[x1:x2], rev(wt.reg$period[x1:x2])), c(sim.prcp$gws.stat[x1:x2, 'Q025'], rev(sim.prcp$gws.stat[x1:x2, 'Q975'])), col="grey") lines(wt.reg$period[x1:x2], wt.reg$gws[x1:x2]) lines(wt.reg$period[x1:x2], sim.prcp$gws.stat[x1:x2, 'MEAN'], lty=2, col="blue") lines(wt.reg$period[x1:x2], wt.reg$gws.sig$signif[x1:x2], col="red", lty=3) boxplot(colMeans(sim.prcp$x), main="Mean") points(mean(clim.reg.yr$PRCP),col="red",pch=19) points(mean(colMeans(sim.prcp$x)), col="blue", pch=19) boxplot(apply(sim.prcp$x, 2, sd), main="Standard Deviation") points(sd(clim.reg.yr$PRCP), col="red", pch=19) points(mean(apply(sim.prcp$x, 2, sd)), col="blue", pch=19) boxplot(apply(sim.prcp$x, 2, skewness), main="Skew") points(skewness(clim.reg.yr$PRCP), col="red", pch=19) points(mean(apply(sim.prcp$x, 2, skewness)),col="blue",pch=19)
This figure shows the simulated and observed regional-average annual precipitation. The simulations were generated using the AR1 model described above.
as.data.frame(sim.prcp$x) %>% mutate(TIMESTEP=row_number()) %>% gather(TRIAL, VALUE, -TIMESTEP) %>% ggplot(aes(TIMESTEP, VALUE)) + geom_line(aes(group=TRIAL), alpha=0.3) + stat_summary(fun.y=mean, geom='line', color='red', size=1) + geom_line(aes(x=TIMESTEP, y=PRCP), data=mutate(clim.reg.yr, TIMESTEP=row_number()), color='blue', size=1) + labs(x="Year", y="Annual Precipitation (mm/yr)", title="WARM Annual Precipitation Simulations\nRed Line=Mean of Simulations, Blue Line=Observd Regional Mean")
Monthly Weather Generator
Select Location
First, select a random location from the complete climate dataset.
clim.locs <- select(clim.yr, LAT, LON) %>% unique() loc <- clim.locs[sample(nrow(clim.locs), size=1),] %>% as.numeric names(loc) <- c('LAT', 'LON') loc
Now extract the monthly and annual climate datasets for this location.
clim.loc.yr <- filter(clim.yr, LAT==loc[['LAT']], LON==loc[['LON']]) %>% spread(VAR, VALUE) %>% arrange(YEAR) clim.loc.mon <- filter(clim.mon, LAT==loc[['LAT']], LON==loc[['LON']]) %>% spread(VAR, VALUE) %>% arrange(YEAR, MONTH)
This figure shows the original climate timeseries for the regional-average dataset and the selected location.
clim.loc.mon %>% select(-LAT, -LON) %>% mutate(SOURCE="Location") %>% rbind(clim.reg.mon %>% mutate(SOURCE="Region")) %>% mutate(SOURCE=factor(SOURCE)) %>% gather(VAR, VALUE, PRCP:TAVG) %>% mutate(DATE=ymd(paste(YEAR, MONTH, 1, sep='-'))) %>% ggplot(aes(DATE, VALUE, color=SOURCE)) + geom_line() + labs(x="Month/Year", y="") + scale_color_manual('', values=c('steelblue', 'orangered')) + facet_wrap(~VAR, ncol=1, scales='free_y')
Nearest Neighbors
Compute the euclidian distance between a single simulated annual precipitation value and each observed year. The simulated annual precipitation is for the first year in the first simulation trial.
j <- 1 # year number samp <- 1 # simulation trial number # compute distances between simulated WARM timeseries and observed regional precip distances <- cbind(clim.reg.yr, DISTANCE=sqrt((sim.prcp$x[j,samp] - clim.reg.yr$PRCP)^2)) ggplot(distances, aes(YEAR, PRCP)) + geom_line() + geom_point(aes(color=DISTANCE), size=3) + geom_hline(yintercept=sim.prcp$x[j,samp], linetype=2, color='red') + geom_text(x=max(distances$YEAR), y=sim.prcp$x[j,samp], label="Sim", hjust=1, vjust=-0.5) + scale_color_gradient(low='green', high='red') + labs(x="Year", y="Annual Precipitation (mm/yr)", title="Regional and Current Simulation Annual Precipitation")
Select the eight years from the observed regional annual simulation that are the most similar to the first simulation trial and year.
distances <- mutate(distances, IN_TOP_8=rank(DISTANCE) %in% seq(1,8)) distances %>% ggplot(aes(YEAR, PRCP)) + geom_line() + geom_point(aes(color=DISTANCE), size=3) + geom_hline(yintercept=sim.prcp$x[j,samp], linetype=2, color='red') + geom_point(mapping=aes(alpha=IN_TOP_8), color='black', shape=1, size=4) + geom_text(aes(x=YEAR, y=PRCP, label=YEAR), data=filter(distances, IN_TOP_8), vjust=-1) + geom_text(x=max(distances$YEAR), y=sim.prcp$x[j,samp], label="Sim", hjust=1, vjust=-0.5) + scale_alpha_manual(values=c("TRUE"=1, "FALSE"=0), guide=FALSE) + scale_color_gradient(low='green', high='red') + labs(x="Year", y="Annual Precipitation (mm)", title="Regional Annual Precipitation with 8 Nearest Neighbors Selected")
Extract the observed regional precipitation for the candidate years and compute the sampling weights as inverse distance-squared weighted.
distances.top <- filter(distances, IN_TOP_8) %>% mutate(WEIGHT=(1/DISTANCE)^2, WEIGHT=WEIGHT/sum(WEIGHT)) distances.top %>% ggplot(aes(PRCP, WEIGHT)) + geom_point() + geom_vline(xintercept=sim.prcp$x[j,samp], linetype=2, color='red') + geom_text(aes(label=YEAR), hjust=-0.1) + geom_text(x=sim.prcp$x[j,samp], y=0.4, label="Sim Precip", hjust=-0.1) + labs(x="Annual Precipitation (mm)", y="Sample Weight")
Now we can randomly sample from the 8 nearest neighbors using the sampling weights based on the inverse distance squared. We'll do this 1000 times and compare the frequency of the samples to the sample weights.
sampled.years <- sample(distances.top$YEAR, size=1000, prob=distances.top$WEIGHT, replace=TRUE) sampled.years.tbl <- prop.table(table(sampled.years)) merge(distances.top, data.frame(YEAR=names(sampled.years.tbl), FREQ=as.vector(sampled.years.tbl)), all.x=TRUE) %>% select(YEAR, WEIGHT, FREQ) %>% ggplot(aes(WEIGHT, FREQ)) + geom_point(size=3) + geom_abline(linetype=2) + labs(x="Sample Weight", y="Sampled Frequency")
Now if we just sample for one year.
sampled.year <- sample(distances.top$YEAR, size=1, prob=distances.top$WEIGHT, replace=TRUE) sampled.year
Then we can extract the observed monthly climate data for the sampled year from the observed dataset of the selected location.
filter(clim.loc.mon, YEAR==sampled.year)
Generate Complete Timeseries
Given the selected location, repeat the nearest neighbor process for each of the num_sim simulations and num_year years, and store the result in a 3-dimensional array.
num_sim <- ncol(sim.prcp$x) num_year <- nrow(sim.prcp$x) # number of composited years num_months <- num_year*12 climate_variables <- select(clim.reg.yr, -YEAR) %>% names() sampled.year <- array(NA, c(num_year, num_sim)) gen.loc.mon.arr <- array(NA, c(num_months, length(climate_variables), num_sim))
# loop through simulations (trials) for (samp in 1:num_sim) { # loop through years for (j in 1:num_year) { # compute distances between simulated WARM timeseries and observed regional precip distances <- mutate(clim.reg.yr, DISTANCE=sqrt((sim.prcp$x[j, samp] - clim.reg.yr$PRCP)^2), IN_TOP_8=rank(DISTANCE) %in% seq(1,8)) distances.top <- filter(distances, IN_TOP_8) %>% mutate(WEIGHT=(1/DISTANCE)^2, WEIGHT=WEIGHT/sum(WEIGHT)) sampled.year[j,samp] <- sample(distances.top$YEAR, size=1, prob=distances.top$WEIGHT, replace=FALSE) row.idx <- ((1+12*(j-1)):(12+12*(j-1))) gen.loc.mon.arr[row.idx, , samp] <- data.matrix(filter(clim.loc.mon, YEAR==sampled.year[j, samp]) %>% select(PRCP, TMAX, TMIN, TAVG)) } }
Now convert the 3-dim array to a list of dataframes with each element in the list containing the simulated data for a single variable.
gen.loc.mon <- sapply(climate_variables, function(v) { i.var <- which(climate_variables==v) df <- as.data.frame(gen.loc.mon.arr[,i.var,]) %>% mutate(TIMESTEP=row_number(), SIM_YEAR=rep(seq(1, num_year), each=12), SIM_MONTH=rep(seq(1, 12), times=num_year), VAR=v) %>% gather(TRIAL, VALUE, -TIMESTEP, -SIM_YEAR, -SIM_MONTH, -VAR) %>% mutate(TRIAL=as.numeric(TRIAL)) return(list(df)) }) str(gen.loc.mon)
Aggregate the generated timeseries to compute annual values or each variable.
gen.loc.yr <- lapply(gen.loc.mon, function(df) { if ('PRCP' %in% unique(df$VAR)) { df.yr <- group_by(df, VAR, TRIAL, SIM_YEAR) %>% summarise(VALUE=sum(VALUE)) } else { df.yr <- group_by(df, VAR, TRIAL, SIM_YEAR) %>% summarise(VALUE=mean(VALUE)) } df.yr <- df.yr %>% ungroup() %>% as.data.frame() return(df.yr) }) str(gen.loc.yr)
This figure shows all the weather generator simulations compared to the observed climate dataset at the selected location.
do.call(rbind, gen.loc.yr) %>% mutate(VAR=factor(VAR)) %>% ggplot() + geom_line(aes(SIM_YEAR, VALUE, group=TRIAL), alpha=0.5) + geom_line(aes(SIM_YEAR, VALUE), data=mutate(clim.reg.yr, SIM_YEAR=row_number()) %>% gather(VAR, VALUE, PRCP:TAVG), color='red') + facet_wrap(~VAR, scales='free_y') + labs(x="Simulation Year", y="Value", title="All Generated Datasets with Observed Location Dataset ")
This figure shows only the first weather generator simulation compared to the observed climate dataset at the selected location.
do.call(rbind, gen.loc.yr) %>% mutate(VAR=factor(VAR)) %>% filter(TRIAL==1) %>% ggplot() + geom_line(aes(SIM_YEAR, VALUE, group=TRIAL), alpha=0.5) + geom_line(aes(SIM_YEAR, VALUE), data=mutate(clim.reg.yr, SIM_YEAR=row_number()) %>% gather(VAR, VALUE, PRCP:TAVG), color='red') + facet_wrap(~VAR, scales='free_y') + labs(x="Simulation Year", y="Value", title="Single Generated Dataset with Observed Location Dataset ")
The following three figures compare the mean, standard deviation, and skewness of the annual precipitation timeseries between the weather generator simulations and the observed climate dataset at the selected location. The boxplots show the distribution of each statistic across the simulations, the red point shows the mean of the statistic for the simulations, and the blue point shows the statistic value for the observed annual timeseries.
do.call(rbind, gen.loc.yr) %>% mutate(VAR=factor(VAR)) %>% group_by(VAR, TRIAL) %>% summarise(MEAN=mean(VALUE)) %>% ggplot(aes(x=1, y=MEAN)) + geom_boxplot() + stat_summary(fun.y=mean, geom="point", color="red", size=3) + geom_point(aes(x=1, y=MEAN), data=gather(clim.loc.yr, VAR, VALUE, PRCP:TAVG) %>% group_by(VAR) %>% summarise(MEAN=mean(VALUE)), color='blue', size=3) + facet_wrap(~VAR, scales='free_y') + labs(x="", y="Mean", title="Mean of Annual Precipitation for Generated (Red) and Observed (Blue)") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())
do.call(rbind, gen.loc.yr) %>% mutate(VAR=factor(VAR)) %>% group_by(VAR, TRIAL) %>% summarise(SD=sd(VALUE)) %>% ggplot(aes(x=1, y=SD)) + geom_boxplot() + stat_summary(fun.y=mean, geom="point", color="red", size=3) + geom_point(aes(x=1, y=SD), data=gather(clim.loc.yr, VAR, VALUE, PRCP:TAVG) %>% group_by(VAR) %>% summarise(SD=sd(VALUE)), color='blue', size=3) + facet_wrap(~VAR, scales='free_y') + labs(x="", y="Standard Deviation", title="StDev of Annual Precipitation for Generated (Red) and Observed (Blue)") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())
do.call(rbind, gen.loc.yr) %>% mutate(VAR=factor(VAR)) %>% group_by(VAR, TRIAL) %>% summarise(SKEW=skewness(VALUE)) %>% ggplot(aes(x=1, y=SKEW)) + geom_boxplot() + stat_summary(fun.y=mean, geom="point", color="red", size=3) + geom_point(aes(x=1, y=SKEW), data=gather(clim.loc.yr, VAR, VALUE, PRCP:TAVG) %>% group_by(VAR) %>% summarise(SKEW=skewness(VALUE)), color='blue', size=3) + facet_wrap(~VAR, scales='free_y') + labs(x="", y="Skewness", title="Skewness of Annual Precipitation for Generated (Red) and Observed (Blue)") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())
weathergen Package
NOTE: The following code no longer works, the monthly simulation was removed from a previous version. I'll leave it here for reference.
knitr::opts_chunk$set(eval=FALSE)
The weathergen::gen_mon_arma() function provides a wrapper around these steps to generate multiple monthly timeseries based on the regional annual precipitation timeseries and the location monthly climate timeseries. This function returns a list object containing the following elements:
$models: list of arma models (only 1 if using ARMA model instead of WARM model)$sim.prcp.yr: output fromweathergen::mc_arimas()containing annual precipitation simulations from the ARMA model$sim.mon: unnamed list of lengthn.iterwith each element containing a data.frame with the generated monthly timeseries for the given location. The columnsYEAR_OBSandPRCP_OBS_YEARcontain the year and corresponding annual precipitation sampled from the monthly timeseries of the location.
wgen_mon <- weathergen::gen_month_arma(x.reg.yr=clim.reg.yr, x.loc.mon=clim.loc.mon, n.iter=50, n.year=nrow(clim.reg.yr)) str(wgen_mon, max.level=1)
Here is an example of the simulated monthly timeseries for the first trial.
head(wgen_mon$sim.mon[[1]])
This figure compares the total annual precipitation for the generated monthly timeseries and the annual arma model simulation. Note that the differences between the two reflect the different annual precipitation in the sampled year for the monthly time series and the simulated year from the ARMA model.
wgen_mon$sim.mon[[1]] %>% group_by(YEAR_SIM) %>% summarise(PRCP=sum(PRCP), YEAR_OBS=mean(YEAR_OBS), PRCP_OBS=mean(PRCP_OBS_YEAR)) %>% gather(VAR, VALUE, PRCP, PRCP_OBS) %>% ggplot(aes(YEAR_SIM, VALUE, color=VAR)) + geom_line() + labs(x="Year", y="Annual Precipitation (mm)")
This figure shows the first simulation trial of the monthly climate timeseries for the selected location.
wgen_mon$sim.mon[[1]] %>% select(YEAR_SIM, MONTH, PRCP:TAVG) %>% gather(VAR, VALUE, PRCP:TAVG) %>% mutate(DATE=ymd(paste(YEAR_SIM, MONTH, 1, sep='-'))) %>% ggplot(aes(DATE, VALUE)) + geom_line() + facet_wrap(~VAR, scales='free_y', ncol=1) + labs(x="Month/Year", y="")
Climate Change Trends
After generating the monthly climate timeseries, we can impose linear trends on the timeseries to perform a climate stress test.
Single Timeseries
The linear trends are incorporated using an additive temperature factor and a multiplicative precipitation factor. For this example, we'll consider an 150% increase in precipitation over the simulation period, and a 2 degC increase in temperature. Note that the same temperature factor is applied to all three temperature variables (TMIN, TMAX, TAVG). We'll impose these trends on only one of the simulated monthly timeseries (chosen at random).
trend.factors <- list(PRCP=1.5, # multiplicative trend factor TEMP=2) # additive trend factor # select one trial at random sim.gen <- wgen_mon$sim.mon[[sample(length(wgen_mon$sim.mon), 1)]] # create linear trend lines trend.lines <- list(PRCP=seq(1, trend.factors[['PRCP']], length.out=nrow(sim.gen)), TEMP=seq(0, trend.factors[['TEMP']], length.out=nrow(sim.gen))) sim.trend <- mutate(sim.gen, PRCP.TREND=PRCP * trend.lines[['PRCP']], TMAX.TREND=TMAX + trend.lines[['TEMP']], TMIN.TREND=TMIN + trend.lines[['TEMP']], TAVG.TREND=TAVG + trend.lines[['TEMP']])
This figure compares the original generate timeseries to the adjusted timeseries for the given trend factors.
sim.trend %>% mutate(PRCP.GEN=PRCP, TMAX.GEN=TMAX, TMIN.GEN=TMIN, TAVG.GEN=TAVG) %>% gather(VAR.GROUP, VALUE, PRCP.GEN:TAVG.GEN, PRCP.TREND:TAVG.TREND) %>% separate(VAR.GROUP, c("VAR", "GROUP")) %>% mutate(DATE=ymd(paste(YEAR_SIM, MONTH, 1, sep='-'))) %>% ggplot(aes(DATE, VALUE, color=GROUP)) + geom_line() + facet_wrap(~VAR, ncol=1, scales='free_y') + labs(x="Month/Year", y="") + scale_color_manual('', values=c(GEN='grey20', TREND='chartreuse3'), labels=c(GEN='w/o trend', TREND='w/ trend'))
Multiple Trends for Stress Test
The weathergen::trends_mon() function applies a range of trend factors to the set of monthly timeseries generated by weathergen::gen_month_arma(). For each pair of temperature and precipitation factors, the function will create a single monthly timeseries by sampling a random iteration from the gen_month_arma() output.
trend.factors <- list(PRCP=seq(0.75, 1.25, by=0.05), TEMP=seq(0, 5, by=0.5)) sim.trends <- trends_mon(sim.mon.list=wgen_mon$sim.mon, temp.factors=trend.factors[['TEMP']], prcp.factors=trend.factors[['PRCP']]) str(sim.trends[[1]])
Now we can merge these trends and add columns containing the trend factors for temperature (TEMP_FACTOR) and precipitaiton (PRCP_FACTOR).
df.trends <- lapply(sim.trends, function(x) { df <- x$df df$PRCP_FACTOR <- x$trend.factors[['PRCP']] df$TEMP_FACTOR <- x$trend.factors[['TEMP']] return(df) }) %>% do.call(rbind, .) summary(df.trends)
This figure shows the monthly precipitation timeseries colored by trend magnitude. This figure only shows precipitation timeseries without any temperature trend. You can see some increase with the 'redder' lines that have higher trend factors.
df.trends %>% filter(TEMP_FACTOR==0) %>% group_by(PRCP_FACTOR) %>% mutate(TIMESTEP=row_number()) %>% ggplot(aes(TIMESTEP, PRCP.TREND, color=PRCP_FACTOR)) + geom_line(alpha=0.2) + scale_color_gradient('Precip Factor', low='green', high='red') + labs(x="Simulation Month", y="Monthly Precipitation (mm)")
This figure shows the generated precipitation timeseries over the complete range of precipitation trend factors including only timeseries with no temperature trend (TEMP_FACTOR==0). This figure clearly shows the increasing trends with higher precipitation trend factors.
df.trends %>% filter(TEMP_FACTOR==0) %>% group_by(TEMP_FACTOR, PRCP_FACTOR, YEAR_SIM) %>% summarise(PRCP=sum(PRCP.TREND)) %>% ggplot(aes(YEAR_SIM, PRCP, color=PRCP_FACTOR, group=PRCP_FACTOR)) + geom_line(alpha=0.7) + scale_color_gradient('Precip Factor', low='green', high='red') + labs(x="Simulation Year", y="Annual Precipitation (mm)")
This figure shows the generated average temperature timeseries over the range of trend factors including only timeseries with no precipitation trend (PRCP_FACTOR==0). This figure clearly shows the increasing trends with higher temperature trend factors.
df.trends %>% filter(PRCP_FACTOR==1) %>% mutate(DATE=ymd(paste(2001, MONTH, 1, sep='-'))) %>% group_by(TEMP_FACTOR, PRCP_FACTOR, YEAR_SIM) %>% summarise(TAVG=sum(TAVG.TREND*days_in_month(DATE))/sum(days_in_month(DATE))) %>% ggplot(aes(YEAR_SIM, TAVG, color=TEMP_FACTOR, group=TEMP_FACTOR)) + geom_line(alpha=0.7) + scale_color_gradient('Precip Factor', low='green', high='red') + labs(x="Simulation Year", y="Annual Mean Temperature (degC)")
Session Info
sessionInfo()
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.