RtsEva package demo

Author

Alois Tilloy

Demo of the Transformed-Stationary Extreme Value Analysis.

Step 1: load an environmental time series

Here we use with discharge data from the HERA hydrological reanalysis.

River network (rivers with an upstream area > 100km2) for which discharge data is available in HERA

RtsEva contains discharge data from the following four rivers:

Ardèche @ Saint Martin:

Total length: 120 km
Mean flow: 50 m³/s
Short description: The Ardèche is mediterranean river mostly known for tourism due to its scenic gorges. The river is known for its extreme autumn floods during the so called “épisodes cévenols”.
More: Wiki

Danube @ Vienna:

Total length: 2850 km
Mean flow: 1800 m³/s
Short description: The Danube is the longest river in the European Union. It is an important waterway for international trade, connecting several countries in Central and Eastern Europe.
More: Wiki

Ebro @ Zaragoza:

Total length: 930 km
Mean flow: 197 m³/s
Short description: The Ebro is Spain’s longest river, with low and high water levels alternating throughout the year, influenced by winter snowmelt and summer evaporation/human usage. The river is vital for agriculture.
More: Wiki

Rhône @ Lyon:

Total length: 813 km
Mean flow: 507 m³/s
Short description: The Rhône is France’s most powerful river, characterized by a significant seasonal variation in flow rates. The Rhône is crucial for transportation, hydropower generation, and irrigation in the region.
More: Wiki

We choose the Ardèche for this example. It is however possible to test different rivers in this Shiny app

Code

data=Ardeche$data
longitude=Ardeche$river$lon
latitude=Ardeche$river$lat
timeStamps=data$time

Location of the selected river:

Plot of the original timeseries

Hydrograph of river Ardeche

Different preprocessing depending on which hazard is studied:

Flood:
- no transformation
- minimum distance between two peaks: 10 days
- analysis performed on daily maximum discharge
Drought:
- reversed transformation
- minimum distance between two peaks: 30 days
- analysis performed on 7-days moving averaged discharged

Step 2: Run TSEVA for flood

Code

haz = "flood"

dt1=min(diff(data$time),na.rm=T)
dt=as.numeric(dt1)
tdim=attributes(dt1)$units
if (tdim=="hours") dt=dt/24
if (dt==1){
  timeDays=data$time
}else{
  timeDays=unique(as.Date(data$time))
}
if (haz=="drought"){
  
  minPeakDistanceInDays=30
  tail="low"
  #7 day average flow for drought analysis
  WindowSize=7
  names(data)=c("date","Qs")
  if (tdim=="hours") dt=dt/24
  nRunMn = ceiling(WindowSize/dt);
  colnames(data)
  data$Q7=rollmean(data$Qs,nRunMn, align = "right", fill=NA)
  
  if (!exists("trans")){trans="rev"}
  print(paste0(trans," transformation used for low flows"))
  
  series=data.frame(data$date,data$Q7)
  
}else if (haz=="flood"){
  tail="high"
  percentile=95
  minPeakDistanceInDays=10
  series <- max_daily_value(data)
  

}
timeAndSeries=series
names(timeAndSeries)=c("timestamp","data")

tsm=1/dt
series=timeAndSeries[,2]
timeWindow = 365.25*30; #time windows in days, the correction is done within the functions
windowSize=366

timeStamps=timeAndSeries$timestamp
transftypes=c("ori","rev")
trendtypes=c("trend","trendPeaks")
trendtrans=expand.grid(transftypes,trendtypes)
trendtrans = data.frame(trendtrans)
names(trendtrans)=c("transformation","trend method")

tt=3
print(trendtrans[tt,])

  transformation trend method
3            ori   trendPeaks

Code

Nonstat<-TsEvaNs(timeAndSeries, timeWindow, transfType=trendtrans[tt,2], ciPercentile= 90, minPeakDistanceInDays = minPeakDistanceInDays, tail, tail = tail, trans=trendtrans[tt,1])
nonStationaryEvaParams=Nonstat[[1]]
stationaryTransformData=Nonstat[[2]]

datex=yday(timeDays)
dtect=c(diff(datex),-1)
last_days <- timeDays[which(dtect<0)]
tindexes=match(last_days,timeDays)
RPgoal=10
timeIndex=1
RLevs100a=ComputeReturnLevels(nonStationaryEvaParams, RPgoal, timeIndex,trans=trendtrans[tt,1])
timeIndex=tindexes[length(tindexes)]-1
RLevs100b=ComputeReturnLevels(nonStationaryEvaParams, RPgoal, timeIndex,trans=trendtrans[tt,1])

