### Process flood data and precipitation data to annual times series selecting the season and forming Rx5d and Rx20d 
library(stringr)
library(MASS)
library(biwavelet)

movingAverage <- function(x, n=1, m=1, centered=FALSE) {
    if (centered) {before <- floor((n-1)/2)
    after  <- ceiling((n-1)/2)
    } else {before <- n-1
    after  <- 0}
    s     <- rep(0, length(x))
    count <- rep(0, length(x))
    new <- x
    count <- count + !is.na(new)
    new[is.na(new)] <- 0
    s <- s + new
    i <- 1
    while (i <= before) {
        new   <- c(rep(NA, i), x[1:(length(x)-i)])
        count <- count + !is.na(new)
        new[is.na(new)] <- 0
        s <- s + new
        i <- i+1}
    i <- 1
    while (i <= after) {
        new   <- c(x[(i+1):length(x)], rep(NA, i))
        count <- count + !is.na(new)
        new[is.na(new)] <- 0
        s <- s + new
        i <- i+1}
    count[count<=m] <- NA
    s/count
    }

indir <- "U:/Documents/papers/19C_Flood/Submitted/CP/CodeInput/InputData/"
outdir <- "U:/Documents/papers/19C_Flood/Submitted/CP/CodeInput/InputData/Output/"

### read in data
lst <- read.table(paste(indir,"filelist.txt",sep=""),header=F)
day <- read.table(paste(indir,"day.txt",sep=""))

### initialize arrays, write time column
#yearround
ap <- array(NA,dim=c(length(1780:2021),dim(lst)[1]+1)) # rx5d raw
apg <- array(NA,dim=c(length(1780:2021),dim(lst)[1]+1)) # rx5d normalized
af <- array(NA,dim=c(length(1780:2021),dim(lst)[1]+1)) # annual peak streamflow raw
afg <- array(NA,dim=c(length(1780:2021),dim(lst)[1]+1)) # annual peak normalized
ap[,1] <- 1780:2021
apg[,1] <- 1780:2021
af[,1] <- 1780:2021
afg[,1] <- 1780:2021

#seasonal
sp <- array(NA,dim=c(length(1780:2021),dim(lst)[1]+1))
spg <- array(NA,dim=c(length(1780:2021),dim(lst)[1]+1))
sf <- array(NA,dim=c(length(1780:2021),dim(lst)[1]+1))
sfg <- array(NA,dim=c(length(1780:2021),dim(lst)[1]+1))
sp[,1] <- 1780:2021
spg[,1] <- 1780:2021
sf[,1] <- 1780:2021
sfg[,1] <- 1780:2021

#yearroundsmoothed
aps <- array(NA,dim=c(length(1780:2021),dim(lst)[1]+1))
apgs <- array(NA,dim=c(length(1780:2021),dim(lst)[1]+1))
afs <- array(NA,dim=c(length(1780:2021),dim(lst)[1]+1))
afgs <- array(NA,dim=c(length(1780:2021),dim(lst)[1]+1))
aps[,1] <- 1780:2021
apgs[,1] <- 1780:2021
afs[,1] <- 1780:2021
afgs[,1] <- 1780:2021

#seasonalsmoothed
sps <- array(NA,dim=c(length(1780:2021),dim(lst)[1]+1))
spgs <- array(NA,dim=c(length(1780:2021),dim(lst)[1]+1))
sfs <- array(NA,dim=c(length(1780:2021),dim(lst)[1]+1))
sfgs <- array(NA,dim=c(length(1780:2021),dim(lst)[1]+1))
sps[,1] <- 1780:2021
spgs[,1] <- 1780:2021
sfs[,1] <- 1780:2021
sfgs[,1] <- 1780:2021

