Load required packages

library(dplyr)
library(SciViews)
library(ggplot2)

Specify start and ending year for the analysis

start=1948
end=2017

Read the daily rainfall data for a specific meteorological station and omit the rows containing “Feb 29”

daily_data<- read.csv("Coxbazar.csv")
#summary(daily_data)
###Omit leap year values
leap<-which(daily_data$Day_month=="29-Feb")
daily_data<-daily_data[-leap,]

Looking at the structure of the daily rainfall dataset

head(daily_data)
##   Day Day_month Month Year TSVALUE_UPDATED
## 1   1    01-Jan   Jan 1948               0
## 2   2    02-Jan   Jan 1948               0
## 3   3    03-Jan   Jan 1948               0
## 4   4    04-Jan   Jan 1948               0
## 5   5    05-Jan   Jan 1948               0
## 6   6    06-Jan   Jan 1948               0

Constructing required tables, At first, Duration: 1 day

First “one_day_table” is created as a NULL column dataframe having 365 values. It is populated by running a loop. During loop operation, each time a year is added as a column. The columns represent the years and there are 365 rows.

“One day max” is a single row dataframe having the maximum daily rainfall values for the individual years.

“one_day” was formed from “rbind”-ing” of “one_day_table” and “one_day_max”

one_day_table<- data.frame(rep(NA,365))

for (i in start:end)
{
  data<-daily_data[daily_data$Year==i,]
  data<-data.frame(data$TSVALUE_UPDATED)
  
  one_day_table<-cbind(one_day_table,data)

}
one_day_table<-one_day_table[-1]
colnames(one_day_table)<-seq(start,end)

one_day_max<-one_day_table %>% summarise_if(is.numeric, max, na.rm=TRUE)
one_day<-rbind(one_day_table,one_day_max)

Duration: 2 day

First “two_day_table” is created as a NULL column dataframe having 365 values. It is populated by running a loop. The columns represent the years and there are 365 rows.

“two_day_lag” is calculated as the cumulative sum of two days rainfall in a rolling manner, all over a specific year. For an individual year (column) there is one null Value at the beginning. Calculation for the 2-day duration is done in this manner within a loop. During loop operation, each time a year is added as a column. Finally, it is appended with the blank “two_day_table”.

“two_day max” is a single row dataframe having the maximum 2-day rainfall values for the individual years.

“two_day” was formed from “rbind”-ing” of “two_day_table” and “two_day_max”

two_day_table <- data.frame(rep(NA,365))

for (i in names(one_day_table)){
  
  two_day_lag<-data.frame(x=RcppRoll::roll_sum(one_day_table[, i],2, fill=NA, align="right"))
  two_day_table<-cbind(two_day_table,two_day_lag)
}
two_day_table<-two_day_table[-1]
colnames(two_day_table)<-seq(start,end)

two_day_max<-two_day_table %>% summarise_if(is.numeric, max, na.rm=TRUE)
two_day<-rbind(two_day_table,two_day_max)

Duration: 3 day

First “three_day_table” is created as a NULL column dataframe having 365 values. It is populated by running a loop. The columns represent the years and there are 365 rows.

“three_day_lag” is calculated as the cumulative sum of three days rainfall in a rolling manner, all over a specific year. For an individual year (column) there are two null Values at the beginning. Calculation for the 3-day duration is done in this manner within a loop. During loop operation, each time a year is added as a column. Finally, it is appended with the blank “three_day_table”.

“three_day max” is a single row dataframe having the maximum 3-day rainfall values for the individual years.

“three_day” was formed from “rbind”-ing” of “three_day_table” and “three_day_max”

three_day_table <- data.frame(rep(NA,365))

for (i in names(one_day_table)){
  
  three_day_lag<-data.frame(x=RcppRoll::roll_sum(one_day_table[, i],3, fill=NA, align="right"))
  three_day_table<-cbind(three_day_table,three_day_lag)
}
three_day_table<-three_day_table[-1]
colnames(three_day_table)<-seq(start,end)

three_day_max<-three_day_table %>% summarise_if(is.numeric, max, na.rm=TRUE)
three_day<-rbind(three_day_table,three_day_max)

Duration: 5 day

First “five_day_table” is created as a NULL column dataframe having 365 values. It is populated by running a loop. The columns represent the years and there are 365 rows.

“five_day_lag” is calculated as the cumulative sum of five days rainfall in a rolling manner, all over a specific year. For an individual year (column) there are four null Values at the beginning. Calculation for the 5-day duration is done in this manner within a loop. During loop operation, each time a year is added as a column. Finally, it is appended with the blank “five_day_table”.

“five_day max” is a single row dataframe having the maximum 5-day rainfall values for the individual years.

“five_day” was formed from “rbind”-ing” of “five_day_table” and “five_day_max”

five_day_table <- data.frame(rep(NA,365))

for (i in names(one_day_table)){
  
  five_day_lag<-data.frame(x=RcppRoll::roll_sum(one_day_table[, i],5, fill=NA, align="right"))
  five_day_table<-cbind(five_day_table,five_day_lag)
}
five_day_table<-five_day_table[-1]
colnames(five_day_table)<-seq(start,end)

