library(dplyr)
library(SciViews)
library(ggplot2)
start=1948
end=2017
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,]
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
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)
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)
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)
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)
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)
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")
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")
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)
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)
return_periods<- c(5,10,30,50,100)
freq_factor <- function(t) {
KT = -sqrt(6)/pi * ( 0.5772 + ln(ln(t/(t-1)) ))
}
kT<- freq_factor(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)
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