### main loop
for (i in 1:length(lst[,1])){

## process precipitation
 x <- read.table(paste(indir,lst[i,2],sep=""),header=F)
 ner <- x[,4]<0 ## remove negative precip, NA 
 ner[is.na(ner)]<- T
 x[ner,4] <- NA
 if (i==21){x[x[,1]<1879,4] <- NA} ### Hohenpeissenberg is inhomogeneous before 1879
## yearround
 xa <- rep(NA,length(day[,1]))
 xa[match(10000*x[,1]+100*x[,2]+x[,3],10000*day[,1]+100*day[,2]+day[,3])]<-x[,4] ### continuous daily series
 xm <- as.matrix(cbind(day,xa)) 
 x5d <- as.numeric(movingAverage(xm[,4],5,centered=T)) ## do 5d mean
 xc <- aggregate(!is.na(xm[,4]),list(xm[,1]),sum,na.rm=T) ## count available observations
 x5m <- aggregate(x5d,list(xm[,1]),max,na.rm=T) ## Rx5d from 5d mean
 x5m[xc[,2]<328,2] <- NA ## kick out years with fewer than 328 days w. measurements per half-year
 ap[,i+1] <- x5m[,2] ### fill matrix of 5d-precipitation output over all stations
 xout <- cbind(1780:2021,NA) # initialize output 
 xmf <- ap[c(1780:2021)>1899&!is.na(ap[,i+1]),i+1] ## gamma-to-normal done only based on data from 1900
 g <- fitdistr(xmf, "gamma")
 xout[!is.na(ap[,i+1]),2] <- qnorm(pgamma(ap[!is.na(ap[,i+1]),i+1],g$estimate[1],g$estimate[2]))
 write.table(xout,paste(outdir,str_replace(lst[i,2], "dailyprecip.txt","Rx5d_gamma_allyear.txt"),sep=""),sep=" ",row.names=F,col.names=F)
 apg[,i+1] <- xout[,2] ### fill matrix of distribution-converted precipitation output over all stations
## season
 if (lst[i,3]=="MJJASO"){selseas <- x[,2]>4&x[,2]<11}else{selseas <- !(x[,2]>4&x[,2]<11)}
 x[!selseas,4] <- NA ### remove out-of-season
 xa <- rep(NA,length(day[,1]))
 xa[match(10000*x[,1]+100*x[,2]+x[,3],10000*day[,1]+100*day[,2]+day[,3])]<-x[,4] ### continuous daily series
 xm <- as.matrix(cbind(day,xa)) 
 x5d <- as.numeric(movingAverage(xm[,4],5,centered=T)) ## do 5d mean
 xc <- aggregate(!is.na(xm[,4]),list(xm[,1]),sum,na.rm=T) ## count available observations
 x5m <- aggregate(x5d,list(xm[,1]),max,na.rm=T) ## Rx5d from 5d mean
 x5m[xc[,2]<328,2] <- NA ## kick out years with fewer than 328 days w. measurements per half-year
 sp[,i+1] <- x5m[,2] ### fill matrix of 5d-precipitation output over all stations
 xout <- cbind(1780:2021,NA) # initialize output 
 xmf <- ap[c(1780:2021)>1899&!is.na(ap[,i+1]),i+1] ## gamma-to-normal done only based on data from 1900
 g <- fitdistr(xmf, "gamma")
 xout[!is.na(ap[,i+1]),2] <- qnorm(pgamma(ap[!is.na(ap[,i+1]),i+1],g$estimate[1],g$estimate[2]))
 write.table(xout,paste(outdir,str_replace(lst[i,2], "dailyprecip.txt","Rx5d_gamma_season.txt"),sep=""),sep=" ",row.names=F,col.names=F)
 spg[,i+1] <- xout[,2] ### fill matrix of distribution-converted precipitation output over all stations

 ### process flood series
 y <- read.table(paste(indir,lst[i,1],sep=""),header=T)
# yearround
 yout <- cbind(1780:2021,NA)
 yout[match(y[,1],1780:2021),2] <- y[,2]
 write.table(yout,str_replace(lst[i,1], ".txt","_flood_yearround_series.txt"),sep=" ",row.names=F,col.names=F)
 af[,i+1] <- yout[,2]
 yg <- yout
 ymf <- yout[c(1780:2021)>1899&!is.na(yout[,2]),2] ## gamma-to-normal done only based on data from 1900
 g <- fitdistr(ymf, "gamma")
 yg[!is.na(yout[,2]),2] <- qnorm(pgamma(yout[!is.na(yout[,2]),2],g$estimate[1],g$estimate[2]))
 write.table(yg,paste(outdir,str_replace(lst[i,1], ".txt","_flood_yearround_gamma.txt"),sep=""),sep=" ",row.names=F,col.names=F)
 afg[,i+1] <- yg[,2]
#season
 if (lst[i,3]=="MJJASO"){selseas <- y[,4]>4&y[,4]<11}else{selseas <- !(y[,4]>4&y[,4]<11)}
 y[!selseas,2] <- NA
 yout <- cbind(1780:2021,NA)
 yout[match(y[,1],1780:2021),2] <- y[,2]
 write.table(yout,str_replace(lst[i,1], ".txt","_flood_season_series.txt"),sep=" ",row.names=F,col.names=F)
 sf[,i+1] <- yout[,2]
 yg <- yout
 ymf <- yout[c(1780:2021)>1899&!is.na(yout[,2]),2] ## gamma-to-normal done only based on data from 1900
 g <- fitdistr(ymf, "gamma")
 yg[!is.na(yout[,2]),2] <- qnorm(pgamma(yout[!is.na(yout[,2]),2],g$estimate[1],g$estimate[2]))
 write.table(yg,paste(outdir,str_replace(lst[i,1], ".txt","_flood_season_gamma.txt"),sep=""),sep=" ",row.names=F,col.names=F)
 sfg[,i+1] <- yg[,2]

 ## smoothing
 aps[,i+1] <- movingAverage(ap[,i+1],30,20,centered=T)
 apgs[,i+1] <- movingAverage(apg[,i+1],30,20,centered=T)
 sps[,i+1] <- movingAverage(sp[,i+1],30,20,centered=T)
 spgs[,i+1] <- movingAverage(spg[,i+1],30,20,centered=T)
 afs[,i+1] <- movingAverage(af[,i+1],30,20,centered=T)
 afgs[,i+1] <- movingAverage(afg[,i+1],30,20,centered=T)
 sfs[,i+1] <- movingAverage(sf[,i+1],30,20,centered=T)
 sfgs[,i+1] <- movingAverage(sfg[,i+1],30,20,centered=T)
 }