FloodRL=data.frame(c(round(RLevs100a$ReturnLevels[3],2),
                     round(RLevs100a$ReturnLevels[3]-RLevs100a$ReturnLevels[5],2),
                     round(RLevs100a$ReturnLevels[3]+RLevs100a$ReturnLevels[5],2)),
                   c(round(RLevs100b$ReturnLevels[3],2),
                     round(RLevs100b$ReturnLevels[3]-RLevs100b$ReturnLevels[5],2),
                     round(RLevs100b$ReturnLevels[3]+RLevs100b$ReturnLevels[5],2)))
colnames(FloodRL)=c("Y1951","Y2020")
rownames(FloodRL)=c("RL","infCI","SupCI")
FloodRL=as.data.frame(FloodRL)
FloodRL$change=round((FloodRL$Y2020-FloodRL$Y1951)/FloodRL$Y1951*100,1)

for more details on the different transformation and trend methods, the reader can refer to the RtsEva package documentation: https://cran.r-project.org/web/packages/RtsEva/RtsEva.pdf

Step 4: Plot results for flood

Timeseries with peaks and evolving GPD density

Return Level curve in 1951

Curve beam (all curves from 1951 to 2020)

Step 5: Run TSEVA for drought

Code

haz = "drought"
rm(trans)

dt1=min(diff(data$time),na.rm=T)
dt=as.numeric(dt1)
tdim=attributes(dt1)$units
if (tdim=="hours") dt=dt/24
if (dt==1){
  timeDays=data$time
}else{
  timeDays=unique(as.Date(data$time))
}
if (haz=="drought"){
  
  minPeakDistanceInDays=30
  tail="low"
  #7 day average flow for drought analysis
  WindowSize=7
  names(data)=c("date","Qs")
  if (tdim=="hours") dt=dt/24
  nRunMn = ceiling(WindowSize/dt);
  colnames(data)
  data$Q7=rollmean(data$Qs,nRunMn, align = "right", fill=NA)
  
  if (!exists("trans")){trans="rev"}
  print(paste0(trans," transformation used for low flows"))
  
  series=data.frame(data$date,data$Q7)
  
}else if (haz=="flood"){
  tail="high"
  percentile=95
  minPeakDistanceInDays=10
  series <- max_daily_value(data)
  

}

[1] "rev transformation used for low flows"

Code

timeAndSeries=series
names(timeAndSeries)=c("timestamp","data")
timeStamps=timeAndSeries$timestamp

tt=4
print(trendtrans[tt,])

  transformation trend method
4            rev   trendPeaks

Code

Nonstat<-TsEvaNs(timeAndSeries, timeWindow, transfType=trendtrans[tt,2], ciPercentile= 90, minPeakDistanceInDays = minPeakDistanceInDays, tail, tail = tail, trans=trendtrans[tt,1])
nonStationaryEvaParams=Nonstat[[1]]
stationaryTransformData=Nonstat[[2]]

stationaryTransformData$timeStampsDay=unique(as.Date(as.character(stationaryTransformData$timeStamps)))
dt1=min(diff(timeStamps),na.rm=T)
dt=as.numeric(dt1)
tdim=attributes(dt1)$units
if (tdim=="hours") dt=dt/24
if (dt==1){
  timeDays=stationaryTransformData$timeStamps
}else{
  timeDays=stationaryTransformData$timeStampsDay
    }
    
datex=yday(timeDays)
dtect=c(diff(datex),-1)
last_days <- timeDays[which(dtect<0)]
tindexes=match(last_days,timeDays)
RPgoal=10
timeIndex=1
RLevs100a=ComputeReturnLevels(nonStationaryEvaParams, RPgoal, timeIndex,trans=trendtrans[tt,1])
timeIndex=tindexes[length(tindexes)]-1
RLevs100b=ComputeReturnLevels(nonStationaryEvaParams, RPgoal, timeIndex,trans=trendtrans[tt,1])

DryRL=data.frame(c(round(RLevs100a$ReturnLevels[3],2),
                     round(RLevs100a$ReturnLevels[3]-RLevs100a$ReturnLevels[5],2),
                     round(RLevs100a$ReturnLevels[3]+RLevs100a$ReturnLevels[5],2)),
                   c(round(RLevs100b$ReturnLevels[3],2),
                     round(RLevs100b$ReturnLevels[3]-RLevs100b$ReturnLevels[5],2),
                     round(RLevs100b$ReturnLevels[3]+RLevs100b$ReturnLevels[5],2)))
colnames(DryRL)=c("Y1951","Y2020")
rownames(DryRL)=c("RL","infCI","SupCI")
DryRL=as.data.frame(DryRL)
DryRL$change=round((DryRL$Y2020-DryRL$Y1951)/DryRL$Y1951*100,1)

Step 6: Plot results for drought

Timeseries with peaks and evolving GPD density

Return Level curve in 1951

Curve bean (all curves from 1951 to 2020)

Summary of the results for the Ardèche

