# R Code to simulate Harris County Intensities (72-hour) # # Notes: (1) Set R to correct directory (gloM and hy35 source) # (2) Put proper L-moments into the simulation calls # (3) Paste this entire file into R and run # # Load the L-moment package from CRAN # Attach as a library ## library(lmomco) # Attach Global Maximum and TP-40 data sets ## gloM<-read.csv(file="global_maxima.csv",header=T); attach(gloM) ## hy35<-read.csv(file="hydro35.csv",header=T); attach(hy35) # Quantile Functions for Depth and Duration simulation # Asquith and others, 2006, Eqns 13, and 14. # sim_depth<-function(f,p1,p2,p3,p4){(p1+(p2/p3)*(1-((1-f^p4)/p4)^p3))} sim_duration<-function(f,p1,p2,p3,p4){(p1+(p2/p3)*(1-((1-f^p4)/p4)^p3))} # # Vector fraction function for finding percentiles vecFrac<-function(arg1,arg2){ # R function to compare element-by-element two vectors and # compute proportion elements of vector2 that are larger # than vector1. # # Theodore G. Cleveland, Ph.D., P.E. 2007_1218. irow<-0 nrow<-length(arg1) for(i in 1:nrow){ if(arg2[i] >= arg1[i]) irow<-irow+1} frac<-irow/nrow return(frac)} # # Code Below For Statistical Simulations # # generate 2500 random probabilities frain<-runif(5000,0,1) fdur<-runif(5000,0,1) # L-moments for each station from Appendix 4, Asquith and others, 2006 # Station 3335, 6-hour inter-event arrival time lmrain<-vec2lmom(c(0.70202, 0.45144, 0.54276, 0.34056 )) lmdur<-vec2lmom(c(6.1336, 3.1726, 0.40652, 0.17418 )) # Station 3335, 72-hour inter-event arrival time lmrain<-vec2lmom(c(0.70202, 0.45144, 0.54276, 0.34056 )) lmdur<-vec2lmom(c(36.722, 22.091, 0.41453, 0.16234 )) # get Kappa parameters from L-moments parrain<-lmom2par(lmrain,type="kap") pardur<-lmom2par(lmdur,type="kap") # generate depths associated with probabilities dep<-sim_depth(frain,parrain$para[1],parrain$para[2],parrain$para[3],parrain$para[4]) # generate durations associated with probabilities dur<-sim_duration(fdur,pardur$para[1],pardur$para[2],pardur$para[3],pardur$para[4]) # Proportion values by Empirical Hyetograph dep1<- 0.216 * dep dep2<- 0.3757 * dep dep3<- 0.5155 * dep dep4<- 0.6304 * dep dep5<- 0.7166 * dep dep6<- 0.7738 * dep dep7<- 0.8089 * dep dep8<- 0.8332 * dep dep9<- 0.8501 * dep dep10<- 0.8635 * dep dep11<- 0.8766 * dep dep12<- 0.8896 * dep dep13<- 0.9018 * dep dep14<- 0.9129 * dep dep15<- 0.9225 * dep dep16<- 0.9305 * dep dep17<- 0.9372 * dep dep18<- 0.9424 * dep dep19<- 0.9464 * dep dep20<- 0.9492 * dep dep21<- 0.9518 * dep dep22<- 0.954 * dep dep23<- 0.957 * dep dep24<- 0.9606 * dep dep25<- 0.9647 * dep dep26<- 0.969 * dep dep27<- 0.9732 * dep dep28<- 0.9768 * dep dep29<- 0.9797 * dep dep30<- 0.9819 * dep dep31<- 0.9838 * dep dep32<- 0.9856 * dep dep33<- 0.9872 * dep dep34<- 0.989 * dep dep35<- 0.9909 * dep dep36<- 0.9929 * dep dep37<- 0.9949 * dep dep38<- 0.997 * dep dep39<- 0.9992 * dep # dur1<-0.025 * dur dur2<-0.05 * dur dur3<-0.075 * dur dur4<-0.1 * dur dur5<-0.125 * dur dur6<-0.15 * dur dur7<-0.175 * dur dur8<-0.2 * dur dur9<-0.225 * dur dur10<-0.25 * dur dur11<-0.275 * dur dur12<-0.3 * dur dur13<-0.325 * dur dur14<-0.35 * dur dur15<-0.375 * dur dur16<-0.4 * dur dur17<-0.425 * dur dur18<-0.45 * dur dur19<-0.475 * dur dur20<-0.5 * dur dur21<-0.525 * dur dur22<-0.55 * dur dur23<-0.575 * dur dur24<-0.6 * dur dur25<-0.625 * dur dur26<-0.65 * dur dur27<-0.675 * dur dur28<-0.7 * dur dur29<-0.725 * dur dur30<-0.75 * dur dur31<-0.775 * dur dur32<-0.8 * dur dur33<-0.825 * dur dur34<-0.85 * dur dur35<-0.875 * dur dur36<-0.9 * dur dur37<-0.925 * dur dur38<-0.95 * dur dur39<-0.975 * dur # # Build the plot vectors depth<-rbind(dep1,dep2,dep3,dep4,dep5,dep6,dep7,dep8,dep9,dep10,dep11,dep12,dep13,dep14,dep15,dep16,dep17,dep18,dep19,dep20,dep21,dep22,dep23,dep24,dep25,dep26,dep27,dep28,dep29,dep30,dep31,dep32,dep33,dep34,dep35,dep36,dep37,dep38,dep39) duration<-rbind(dur1,dur2,dur3,dur4,dur5,dur6,dur7,dur8,dur9,dur10,dur11,dur12,dur13,dur14,dur15,dur16,dur17,dur18,dur19,dur20,dur21,dur22,dur23,dur24,dur25,dur26,dur27,dur28,dur29,dur30,dur31,dur32,dur33,dur34,dur35,dur36,dur37,dur38,dur39) # build the plot base plot(duration,depth/duration,log="xy",col=rgb(0,0,1,0.5),xlab="Duration (hours)",ylab="Intensity (in/hr)",xaxp=c(.1,1000,3),yaxp=c(.1,1000,3),xlim=c(1/60,240),ylim=c(0.02,120),tck=.01,pch=21,type="p") ## plot(dur,dep/dur,log="xy",col=rgb(0,0,1,0.25),xlab="Duration (hours)",ylab="Intensity (in/hr)",xaxp=c(.1,1000,3),yaxp=c(.1,1000,3),xlim=c(1/60,240),ylim=c(0.02,120),tck=1,pch=16) # add the TP-40 lines(h35_dur,h35_dep/h35_dur,type="p",pch=21,col=rgb(1,0,0,1),cex=1.5) # add Global max lines(glo_dur,glo_dep/glo_dur,type="p",pch=16,col=rgb(1,0,0,1),cex=1.5) # # Find the equivalent percentile from TxDOT EBDLKUP.xls # Need to monkey with time units (minutes 2 hours) timeD<-seq(1,240*60,1)/60 #jPred=(91/60)/((duration+.1316667)^.706); vecFrac(depth/duration,jPred) jPred=(91)/((duration*60+7.9)^.706) vecFrac(depth/duration,jPred) intTxDOT=(91)/( ((timeD*60)+7.9) ^.706 ) lines(timeD,intTxDOT,col=rgb(1,0,0,1),lwd=3) gloMax=(16.6/(timeD+0.0)^(1-0.475)) lines(timeD,gloMax,col=rgb(0,1,0,1),lwd=3) emp99<-(6/(timeD+0.2)^(1)) emp50<-(.4/(timeD+0.2)^(1)) jPred=(6)/((duration)^1) vecFrac(depth/duration,jPred) jPred=(.4)/((duration)^1) vecFrac(depth/duration,jPred) lines(timeD,emp99,col=rgb(0,1,1,1),lwd=3) lines(timeD,emp50,col=rgb(0,1,1,1),lwd=3)