# write files 
write.table(ap,paste(outdir,"all_sites_rx5d.yearround.txt",sep=""),sep=" ",row.names=F,col.names=c("yr",lst[,1]))
write.table(apg,paste(outdir,"all_sites_rx5d.yearround.gamma.txt",sep=""),sep=" ",row.names=F,col.names=c("yr",lst[,1]))
write.table(af,paste(outdir,"all_sites_floods.yearround.txt",sep=""),sep=" ",row.names=F,col.names=c("yr",lst[,1]))
write.table(afg,paste(outdir,"all_sites_floods.yearround.gamma.txt",sep=""),sep=" ",row.names=F,col.names=c("yr",lst[,1]))
write.table(sp,paste(outdir,"all_sites_rx5d.season.txt",sep=""),sep=" ",row.names=F,col.names=c("yr",lst[,1]))
write.table(spg,paste(outdir,"all_sites_rx5d.season.gamma.txt",sep=""),sep=" ",row.names=F,col.names=c("yr",lst[,1]))
write.table(sf,paste(outdir,"all_sites_floods.season.txt",sep=""),sep=" ",row.names=F,col.names=c("yr",lst[,1]))
write.table(sfg,paste(outdir,"all_sites_floods.season.gamma.txt",sep=""),sep=" ",row.names=F,col.names=c("yr",lst[,1]))

cor.indiv <- matrix(NA,nrow=dim(apg)[2]-1,ncol=4)
for (i in 2:dim(apg)[2]){
cor.indiv[i-1,1] <- cor(apg[,i],afg[,i],use="complete.obs")
cor.indiv[i-1,2] <- cor.test(apg[,i],afg[,i])$p.value
cor.indiv[i-1,3] <- cor(spg[,i],sfg[,i],use="complete.obs")
cor.indiv[i-1,4] <- cor.test(spg[,i],sfg[,i])$p.value
}
write.table(cor.indiv,paste(outdir,"cor.Rx5d.streamflow.per.river.txt",sep=""),sep=" ",col.names=c("cor.annual","p-value.annual","cor.seasonal","p-value.seasonal"),row.names=lst[,1])

#### Analyses and figures
### Fig 1a, plot longest series
yrs <- 1800:2020
lst <- read.table(paste(indir,"filelist_long.txt",sep=""),header=F)
ynorm <- matrix(NA,nrow=length(yrs),ncol=length(lst[,1])+1) 
ynorm[,1] <- yrs
ynormmv <- ynorm
plotoffset <- c(0,0,1,1,1,2,2,3,3,4,5,6,7,7) # only for plotting
for (i in 1:length(lst[,1])){
    y <- read.table(paste(indir,lst[i,1],sep=""),header=T)
    yscale <- y[,2]/mean(y[y[,1]>1899,2],na.rm=T)
    ymatch <- rep(NA,length(yrs))
    ymatch[match(y[,1],yrs)]<- yscale
    if (i==1){plot(yrs,ylim=c(0,17),ymatch,type="l")}else{lines(yrs,ymatch+(plotoffset[i]*2),col=i)} # plot
    g <- fitdistr(ymatch[yrs>1899&!is.na(ymatch)], "gamma")
    ynorm[,i+1] <- qnorm(pgamma(ymatch,g$estimate[1],g$estimate[2]))
}