The table shows the 10 years return level of high flows (daily maximum flow) and low flows (7 days averaged flow) for 1951 and 2020.

Changes in high and low flows
RL in 1951 (m3/s) [68% CI]	RL in 2020 (m3/s) [68% CI]	Relative change (%) [68% CI]
831.39 [762.24 -900.53]	799.23 [732.18 - 866.27]	-3.9 [-3.9 - -3.8]
19.05 [17.24 - 20.87]	18.42 [16.64 - 20.2]	-3.3 [-3.5 - -3.2]

Although the river displays a non-linear trend, the transformed-stationary method estimates a slight decrease in both extreme high an low flows on the Ardèche between 1951 and 2020.

--- title: "RtsEva package demo" author: "Alois Tilloy" format: html: html-math-method: katex code-tools: true self-contained: true execute: warning: false editor: visual toc: true editor_options: chunk_output_type: console --- ::: grid ::: {.headline .g-col-lg-9 .g-col-12 .g-col-md-12} ::: h1 Demo of the Transformed-Stationary Extreme Value Analysis. ::: ::: ::: {.g-col-lg-3 .g-col-12 .g-col-md-12} [![](RtsEva.png){fig.align="right" height="200px"}](https://cran.r-project.org/web/packages/RtsEva/index.html) ::: ::: ```{r} #| echo: false #| warning: false #Loading RtsEva package library(RtsEva) library(zoo) library(lubridate) library(knitr) ComputeReturnLevels<-function(nonStationaryEvaParams, RPgoal, timeIndex,trans=NA){ #GEV epsilonGEV <- nonStationaryEvaParams[[1]]$parameters$epsilon sigmaGEV <- mean(nonStationaryEvaParams[[1]]$parameters$sigma[timeIndex]) muGEV <- mean(nonStationaryEvaParams[[1]]$parameters$mu[timeIndex]) dtSampleYears <- nonStationaryEvaParams[[1]]$parameters$timeDeltaYears #GPD epsilonGPD <- nonStationaryEvaParams[[2]]$parameters$epsilon sigmaGPD <- mean(nonStationaryEvaParams[[2]]$parameters$sigma[timeIndex]) thresholdGPD <- mean(nonStationaryEvaParams[[2]]$parameters$threshold[timeIndex]) nPeaks <- nonStationaryEvaParams[[2]]$parameters$nPeaks thStart <- nonStationaryEvaParams[[2]]$parameters$timeHorizonStart thEnd <- nonStationaryEvaParams[[2]]$parameters$timeHorizonEnd sampleTimeHorizon <- as.numeric((thEnd - thStart)/365.2425) if (nonStationaryEvaParams[[1]]$method=="No fit"){ print("could not fit EVD to this pixel") ParamGEV=c(epsilonGEV,sigmaGEV,muGEV,NA, NA, NA) names(ParamGEV)=c("epsilonGEV","sigmaGEV","muGEV","epsilonStdErrGEV","sigmaStdErrGEV","muStdErrGEV") ParamGPD=c(epsilonGPD,sigmaGPD,thresholdGPD,NA,NA, NA,nPeaks,sampleTimeHorizon) names(ParamGPD)=c("epsilonGPD","sigmaGPD","thresholdGPD","epsilonStdErrGPD","sigmaStdErrGPD","thresholdStdErrGPD","nPeaks","SampleTimeHorizon") return(list(Fit="No fit",Params=c(ParamGEV,ParamGPD))) }else{ #GEV # epsilonGEV <- nonStationaryEvaParams[[1]]$parameters$epsilon # sigmaGEV <- mean(nonStationaryEvaParams[[1]]$parameters$sigma[timeIndex]) # muGEV <- mean(nonStationaryEvaParams[[1]]$parameters$mu[timeIndex]) # dtSampleYears <- nonStationaryEvaParams[[1]]$parameters$timeDeltaYears epsilonStdErrGEV <- nonStationaryEvaParams[[1]]$paramErr$epsilonErr sigmaStdErrGEV <- mean(nonStationaryEvaParams[[1]]$paramErr$sigmaErr[timeIndex]) muStdErrGEV <- mean(nonStationaryEvaParams[[1]]$paramErr$muErr[timeIndex]) #GPD # epsilonGPD <- nonStationaryEvaParams[[2]]$parameters$epsilon # sigmaGPD <- mean(nonStationaryEvaParams[[2]]$parameters$sigma[timeIndex]) # thresholdGPD <- mean(nonStationaryEvaParams[[2]]$parameters$threshold[timeIndex]) # nPeaks <- nonStationaryEvaParams[[2]]$parameters$nPeaks epsilonStdErrGPD <- nonStationaryEvaParams[[2]]$paramErr$epsilonErr sigmaStdErrGPD <- mean(nonStationaryEvaParams[[2]]$paramErr$sigmaErr[timeIndex]) thresholdStdErrGPD <- mean(nonStationaryEvaParams[[2]]$paramErr$thresholdErr[timeIndex]) # thStart <- nonStationaryEvaParams[[2]]$parameters$timeHorizonStart # thEnd <- nonStationaryEvaParams[[2]]$parameters$timeHorizonEnd # sampleTimeHorizon <- as.numeric((thEnd - thStart)/365.2425) if (trans=="rev"){ sigmaGEV=-sigmaGEV muGEV=-muGEV sigmaGPD=-sigmaGPD thresholdGPD=-thresholdGPD } returnLevelsGEV <- tsEvaComputeReturnLevelsGEV(epsilonGEV, sigmaGEV, muGEV, epsilonStdErrGEV, sigmaStdErrGEV, muStdErrGEV, RPgoal) returnLevelsGPD <- tsEvaComputeReturnLevelsGPD(epsilonGPD, sigmaGPD, thresholdGPD, epsilonStdErrGPD, sigmaStdErrGPD, thresholdStdErrGPD, nPeaks, sampleTimeHorizon, RPgoal) rlevGEV=returnLevelsGEV$returnLevels rlevGPD=returnLevelsGPD$returnLevels errGEV=returnLevelsGEV$returnLevelsErr errGPD=returnLevelsGPD$returnLevelsErr ParamGEV=c(epsilonGEV,sigmaGEV,muGEV,epsilonStdErrGEV, sigmaStdErrGEV, muStdErrGEV) names(ParamGEV)=c("epsilonGEV","sigmaGEV","muGEV","epsilonStdErrGEV","sigmaStdErrGEV","muStdErrGEV") ParamGPD=c(epsilonGPD,sigmaGPD,thresholdGPD,epsilonStdErrGPD,sigmaStdErrGPD, thresholdStdErrGPD,nPeaks,sampleTimeHorizon) names(ParamGPD)=c("epsilonGPD","sigmaGPD","thresholdGPD","epsilonStdErrGPD","sigmaStdErrGPD","thresholdStdErrGPD","nPeaks","SampleTimeHorizon") return(list(Fit="Fitted",ReturnLevels=c(ReturnPeriod=RPgoal, GEV=as.numeric(rlevGEV),GPD=as.numeric(rlevGPD),errGEV=as.numeric(errGEV),errGPD=as.numeric(errGPD)),Params=c(ParamGEV,ParamGPD))) } } ``` ## Step 1: load an environmental time series ### Here we use with discharge data from the [HERA hydrological reanalysis](https://data.jrc.ec.europa.eu/dataset/a605a675-9444-4017-8b34-d66be5b18c95). ![River network (rivers with an upstream area > 100km2) for which discharge data is available in HERA](HERA_network.png) ### RtsEva contains discharge data from the following four rivers: #### Ardèche \@ Saint Martin: - Total length: 120 km - Mean flow: 50 m^3^/s - Short description: The Ardèche is mediterranean river mostly known for tourism due to its scenic gorges. The river is known for its extreme autumn floods during the so called "épisodes cévenols". - More: [Wiki](https://en.wikipedia.org/wiki/Ard%C3%A8che_(river)) #### Danube \@ Vienna: - Total length: 2850 km - Mean flow: 1800 m^3^/s - Short description: The Danube is the longest river in the European Union. It is an important waterway for international trade, connecting several countries in Central and Eastern Europe. - More: [Wiki](https://en.wikipedia.org/wiki/Danube) #### Ebro \@ Zaragoza: - Total length: 930 km - Mean flow: 197 m^3^/s - Short description: The Ebro is Spain's longest river, with low and high water levels alternating throughout the year, influenced by winter snowmelt and summer evaporation/human usage. The river is vital for agriculture. - More: [Wiki](https://en.wikipedia.org/wiki/Ebro) #### Rhône \@ Lyon: - Total length: 813 km - Mean flow: 507 m^3^/s - Short description: The Rhône is France's most powerful river, characterized by a significant seasonal variation in flow rates. The Rhône is crucial for transportation, hydropower generation, and irrigation in the region. - More: [Wiki](https://en.wikipedia.org/wiki/Rh%C3%B4ne) ```{r} #| echo: false #| warning: false Ebro=list(Nsq=42, river=data.frame("lon"=-0.825,"lat"=41.608), catch="Ebro @ Zaragoza", data=EbroZaragoza) Ardeche=list(Nsq=42, river=data.frame("lon"=4.575,"lat"=44.297), catch="Ardeche @ Saint-Martin", data=ArdecheStMartin) Rhone=list(Nsq=52, river=data.frame("lon"=4.891,"lat"=45.772), catch="Rhone @ Lyon", data=RhoneLyon) Danube=list(Nsq=63, river=data.frame("lon"=16.64,"lat"=48.13), catch="Danube @ Vienna", data=DanubeVienna) ``` #### We choose the Ardèche for this example. It is however possible to test different rivers in this [Shiny app](https://alowis.shinyapps.io/RtsEva_demo/) ```{r} #| code-fold: true data=Ardeche$data longitude=Ardeche$river$lon latitude=Ardeche$river$lat timeStamps=data$time ``` ### Location of the selected river: ```{r} #| label: map #| echo: false library(leaflet) leaflet() |> setView(longitude, latitude, zoom = 12) |> addTiles() |> addMarkers(longitude, latitude, popup = "Maungawhau") ``` ### Plot of the original timeseries ```{r} #| label: hydrograph #| fig-cap: "Hydrograph of river Ardeche" #| echo: false library(dygraphs) dygraph(data) %>% dyRangeSelector(dateWindow = c("1951-01-01", "1960-01-01")) ``` Different preprocessing depending on which hazard is studied: - Flood: - no transformation - minimum distance between two peaks: 10 days - analysis performed on daily maximum discharge - Drought: - reversed transformation - minimum distance between two peaks: 30 days - analysis performed on 7-days moving averaged discharged ## Step 2: Run TSEVA for flood ```{r} #| code-fold: true haz = "flood" dt1=min(diff(data$time),na.rm=T) dt=as.numeric(dt1) tdim=attributes(dt1)$units if (tdim=="hours") dt=dt/24 if (dt==1){ timeDays=data$time }else{ timeDays=unique(as.Date(data$time)) } if (haz=="drought"){ minPeakDistanceInDays=30 tail="low" #7 day average flow for drought analysis WindowSize=7 names(data)=c("date","Qs") if (tdim=="hours") dt=dt/24 nRunMn = ceiling(WindowSize/dt); colnames(data) data$Q7=rollmean(data$Qs,nRunMn, align = "right", fill=NA) if (!exists("trans")){trans="rev"} print(paste0(trans," transformation used for low flows")) series=data.frame(data$date,data$Q7) }else if (haz=="flood"){ tail="high" percentile=95 minPeakDistanceInDays=10 series <- max_daily_value(data) } timeAndSeries=series names(timeAndSeries)=c("timestamp","data") tsm=1/dt series=timeAndSeries[,2] timeWindow = 365.25*30; #time windows in days, the correction is done within the functions windowSize=366 timeStamps=timeAndSeries$timestamp transftypes=c("ori","rev") trendtypes=c("trend","trendPeaks") trendtrans=expand.grid(transftypes,trendtypes) trendtrans = data.frame(trendtrans) names(trendtrans)=c("transformation","trend method") tt=3 print(trendtrans[tt,]) Nonstat<-TsEvaNs(timeAndSeries, timeWindow, transfType=trendtrans[tt,2], ciPercentile= 90, minPeakDistanceInDays = minPeakDistanceInDays, tail, tail = tail, trans=trendtrans[tt,1]) nonStationaryEvaParams=Nonstat[[1]] stationaryTransformData=Nonstat[[2]] datex=yday(timeDays) dtect=c(diff(datex),-1) last_days <- timeDays[which(dtect<0)] tindexes=match(last_days,timeDays) RPgoal=10 timeIndex=1 RLevs100a=ComputeReturnLevels(nonStationaryEvaParams, RPgoal, timeIndex,trans=trendtrans[tt,1]) timeIndex=tindexes[length(tindexes)]-1 RLevs100b=ComputeReturnLevels(nonStationaryEvaParams, RPgoal, timeIndex,trans=trendtrans[tt,1]) FloodRL=data.frame(c(round(RLevs100a$ReturnLevels[3],2), round(RLevs100a$ReturnLevels[3]-RLevs100a$ReturnLevels[5],2), round(RLevs100a$ReturnLevels[3]+RLevs100a$ReturnLevels[5],2)), c(round(RLevs100b$ReturnLevels[3],2), round(RLevs100b$ReturnLevels[3]-RLevs100b$ReturnLevels[5],2), round(RLevs100b$ReturnLevels[3]+RLevs100b$ReturnLevels[5],2))) colnames(FloodRL)=c("Y1951","Y2020") rownames(FloodRL)=c("RL","infCI","SupCI") FloodRL=as.data.frame(FloodRL) FloodRL$change=round((FloodRL$Y2020-FloodRL$Y1951)/FloodRL$Y1951*100,1) ``` for more details on the different transformation and trend methods, the reader can refer to the RtsEva package documentation: <https://cran.r-project.org/web/packages/RtsEva/RtsEva.pdf> ## Step 4: Plot results for flood ### Timeseries with peaks and evolving GPD density ```{r} #| label: plot1 #| warning: false #| fig.width: 14 #| fig.height: 7 #| echo: false trans=trendtrans[tt,1] ExRange= c(min(nonStationaryEvaParams$potObj$parameters$peaks),max(nonStationaryEvaParams$potObj$parameters$peaks)) if (haz=="flood") wr2 <- c(seq(0.8*min(ExRange),1.2*max(ExRange),length.out=300)) if (haz=="drought") wr2 <- c(seq(1.2*min(ExRange),0.2*max(ExRange),length.out=300)) Plot1= tsEvaPlotGPDImageScFromAnalysisObj(wr2, nonStationaryEvaParams, stationaryTransformData,trans=trans) Plot1 ``` ### Return Level curve in 1951 ```{r} #| label: plotRL #| warning: false #| fig.width: 12 #| fig.height: 6 #| echo: false timeIndex=1 epsilon <- nonStationaryEvaParams[[2]]$parameters$epsilon sigma <- mean(nonStationaryEvaParams[[2]]$parameters$sigma[timeIndex]) threshold <- mean(nonStationaryEvaParams[[2]]$parameters$threshold[timeIndex]) thStart <- nonStationaryEvaParams[[2]]$parameters$timeHorizonStart thEnd <- nonStationaryEvaParams[[2]]$parameters$timeHorizonEnd timeHorizonInYears <- round(as.numeric((thEnd - thStart) / 365.2425)) nPeaks <- nonStationaryEvaParams[[2]]$parameters$nPeaks peax <- nonStationaryEvaParams[[2]]$parameters$peaks peaxID <- nonStationaryEvaParams[[2]]$parameters$peakID timeStamps <- stationaryTransformData$timeStamps dt1 <- min(diff(timeStamps), na.rm = T) dt <- as.numeric(dt1) tdim <- attributes(dt1)$units if (tdim == "hours") dt <- dt / 24 if (dt >= 1) { timeStamps <- stationaryTransformData$timeStamps trendPeaks <- stationaryTransformData$trendSeries[peaxID] stdPeaks <- stationaryTransformData$stdDevSeries[peaxID] } else { timeStamps <- stationaryTransformData$timeStampsDay trendPeaks <- stationaryTransformData$trendSeriesOr[peaxID] stdPeaks <- stationaryTransformData$stdDevSeriesOr[peaxID] dt <- 1 } peaksCor <- (peax - trendPeaks) / stdPeaks tstamps <- timeStamps[timeIndex] epsilonStdErr <- nonStationaryEvaParams[[2]]$paramErr$epsilonErr sigmaStdErr <- mean(nonStationaryEvaParams[[2]]$paramErr$sigmaErr[timeIndex]) thresholdStdErr <- mean(nonStationaryEvaParams[[2]]$paramErr$thresholdErr[timeIndex]) peaxday <- lubridate::yday(stationaryTransformData$timeStamps[peaxID]) xday <- lubridate::yday(stationaryTransformData$timeStamps[timeIndex]) * (2 * pi / 365.25) seasonalityTimeWindow <- 30.4 * 12 maxD <- sqrt(2 - 2 * cos((2 * pi / 365.25) - seasonalityTimeWindow * (pi / 365.25))) seasonRatio <- seasonalityTimeWindow / 365.25 theta <- peaxday * (2 * pi / 365.25) D <- sqrt(2 - 2 * cos(xday - theta)) sel <- which(D < (maxD)) nYears <- round(length(timeStamps) / 365.25 * dt) rlvlpeax <- empdis(peaksCor, nYears) rlvlpeax$QNS <- peax[order(peax)] rlvlpeax$Idt <- stationaryTransformData$timeStamps[peaxID][order(peax)] rlvlpeam <- empdis(peaksCor[sel], nYears * seasonRatio) rlvlpeam$QNS <- peax[sel][order(peax[sel])] rlvlpeam$Idt <- stationaryTransformData$timeStamps[peaxID[sel]][order(peax[sel])] if (nonStationaryEvaParams[[1]]$parameters$timeDeltaYears < 1) { rlvmax <- rlvlpeam } else { rlvmax <- rlvlpeax } varargin=NA RLtstep <- tsEvaPlotReturnLevelsGPD(epsilon, sigma, threshold, epsilonStdErr, sigmaStdErr, thresholdStdErr, nPeaks = nPeaks, timeHorizonInYears = timeHorizonInYears, rlvmax, tstamps, trans,varargin) RLtstep ``` ### Curve beam (all curves from 1951 to 2020) ```{r} #| label: plotBeam #| warning: false #| fig.width: 12 #| fig.height: 6 #| echo: false Plot2 = tsEvaPlotReturnLevelsGPDFromAnalysisObj(nonStationaryEvaParams, stationaryTransformData, timeIndex, trans=trans,ylabel="Discharge (m3/s)") Plot2 ``` ## Step 5: Run TSEVA for drought ```{r} #| label: TSEVAdrought #| code-fold: true haz = "drought" rm(trans) dt1=min(diff(data$time),na.rm=T) dt=as.numeric(dt1) tdim=attributes(dt1)$units if (tdim=="hours") dt=dt/24 if (dt==1){ timeDays=data$time }else{ timeDays=unique(as.Date(data$time)) } if (haz=="drought"){ minPeakDistanceInDays=30 tail="low" #7 day average flow for drought analysis WindowSize=7 names(data)=c("date","Qs") if (tdim=="hours") dt=dt/24 nRunMn = ceiling(WindowSize/dt); colnames(data) data$Q7=rollmean(data$Qs,nRunMn, align = "right", fill=NA) if (!exists("trans")){trans="rev"} print(paste0(trans," transformation used for low flows")) series=data.frame(data$date,data$Q7) }else if (haz=="flood"){ tail="high" percentile=95 minPeakDistanceInDays=10 series <- max_daily_value(data) } timeAndSeries=series names(timeAndSeries)=c("timestamp","data") timeStamps=timeAndSeries$timestamp tt=4 print(trendtrans[tt,]) Nonstat<-TsEvaNs(timeAndSeries, timeWindow, transfType=trendtrans[tt,2], ciPercentile= 90, minPeakDistanceInDays = minPeakDistanceInDays, tail, tail = tail, trans=trendtrans[tt,1]) nonStationaryEvaParams=Nonstat[[1]] stationaryTransformData=Nonstat[[2]] stationaryTransformData$timeStampsDay=unique(as.Date(as.character(stationaryTransformData$timeStamps))) dt1=min(diff(timeStamps),na.rm=T) dt=as.numeric(dt1) tdim=attributes(dt1)$units if (tdim=="hours") dt=dt/24 if (dt==1){ timeDays=stationaryTransformData$timeStamps }else{ timeDays=stationaryTransformData$timeStampsDay } datex=yday(timeDays) dtect=c(diff(datex),-1) last_days <- timeDays[which(dtect<0)] tindexes=match(last_days,timeDays) RPgoal=10 timeIndex=1 RLevs100a=ComputeReturnLevels(nonStationaryEvaParams, RPgoal, timeIndex,trans=trendtrans[tt,1]) timeIndex=tindexes[length(tindexes)]-1 RLevs100b=ComputeReturnLevels(nonStationaryEvaParams, RPgoal, timeIndex,trans=trendtrans[tt,1]) DryRL=data.frame(c(round(RLevs100a$ReturnLevels[3],2), round(RLevs100a$ReturnLevels[3]-RLevs100a$ReturnLevels[5],2), round(RLevs100a$ReturnLevels[3]+RLevs100a$ReturnLevels[5],2)), c(round(RLevs100b$ReturnLevels[3],2), round(RLevs100b$ReturnLevels[3]-RLevs100b$ReturnLevels[5],2), round(RLevs100b$ReturnLevels[3]+RLevs100b$ReturnLevels[5],2))) colnames(DryRL)=c("Y1951","Y2020") rownames(DryRL)=c("RL","infCI","SupCI") DryRL=as.data.frame(DryRL) DryRL$change=round((DryRL$Y2020-DryRL$Y1951)/DryRL$Y1951*100,1) ``` ## Step 6: Plot results for drought ### Timeseries with peaks and evolving GPD density ```{r} #| label: plot1D #| warning: false #| fig.width: 14 #| fig.height: 7 #| echo: false trans=trendtrans[tt,1] ExRange= c(min(nonStationaryEvaParams$potObj$parameters$peaks),max(nonStationaryEvaParams$potObj$parameters$peaks)) if (haz=="flood") wr2 <- c(seq(0.8*min(ExRange),1.2*max(ExRange),length.out=300)) if (haz=="drought") wr2 <- c(seq(1.2*min(ExRange),0.2*max(ExRange),length.out=300)) Plot1= tsEvaPlotGPDImageScFromAnalysisObj(wr2, nonStationaryEvaParams, stationaryTransformData,trans=trans) Plot1 ``` ### Return Level curve in 1951 ```{r} #| label: plotRLD #| warning: false #| fig.width: 12 #| fig.height: 6 #| echo: false timeIndex=1 epsilon <- nonStationaryEvaParams[[2]]$parameters$epsilon sigma <- mean(nonStationaryEvaParams[[2]]$parameters$sigma[timeIndex]) threshold <- mean(nonStationaryEvaParams[[2]]$parameters$threshold[timeIndex]) thStart <- nonStationaryEvaParams[[2]]$parameters$timeHorizonStart thEnd <- nonStationaryEvaParams[[2]]$parameters$timeHorizonEnd timeHorizonInYears <- round(as.numeric((thEnd - thStart) / 365.2425)) nPeaks <- nonStationaryEvaParams[[2]]$parameters$nPeaks peax <- nonStationaryEvaParams[[2]]$parameters$peaks peaxID <- nonStationaryEvaParams[[2]]$parameters$peakID timeStamps <- stationaryTransformData$timeStamps dt1 <- min(diff(timeStamps), na.rm = T) dt <- as.numeric(dt1) tdim <- attributes(dt1)$units if (tdim == "hours") dt <- dt / 24 if (dt >= 1) { timeStamps <- stationaryTransformData$timeStamps trendPeaks <- stationaryTransformData$trendSeries[peaxID] stdPeaks <- stationaryTransformData$stdDevSeries[peaxID] } else { timeStamps <- stationaryTransformData$timeStampsDay trendPeaks <- stationaryTransformData$trendSeriesOr[peaxID] stdPeaks <- stationaryTransformData$stdDevSeriesOr[peaxID] dt <- 1 } peaksCor <- (peax - trendPeaks) / stdPeaks tstamps <- timeStamps[timeIndex] epsilonStdErr <- nonStationaryEvaParams[[2]]$paramErr$epsilonErr sigmaStdErr <- mean(nonStationaryEvaParams[[2]]$paramErr$sigmaErr[timeIndex]) thresholdStdErr <- mean(nonStationaryEvaParams[[2]]$paramErr$thresholdErr[timeIndex]) peaxday <- lubridate::yday(stationaryTransformData$timeStamps[peaxID]) xday <- lubridate::yday(stationaryTransformData$timeStamps[timeIndex]) * (2 * pi / 365.25) seasonalityTimeWindow <- 30.4 * 12 maxD <- sqrt(2 - 2 * cos((2 * pi / 365.25) - seasonalityTimeWindow * (pi / 365.25))) seasonRatio <- seasonalityTimeWindow / 365.25 theta <- peaxday * (2 * pi / 365.25) D <- sqrt(2 - 2 * cos(xday - theta)) sel <- which(D < (maxD)) nYears <- round(length(timeStamps) / 365.25 * dt) rlvlpeax <- empdis(peaksCor, nYears) rlvlpeax$QNS <- peax[order(peax)] rlvlpeax$Idt <- stationaryTransformData$timeStamps[peaxID][order(peax)] rlvlpeam <- empdis(peaksCor[sel], nYears * seasonRatio) rlvlpeam$QNS <- peax[sel][order(peax[sel])] rlvlpeam$Idt <- stationaryTransformData$timeStamps[peaxID[sel]][order(peax[sel])] if (nonStationaryEvaParams[[1]]$parameters$timeDeltaYears < 1) { rlvmax <- rlvlpeam } else { rlvmax <- rlvlpeax } varargin=NA RLtstep <- tsEvaPlotReturnLevelsGPD(epsilon, sigma, threshold, epsilonStdErr, sigmaStdErr, thresholdStdErr, nPeaks = nPeaks, timeHorizonInYears = timeHorizonInYears, rlvmax, tstamps, trans,varargin) RLtstep ``` ### Curve bean (all curves from 1951 to 2020) ```{r} #| label: plotBeamD #| warning: false #| fig.width: 12 #| fig.height: 6 #| echo: false Plot2 = tsEvaPlotReturnLevelsGPDFromAnalysisObj(nonStationaryEvaParams, stationaryTransformData, timeIndex, trans=trans,ylabel="Discharge (m3/s)") Plot2 ``` ## Summary of the results for the Ardèche The table shows the 10 years return level of high flows (daily maximum flow) and low flows (7 days averaged flow) for 1951 and 2020. ```{r} #| label: tbl #| tbl-cap: "Changes in high and low flows" #| echo: false df <- rbind( data.frame( Y1951 = paste0(FloodRL$Y1951[1], " [", FloodRL$Y1951[2], " -", FloodRL$Y1951[3], "]"), Y2020 = paste0(FloodRL$Y2020[1], " [", FloodRL$Y2020[2], " - ", FloodRL$Y2020[3], "]"), change = paste0(FloodRL$change[1], " [", FloodRL$change[2], " - ", FloodRL$change[3], "]") ), data.frame( Y1951 = paste0(DryRL$Y1951[1], " [", DryRL$Y1951[2], " - ", DryRL$Y1951[3], "]"), Y2020 = paste0(DryRL$Y2020[1], " [", DryRL$Y2020[2], " - ", DryRL$Y2020[3], "]"), change = paste0(DryRL$change[1], " [", DryRL$change[2], " - ", DryRL$change[3], "]") ) ) names(df)=c("RL in 1951 (m3/s) [68% CI]","RL in 2020 (m3/s) [68% CI]","Relative change (%) [68% CI]") # Create the table kable(df, caption = "Summary table") ``` Although the river displays a non-linear trend, the transformed-stationary method estimates a slight decrease in both extreme high an low flows on the Ardèche between 1951 and 2020.