This assignment was designed to incorporate our knowledge of Hydrology and data manipulation skills with a free online software package ‘RStudio’. The task was to generate a design storm for the upper catchment of the Mary River by importing IFD data from the Bureau of Meteorology and by using various techniques to apply the necessary temporal patterns and areal reduction factors to the data for the final product to be produced.
The details of each step will be outlined below through Rstudio’s KnitR package, which will allow the use of imbedded code and graphs from R, and a final reproducible html file format for the report.
For the Design Storm to be created, the following steps were performed:
Although there were some problems (to be expected with a free software) Rstudio proved to be a highly efficient and effective program for data interpretation and analysis.
This report will look at the methods and methodology behind each step, and will explore in detail, the coding implemented.
The following ‘add-on’ packages were used throughout:
The first step in the process was to produce a map with a visible outline of the catchment with a clearly visible centroid for visualisation. The appropriate coding was provided through the assignment instructions, and with a few changes, the final code has been provided below.
By accessing the Queensland governments Wetlands Information Server, a KML file of the Mary Drainage Basin (outline) could be obtained, and when converted into a zipped KMZ file, the outline could be extracted and plotted over the top of the projected map.
And by plotting the centroid point on the final product, the map was created. Reading from the sc.rg package (Spatial Polygons Data Frame) in the Environment, the Centroid of the Catchment was found. The map of the area has been provided below.
library("OpenStreetMap")
## Loading required package: rJava
## Loading required package: raster
## Loading required package: sp
## Loading required package: rgdal
## rgdal: version: 0.8-16, (SVN revision 498)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgdal/proj
library("rgdal")
setwd("/Users/jordanandrews/Desktop/Subcatchments/")
sc.rg<-readOGR(".","Subcatchments",verbose=FALSE)
mary_map<-openmap(c(-25.29, 151.97),c(-26.86, 153.2),type="esri-topo")
new_mary_map<-openproj(mary_map,projection= proj4string(sc.rg))
plot(new_mary_map)
plot(sc.rg,add=TRUE, border="red",lwd=1.5)
dir_kml<-"/Users/jordanandrews/Desktop/Subcatchments/"
fkml<-"basin-mary.kml"
kmlname<-"Mary drainage basin"
kmlplot<-readOGR(paste(dir_kml,fkml,sep="/"),kmlname,verbose=FALSE)
mrb<-spTransform(kmlplot,CRS(proj4string(sc.rg)))
plot(mrb,add=T,border="blue")
points(470599.3,7054797.8,col="red",pch=20)
The following is a large chunk of code which was also provided in the assignment outline, it is used to produce the IFD chart which is produced by the Bureau of Meteorology Site when entering the centroid of the catchment.
By extracting the table of rainfall data from the site into a text file in R, the code below was used to read the table and plot the duration over the log of intensity.
The first section is used to manipulate the table of data extracted to make it easier to plot using the ggplot function below, while the last line of code contains many elements of the ggplot function, enabling the graph to be a replica of the Bureau.
A plot of the final product has been provided below. This IFD chart represents the log of intensity of rainfall over the specified durations, for each ARI for the Mar River Catchment.
library("reshape")
library("ggplot2")
dirv<-"/Users/jordanandrews/Desktop/Subcatchments/"
fname<-"IFDData"
ifd<-read.table(paste(dirv,fname,sep="/"),sep=",",header=T)
ifd[1:5,9]<-as.numeric(gsub("Mins","",ifd[1:5,1]))
ifd[6,9]<-as.numeric(gsub("Hr","",ifd[6,1]))*60
ifd[7:13,9]<-as.numeric(gsub("Hrs","",ifd[7:13,1]))*60
names(ifd)<-c("xlab","ari1","ari2","ari5","ari10","ari20","ari50","ari100","mins")
xlab<-ifd$xlab
bmins<-ifd$mins
seq_v<-c(1,2,3,4,6,7,8,10)
ylab<-c(seq_v[3:9]/10,seq_v,seq_v*10,seq_v[1:6]*100)
brain<-ylab
ifdl<-melt(ifd[,2:8])
## Using as id variables
ifdl$mins<-ifd[,9]
attach(ifdl)
colvec=c("dark blue","dark green","brown","light blue","red","purple","green")
labvec=c("1 Year(lower curve)","2 Years","5 Years","10 Years","20 Years","50 Years","100 Years(upper curve)")
ifdl $ variable <-as.factor(ifdl$variable)
ifdp<-ggplot(ifdl)+geom_line(aes(x=mins,y=value,group=variable,colour=variable))+
scale_x_log10(limits=c(5,4320),breaks=bmins,labels=xlab)+scale_y_log10(limits=c(0.3,1000),breaks=brain,labels=ylab)+
xlab("DURATION IN MINUTES OR HOURS")+ylab("RAINFALL INTENSITY IN MILLIMETRES PER HOUR")+
scale_color_manual(name="AVERAGE RECURRENCE INTERVAL",values=colvec,labels=labvec)+
theme(legend.justification = c(0, 0),legend.position = c(0, 0),panel.background = element_rect(colour = "black",fill="white"),panel.grid.major = element_line(colour = 'grey'),panel.grid.minor = element_line(colour = 'grey'),plot.title = element_text(lineheight=.8, face="bold"))+
guides(col = guide_legend(reverse = TRUE))+
ggtitle("DESIGN RAINFALL INTENSITY CHART\n Mary River Basin BoM weather station location")
ifdp
The next step was to import the temporal pattern data for the catchment which was provided as a large Output RIFD comma separated value file. This was imported into Rstudio and the following lines of code were used to make it usable to plot, similar to the process used for the IFD table. By sub setting a column with either 3 _ 1 or 3 _ 2 values the ggplot function was able to plot the perc_v column with the two sub setted columns, creating the initial plot of the Temporal Pattern.
As seen from the coding above, a bar and line plot have been created for each ARI: Below is the bar graph and line graph for the ARI1 (>30)
library("ggplot2")
dirn<-"/Users/jordanandrews/Desktop/Subcatchments/"
fname<-"Output4RIFD.csv"
tp<-read.csv(paste(dirn,fname,sep=""),sep=",",header=T)
tp$groupv <-paste(tp$zone,tp$ari,sep="_")
tpss1<-subset(tp,tp$ groupv=="3_1"&tp$dur_min==720)
tpp1<-ggplot(data=tpss1,aes(x=period_num, y=perc_v))+geom_bar(stat="identity")
tppl1<-ggplot(data=tpss1,aes(x=period_num, y=perc_v))+geom_line()
tpp1
tppl1
By now applying the temporal pattern data to the perc_v column, the intensity of rainfall can now be plotted to show the amount of rainfall that fell for each time period over the period of the storm. This plot has not yet had the areal reduction factor applied to it yet.
The formula below (tpss1 $ perc_ v/100 * ifdv * tpss1$time_step/60) is used to apply the temporal pattern (percentage of rain that fell) to the intensity of the rain (from the ifd data) but by also changing the time step value into minutes by dividing by 60, and by multiplying by a set ifdv value of 100. The code has been provided below and the two plots for ARI1 and ARI2. The tables were simply plotted as they are still to be reduced by the areal reduction factor, they are for viewing purposes only.
ifdv<-100
ds1<-tpss1$perc_v/100*ifdv*tpss1$time_step/60
barplot(ds1)
tpss2<-subset(tp,tp$groupv=="3_2"&tp$dur_min==720)
ds2<-tpss2$perc_v/100*ifdv*tpss2 $ time_step/60
barplot(ds2)
The Areal Reduction Factor can now be applied for long duration to the temporal pattern plots above, the tables are reduced slighty due to the size of the catchment. Larger sized catchments will be reduced by a lot more, but as the area is only approx.. 18.22km2, it will only be a factor of approx.. 0.9. The area was also read from the sc.rg package in the environment, where $AreaM2 = 18222759.
The equation for reducing a catchment over a long duration is:
ARF = Min{1,[1 – 0.2257(Area0.1685 – 0.8306log10Duration)Duration-0.3994]}
library("knitcitations")
library("RefManageR")
citep("10.7158/w11-858.2012.16.1")
## [1] "(French and Jones, 2012)"
Where the Area = 18.22km2 and Duration = 12hrs.
The code has been provided below and the (ds1red) and (ds2red) tables have been bar plotted as the final design storms. The plots were given axis and heading titles using the code below, and 2 plots have been included in this report for ARI>30 and ARI<30.
dirn<-"/Users/jordanandrews/Desktop/Subcatchments/"
fname<-"Output4RIFD.csv"
tp<-read.csv(paste(dirn,fname,sep=""),sep=",",header=T)
tp$groupv<-paste(tp$zone,tp$ari,sep="_")
tpss1<-subset(tp,tp$groupv=="3_1"&tp$dur_min==720)
ifdv<-100
ds1<-tpss1$perc_v/100*ifdv*tpss1$time_step/60
ds1red<-ds1*1-0.2257*(18.22^(0.1685)-0.8306*log10(12))*12^0.3994
barplot(ds1red,xlab="time", ylab="rainfall(cm)",main="Design Storm for Mary River Catchment ARI > 30 for Duration=12hr", names.arg=("0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5 10 10.5 11 11.5 12"))
dirn<-"/Users/jordanandrews/Desktop/Subcatchments/"
fname<-"Output4RIFD.csv"
tp<-read.csv(paste(dirn,fname,sep=""),sep=",",header=T)
tp$groupv<-paste(tp$zone,tp$ari,sep="_")
tpss2<-subset(tp,tp$ groupv=="3_2"&tp$dur_min==720)
ifdv<-100
ds2<-tpss2$perc_v/100*ifdv*tpss2$time_step/60
ds2red<-ds2*1-0.2257*(18.22^(0.1685)-0.8306*log10(12))*12^-0.3994
barplot(ds2red, xlab = "time (hrs)", ylab = "rainfall (cm)", main = "Design Storm for Mary River Catchment ARI < 30 for Duration=12hr", names.arg=("0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5 10 10.5 11 11.5 12"))
Although not used for this Design Storm, filtering is normally required as part of calibration and verification of results. Filtering is often used to modify the results and data due to the type of catchment being dealt with, as not all rainfall is ‘direct.’ (Thompson, R. 2013) This may be due to outside factors influencing the water measurement results, such as non-vertical rain, runoff from outside the catchment or foliage reducing water catchment.
**Bibliography**
R. French and M. Jones, 2012. “Design rainfall temporal patterns in Australian Rainfall and Runoff: Durations exceeding one hour”. In: \(\lbrace\)AJWR\(\rbrace\) 16.1 (2012). DOI: 10.7158/w11-858.2012.16.1. .
R. Thompson, 2013. “Direct Rainfall Method: Australian Rainfall and Runoff/Runoff Estimation.” From: http://www.arr.org.au/wp-content/uploads/2013/Presentations/Direct%20Rainfall%20Presentation.pd