### Fig 1b, plot moving averages 
cold <- c(1:10,13,14) ## cold season rivers
cold.ave <- apply(ynorm[,cold+1],1,mean,na.rm=T)
cold.n <- apply(!is.na(ynorm[,cold+1]),1,sum,na.rm=T)
cold.ave[cold.n<length(cold)/2] <- NA
cold.mv <- movingAverage(cold.ave,30,20,centered=T)
plot(yrs,cold.mv,type="l",col="blue")
warm <- c(11:12) ## warm season rivers
warm.ave <- apply(ynorm[,warm+1],1,mean,na.rm=T)
warm.n <- apply(!is.na(ynorm[,warm+1]),1,sum,na.rm=T)
warm.ave[warm.n<=length(warm)/2] <- NA
warm.mv <- movingAverage(warm.ave,30,20,centered=T)
lines(yrs,warm.mv,col="orange")

### Fig 2, plot alpine catchments 
#basel/domatems
plot(afgs[,1],afgs[,28],type="l",ylim=c(-1,1),col="blue")
lines(afgs[,1],afgs[,30],col="lightblue")
#achleiten/martinsbruck
plot(afgs[,1],afgs[,5],type="l",ylim=c(-1,1),col="blue")
lines(afgs[,1],afgs[,19],col="lightblue")
#chancy/porteduscex
plot(afgs[,1],afgs[,37],type="l",ylim=c(-1,1),col="blue")
lines(afgs[,1],afgs[,38],col="lightblue")

### Fig 3, all clusters 
#thames/lea
plot(afgs[,1],afgs[,21],type="l",ylim=c(-1,1.5),col="grey")
lines(afgs[,1],afgs[,41],col="black")
#S-Norway
plot(afgs[,1],afgs[,4],type="l",ylim=c(-1,1.5),col="black")
lines(afgs[,1],afgs[,26],col="blue")
lines(afgs[,1],afgs[,43],col="orange")
#Danube
plot(afgs[,1],afgs[,5],type="l",ylim=c(-1,1.5),col="red")
lines(afgs[,1],afgs[,6],col="purple")
lines(afgs[,1],afgs[,7],col="black")
lines(afgs[,1],afgs[,17],col="purple")
lines(afgs[,1],afgs[,20],col="blue")
lines(afgs[,1],afgs[,22],col="yellow")
lines(afgs[,1],afgs[,39],col="orange")
lines(afgs[,1],afgs[,40],col="green")
#UpperRhine
plot(afgs[,1],afgs[,2],type="l",ylim=c(-1,1.5),col="green")
lines(afgs[,1],afgs[,27],col="black")
lines(afgs[,1],afgs[,28],col="purple")
lines(afgs[,1],afgs[,42],col="red")
#Rhone
plot(afgs[,1],afgs[,36],type="l",ylim=c(-1,1.5),col="black")
lines(afgs[,1],afgs[,37],col="orange")
#elbe
plot(afgs[,1],afgs[,8],type="l",ylim=c(-1,1.5),col="black")
lines(afgs[,1],afgs[,9],col="green")
lines(afgs[,1],afgs[,10],col="orange")
lines(afgs[,1],afgs[,11],col="brown")
lines(afgs[,1],afgs[,12],col="red")
lines(afgs[,1],afgs[,13],col="blue")
#central germany
plot(afgs[,1],afgs[,3],type="l",ylim=c(-1,1.5),col="lightgreen")
lines(afgs[,1],afgs[,18],col="orange")
lines(afgs[,1],afgs[,23],col="blue")
lines(afgs[,1],afgs[,24],col="red")
lines(afgs[,1],afgs[,44],col="black")
lines(afgs[,1],afgs[,45],col="purple")
lines(afgs[,1],afgs[,46],col="brown")
lines(afgs[,1],afgs[,47],col="green")
#lower rhine
plot(afgs[,1],afgs[,14],type="l",ylim=c(-1,1),col="red")
lines(afgs[,1],afgs[,25],col="black")
lines(afgs[,1],afgs[,29],col="purple")
lines(afgs[,1],afgs[,31],col="green")
lines(afgs[,1],afgs[,32],col="orange")
lines(afgs[,1],afgs[,34],col="brown")
#excluded series (supplementary figures)
plot(afgs[,1],afgs[,15],type="l",ylim=c(-1,1),col="lightblue")
lines(afgs[,1],afgs[,30],col="red")
lines(afgs[,1],afgs[,33],col="orange")
lines(afgs[,1],afgs[,35],col="yellow")
lines(afgs[,1],afgs[,19],col="green")
lines(afgs[,1],afgs[,38],col="purple")