five_day_max<-five_day_table %>% summarise_if(is.numeric, max, na.rm=TRUE)
five_day<-rbind(five_day_table,five_day_max)

Duration: 10 day

First “ten_day_table” is created as a NULL column dataframe having 365 values. It is populated by running a loop. The columns represent the years and there are 365 rows.

“ten_day_lag” is calculated as the cumulative sum of ten days rainfall in a rolling manner, all over a specific year. For an individual year (column) there are nine null Values at the beginning. Calculation for the 10-day duration is done in this manner within a loop. During loop operation, each time a year is added as a column. Finally, it is appended with the blank “ten_day_table”.

“ten_day max” is a single row dataframe having the maximum 10-day rainfall values for the individual years.

“ten_day” was formed from “rbind”-ing” of “ten_day_table” and “ten_day_max”

ten_day_table <- data.frame(rep(NA,365))

for (i in names(one_day_table)){
  
  ten_day_lag<-data.frame(x=RcppRoll::roll_sum(one_day_table[, i],10, fill=NA, align="right"))
  ten_day_table<-cbind(ten_day_table,ten_day_lag)
}
ten_day_table<-ten_day_table[-1]
colnames(ten_day_table)<-seq(start,end)

ten_day_max<-ten_day_table %>% summarise_if(is.numeric, max, na.rm=TRUE)
ten_day<-rbind(ten_day_table,ten_day_max)

Constructing a table of maximum values of rainfall.

The rows represent each year, starting from start to end year, top to bottom; and the columns represent each duration.

final_max_table<- data.frame(rep(NA, end-start+1))  

final_max_table<-cbind(final_max_table,t(one_day[366,]))
final_max_table<-cbind(final_max_table,t(two_day[366,]))
final_max_table<-cbind(final_max_table,t(three_day[366,]))
final_max_table<-cbind(final_max_table,t(five_day[366,]))
final_max_table<-cbind(final_max_table,t(ten_day[366,]))

final_max_table<-final_max_table[-1]

colnames(final_max_table)<-c("one day","two day","three day","five day","ten day")

Converting values from “day” scale to “hourly” scale.

final_max_table_hr<- cbind(final_max_table$`one day`/24,final_max_table$`two day`/(24*2),
                           final_max_table$`three day`/(24*3), final_max_table$`five day`/(24*5),
                           final_max_table$`ten day`/(24*10))

colnames(final_max_table_hr)<-c("one day","two day","three day","five day","ten day")

Calculate the Standard Deviation of the Maximum rainfall values.

Each column of the “final_max_table_hr” represents a vector of Maximum values of 1,2,3,5 and 10 days rainfall (converted to hours)

sd_vec <- apply(final_max_table_hr, 2, sd)

Calculate the Mean of the Maximum rainfall values.

Each column of the “final_max_table_hr” represents a vector of Maximum values of 1,2,3,5 and 10 days rainfall (converted to hours)

mean_vec <- apply(final_max_table_hr, 2, mean)

Create a vector of Return Period.

return_periods<- c(5,10,30,50,100)

Calculate Frequency Factor from individual Return Periods.

freq_factor <- function(t) {
  KT =  -sqrt(6)/pi * ( 0.5772 + ln(ln(t/(t-1)) ))
  
}

kT<- freq_factor(return_periods)

Calculate XTs - (value of the random variable associated with given return periods

one_day_XT<- mean_vec[1] + kT * sd_vec[1]
two_day_XT<- mean_vec[2] + kT * sd_vec[2] 
three_day_XT<- mean_vec[3] + kT * sd_vec[3] 
five_day_XT<- mean_vec[4] + kT * sd_vec[4] 
ten_day_XT<- mean_vec[5] + kT * sd_vec[5] 

final<- rbind(one_day_XT, two_day_XT, three_day_XT, five_day_XT, ten_day_XT)
final<- data.frame(final)

Plot IDF Curve

days<-c(1,2,3,5,10)

final<- cbind(days,final)
colnames(final)<- c("days",  "RP5year","RP10year","RP30year","RP50year","RP100year")


IDF_curve<- ggplot(data = final, aes(x = days)) +
  geom_line(aes(y = RP5year, colour = "RP 5 year"),lwd=2) +
  geom_line(aes(y = RP10year, colour = "RP 10 year"),lwd=2) +
  geom_line(aes(y = RP30year, colour = "RP 30 year"),lwd=2) +
  geom_line(aes(y = RP50year, colour = "RP 50 year"),lwd=2) +
  geom_line(aes(y = RP100year, colour = "RP 100 year"),lwd=2) +
  
  
  scale_colour_manual("", 
                      breaks = c("RP 5 year", "RP 10 year","RP 30 year","RP 50 year","RP 100 year"),
                      values = c("darkblue", "red", "orange", "black", "green")) +
  labs(y = "Rainfall Intensity", title = "IDF Curve")  

IDF_curve