## Figure 4 and analyses
### form regional averages for peak streamflow and Rx5d
regap <- matrix(NA,nrow=length(afg[,1]),ncol=9)
regaf <- matrix(NA,nrow=length(afg[,1]),ncol=9)
regaps <- matrix(NA,nrow=length(afg[,1]),ncol=9)
regafs <- matrix(NA,nrow=length(afg[,1]),ncol=9)
regsp <- matrix(NA,nrow=length(afg[,1]),ncol=9)
regsf <- matrix(NA,nrow=length(afg[,1]),ncol=9)
regsps <- matrix(NA,nrow=length(afg[,1]),ncol=9)
regsfs <- matrix(NA,nrow=length(afg[,1]),ncol=9)
regap[,1] <- 1780:2021
regaf[,1] <- 1780:2021
regaps[,1] <- 1780:2021
regafs[,1] <- 1780:2021
regsp[,1] <- 1780:2021
regsf[,1] <- 1780:2021
regsps[,1] <- 1780:2021
regsfs[,1] <- 1780:2021
sel.columns <- list(c(21,41),c(4,26,43),c(5,6,7,17,20,22,39,40),c(2,27,28,42),
                     c(36,37),8:13,c(3,18,23,24,44:47),c(14,25,29,31,32,34))
for (i in 1:8){
 sel <- sel.columns[[i]]
 regaf[,i+1] <- apply(afg[,sel],1,mean,na.rm=T)
 regaf[apply(!is.na(afg[,sel]),1,sum,na.rm=T)<1,i+1] <- NA
 regap[,i+1] <- apply(apg[,sel],1,mean,na.rm=T)
 regap[apply(!is.na(apg[,sel]),1,sum,na.rm=T)<1,i+1] <- NA
 regafs[,i+1] <- movingAverage(regaf[,i+1],30,20,centered=T)
 regaps[,i+1] <- movingAverage(regap[,i+1],30,20,centered=T)
 regsf[,i+1] <- apply(sfg[,sel],1,mean,na.rm=T)
 regsf[apply(!is.na(sfg[,sel]),1,sum,na.rm=T)<1,i+1] <- NA
 regsp[,i+1] <- apply(spg[,sel],1,mean,na.rm=T)
 regsp[apply(!is.na(spg[,sel]),1,sum,na.rm=T)<1,i+1] <- NA
 regsfs[,i+1] <- movingAverage(regsf[,i+1],30,20,centered=T)
 regsps[,i+1] <- movingAverage(regsp[,i+1],30,20,centered=T)
}

# read in Blöschl series
b <- read.table(paste(indir,"Bloeschl_regions.txt",sep=""),header=T)
sel <- c(rep(c(F,F,T,F),60),F,F)
corout <- matrix(NA,nrow=4,ncol=10)
sx <- b[,1]
for (i in 1:8){
 aggregf <- movingAverage(regaf[,i+1],4,centered=T)
 sx <- cbind(sx,aggregf[sel])
 corout[1,i] <- cor(aggregf[sel],b[,i+1],use="complete.obs")
 corout[2,i] <- cor.test(aggregf[sel],b[,i+1])$p.value
 corout[3,i] <- cor.test(aggregf[sel][-1],aggregf[sel][-60])$p.value
 corout[4,i] <- cor.test(b[,i+1][-1],b[,i+1][-60])$p.value
}
corout[1,9] <- cor(apply(sx[,2:3],1,mean),apply(b[,2:3],1,mean),use="complete.obs")
corout[2,9] <- cor.test(apply(sx[,2:3],1,mean),apply(b[,2:3],1,mean))$p.value
corout[3,9] <- cor.test(apply(sx[,2:3],1,mean)[-1],apply(sx[,2:3],1,mean)[-60])$p.value
corout[4,9] <- cor.test(apply(b[,2:3],1,mean)[-1],apply(b[,2:3],1,mean)[-60])$p.value
corout[1,10] <- cor(apply(sx[,4:9],1,mean),apply(b[,4:9],1,mean),use="complete.obs")
corout[2,10] <- cor.test(apply(sx[,4:9],1,mean),apply(b[,4:9],1,mean))$p.value
corout[3,10] <- cor.test(apply(sx[,4:9],1,mean)[-1],apply(sx[,4:9],1,mean)[-60])$p.value
corout[4,10] <- cor.test(apply(b[,4:9],1,mean)[-1],apply(b[,4:9],1,mean)[-60])$p.value
corout 
### correlation for 4-yr aggregation, for each of the 8 regions plus Northwestern Europe (#1,2) and central Europe (#3 to 8) (columns) this gives correlation (1st row),  p-value (2nd), 
## 1st order autocorr of streamflow (3rd), 1st order autocorr of Blöschl (4th) 


bx <- cbind(b[,1],scale(as.numeric(filter(apply(b[,2:3],1,mean),c(1,2,3,4,5,4,3,2,1)/25,method="convolution",sides=2))))
plot(bx[,1],bx[,2],col="purple",type="l",ylim=c(-2.5,2.5))
lines(regafs[,1],scale(apply(regafs[,2:3],1,mean,na.rm=T)),col="black")

bx <- cbind(b[,1],scale(as.numeric(filter(apply(b[,4:9],1,mean),c(1,2,3,4,5,4,3,2,1)/25,method="convolution",sides=2))))
plot(bx[,1],bx[,2],col="purple",type="l")
lines(regafs[43:242,1],scale(apply(regafs[43:242,4:9],1,mean,na.rm=T)),col="black")


### plot panels in Fig. 4 and suppl.
corout4.all <- matrix(NA,nrow=4,ncol=8)
corout4.seas <- matrix(NA,nrow=4,ncol=8)
corout1.all <- matrix(NA,nrow=2,ncol=8)
corout1.seas <- matrix(NA,nrow=2,ncol=8)

dxdq.s.m <- read.table(paste(indir,"CONV5d.seas.30mv.SCALED.MEAN.txt",sep=""),header=T)
dxdq.s.s <- read.table(paste(indir,"CONV5d.seas.30mv.SCALED.SD.txt",sep=""),header=T)
dxdq.a.m <- read.table(paste(indir,"CONV5d.ann.30mv.SCALED.MEAN.txt",sep=""),header=T)
dxdq.a.s <- read.table(paste(indir,"CONV5d.ann.30mv.SCALED.SD.txt",sep=""),header=T)

sel <- c(rep(c(F,F,T,F),60),F,F)
for (i in 1:8){
### plot CONV5d/stramflow/Rx5d in Fig. 4
 dt <- dxdq.s.m[,1]
 dm <- dxdq.s.m[,i+1]
 ds <- dxdq.s.s[,i+1]
 dmin <- -dm+ds
 dmax <- -dm-ds
# postscript(paste(outdir,"CONV5d.seas.region",i,".ps",sep=""),width=4,height=4.8,horizontal=F)
 plot(dt,-dm,type="l",ylim=c(-3,3))
 polygon(c(dt[6:204], rev(dt[6:204])), c(dmax[6:204], rev(dmin[6:204])),col = "grey90", border = NA)
 lines(dt,-dm,col="grey50")
 lines(regsps[,1],scale(regsps[,i+1]),col="blue")
 lines(regsfs[,1],scale(regsfs[,i+1]),col="black")
# dev.off() 
### plot supplement (year round, plus Blöschl)
 dt <- dxdq.a.m[,1]
 dm <- dxdq.a.m[,i+1]
 dss <- dxdq.a.s[,i+1]
 dmin <- -dm-dss
 dmax <- -dm+dss
# postscript(paste(outdir,"CONV5d.all.region",i,".ps",sep=""),width=4,height=4.8,horizontal=F)
 plot(dt,-dm,type="l",ylim=c(-3,3))
 polygon(c(dt[6:204], rev(dt[6:204])), c(dmax[6:204], rev(dmin[6:204])),col = "grey90", border = NA)
 lines(dt,-dm,col="grey50")
 lines(regsps[,1],scale(regsps[,i+1]),col="blue")
 lines(regsfs[,1],scale(regsfs[,i+1]),col="black")
 # dev.off() 
 #### corrleation peaks streamflow / rx5d  
 aggregf <- movingAverage(regaf[,i+1],4,centered=T)
 aggregp <- movingAverage(regap[,i+1],4,centered=T)
 corout1.all[1,i] <- cor(regaf[,i+1],regap[,i+1],use="complete.obs")
 corout1.all[2,i] <- cor.test(regaf[,i+1],regap[,i+1])$p.value
 corout1.seas[1,i] <- cor(regsf[,i+1],regsp[,i+1],use="complete.obs")
 corout1.seas[2,i] <- cor.test(regsf[,i+1],regsp[,i+1])$p.value
 aggregf <- movingAverage(regaf[,i+1],4,centered=T)
 aggregp <- movingAverage(regap[,i+1],4,centered=T)
 corout4.all[1,i] <- cor(aggregf,aggregp,use="complete.obs")
 corout4.all[2,i] <- cor.test(aggregf,aggregp)$p.value
 corout4.all[3,i] <- cor.test(aggregf[-1],aggregf[-60])$p.value
 corout4.all[4,i] <- cor.test(aggregp[-1],aggregp[-60])$p.value
 aggregf <- movingAverage(regsf[,i+1],4,centered=T)
 aggregp <- movingAverage(regsp[,i+1],4,centered=T)
 corout4.seas[1,i] <- cor(aggregf,aggregp,use="complete.obs")
 corout4.seas[2,i] <- cor.test(aggregf,aggregp)$p.value
 corout4.seas[3,i] <- cor.test(aggregf[-1],aggregf[-60])$p.value
 corout4.seas[4,i] <- cor.test(aggregp[-1],aggregp[-60])$p.value
}
corout1.all
corout1.seas
corout4.all
corout4.seas
### for each of the 8 reagions (columns) this gives correlation (1st row),  p-value (2nd), 
## 1st order autocorr of streamflow (3rd), 1st order autocorr of Rx5d (4th) 

## Fig. 5 and analyses
## mean of regional means, seasonal
conv5d <- read.table(paste(indir,"CONV5d.seas.all.txt",sep=""),header=T)[,2]
conv5ds <- read.table(paste(indir,"CONV5d.seas.all.30mv.txt",sep=""),header=T)[,2]
Rx5d <- apply(regsps[,2:9],1,mean,na.rm=T)
streamflow <- apply(regsfs[,2:9],1,mean,na.rm=T)
Rx5d.sd <- apply(regsps[,2:9],1,sd,na.rm=T)
Rx5d.n <- apply(!is.na(regsps[,2:9]),1,sum)
Rx5d.sd[Rx5d.n<4] <- NA
Rx5d[Rx5d.n<4] <- NA
streamflow.sd <- apply(regsfs[,2:9],1,sd,na.rm=T)
streamflow.n <- apply(!is.na(regsfs[,2:9]),1,sum)
streamflow.sd[streamflow.n<4] <- NA
streamflow[streamflow.n<4] <- NA
### plot panel 5a
plot(1806:2015,(-conv5ds*25000)-0.6,type="l",ylim=c(-0.5,1))
lines(regsps[,1],Rx5d,type="l",ylim=c(-0.8,0.8),col="lightblue")
lines(regsps[,1],Rx5d-(Rx5d.sd/sqrt(Rx5d.n)),col="lightblue")
lines(regsps[,1],Rx5d+(Rx5d.sd/sqrt(Rx5d.n)),col="lightblue")
lines(regsfs[,1],streamflow,col="blue")
lines(regsfs[,1],streamflow-(streamflow.sd/sqrt(streamflow.n)),col="blue")
lines(regsfs[,1],streamflow+(streamflow.sd/sqrt(streamflow.n)),col="blue")
### plot wavelets (supplement)
rx <- cbind(1806:2015,scale(Rx5d[27:236]))
st <- cbind(1806:2015,scale(streamflow[27:236]))
dx <- cbind(1806:2015,scale(conv5d))
xwt.prec.flood = xwt(rx[!is.na(rx[,2])&!is.na(st[,2]),], st[!is.na(rx[,2])&!is.na(st[,2]),]) 
xwt.prec.conv = xwt(rx[!is.na(rx[,2])&!is.na(dx[,2]),], dx[!is.na(rx[,2])&!is.na(dx[,2]),]) 
xwt.conv.flood = xwt(dx[!is.na(st[,2])&!is.na(dx[,2]),], st[!is.na(st[,2])&!is.na(dx[,2]),]) 
plot(xwt.prec.flood, plot.cb=TRUE, main = "Cross-Wavelet: Rx5Day vs. Streamflow")
plot(xwt.prec.conv, plot.cb=TRUE, main = "Cross-Wavelet: Rx5Day vs. CONV5d")
plot(xwt.conv.flood, plot.cb=TRUE, main = "Cross-Wavelet: CONV5d vs. Streamflow")

## mean of regional means, annual
conv5d <- read.table(paste(indir,"CONV5d.ann.all.txt",sep=""),header=T)[,2]
conv5ds <- read.table(paste(indir,"CONV5d.ann.all.30mv.txt",sep=""),header=T)[,2]
Rx5d <- apply(regaps[,2:9],1,mean,na.rm=T)
streamflow <- apply(regafs[,2:9],1,mean,na.rm=T)
Rx5d.sd <- apply(regaps[,2:9],1,sd,na.rm=T)
Rx5d.n <- apply(!is.na(regaps[,2:9]),1,sum)
Rx5d.sd[Rx5d.n<4] <- NA
Rx5d[Rx5d.n<4] <- NA
streamflow.sd <- apply(regafs[,2:9],1,sd,na.rm=T)
streamflow.n <- apply(!is.na(regafs[,2:9]),1,sum)
streamflow.sd[streamflow.n<4] <- NA
streamflow[streamflow.n<4] <- NA
### plot panel 5a
plot(1806:2015,(-conv5ds*25000)-0.6,type="l",ylim=c(-1.5,1))
lines(regsps[,1],Rx5d,type="l",ylim=c(-0.8,0.8),col="lightblue")
lines(regsps[,1],Rx5d-(Rx5d.sd/sqrt(Rx5d.n)),col="lightblue")
lines(regsps[,1],Rx5d+(Rx5d.sd/sqrt(Rx5d.n)),col="lightblue")
lines(regsfs[,1],streamflow,col="blue")
lines(regsfs[,1],streamflow-(streamflow.sd/sqrt(streamflow.n)),col="blue")
lines(regsfs[,1],streamflow+(streamflow.sd/sqrt(streamflow.n)),col="blue")
### plot wavelets (supplement)
rx <- cbind(1806:2015,scale(Rx5d[27:236]))
st <- cbind(1806:2015,scale(streamflow[27:236]))
dx <- cbind(1806:2015,scale(conv5d))
xwt.prec.flood = xwt(rx[!is.na(rx[,2])&!is.na(st[,2]),], st[!is.na(rx[,2])&!is.na(st[,2]),]) 
xwt.prec.conv = xwt(rx[!is.na(rx[,2])&!is.na(dx[,2]),], dx[!is.na(rx[,2])&!is.na(dx[,2]),]) 
xwt.conv.flood = xwt(dx[!is.na(st[,2])&!is.na(dx[,2]),], st[!is.na(st[,2])&!is.na(dx[,2]),]) 
plot(xwt.prec.flood, plot.cb=TRUE, main = "Cross-Wavelet: Rx5Day vs. Streamflow")
plot(xwt.prec.conv, plot.cb=TRUE, main = "Cross-Wavelet: Rx5Day vs. CONV5d")
plot(xwt.conv.flood, plot.cb=TRUE, main = "Cross-Wavelet: CONV5d vs. Streamflow")

# plot panel 5b
fig5b <- read.table(paste(indir,"CONV5d.seas.all.partitioned.30mv.txt",sep=""),header=T)
plot(fig5b[,1],-fig5b[,2],ylim=c(-0.00003,0.00006),type="l")
lines(fig5b[,1],-fig5b[,2]-fig5b[,3],col="grey")
lines(fig5b[,1],-fig5b[,2]+fig5b[,3],col="grey")
lines(fig5b[,1],-fig5b[,4],col="green")
lines(fig5b[,1],-fig5b[,4]+fig5b[,5],col="lightgreen")
lines(fig5b[,1],-fig5b[,4]-fig5b[,5],col="lightgreen")
lines(fig5b[,1],-fig5b[,6],col="blue")
lines(fig5b[,1],-fig5b[,6]+fig5b[,7],col="lightblue")
lines(fig5b[,1],-fig5b[,6]-fig5b[,7],col="lightblue")
lines(fig5b[,1],-fig5b[,8],col="purple")
lines(fig5b[,1],-fig5b[,8]+fig5b[,9],col="pink")
lines(fig5b[,1],-fig5b[,8]-fig5b[,9],col="pink")

# Trends
conv5tr <- read.table(paste(indir,"CONV5d.seas.all.partitioned.txt",sep=""),header=T)[158:210,]
tr.m <- lm(conv5tr[,2]~conv5tr[,1])
tr.c <- lm(conv5tr[,4]~conv5tr[,1])
tr.v <- lm(conv5tr[,6]~conv5tr[,1])
rbind(c(summary.lm(tr.m)$coefficients[[2]],summary.lm(tr.m)$coefficients[[8]]),
c(summary.lm(tr.c)$coefficients[[2]],summary.lm(tr.c)$coefficients[[8]]),
c(summary.lm(tr.v)$coefficients[[2]],summary.lm(tr.v)$coefficients[[8]]))

# Fig. S2
ff <- read.table(paste(indir,"flood_frequency.txt",sep=""),header=T)
plot(ff[,1],ff[,5]-1.25,ylim=c(-0.5,0.5),type="l")
lines(ff[,1],ff[,4]-1.25,col="blue")
lines(afgs[,1],sfgs[,36],type="l",ylim=c(-1,1),col="lightblue")
lines(afgs[,1],sfgs[,28],type="l",ylim=c(-1,1),col="grey")


