Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- library(tidyverse)
- rain_design <- 90.9
- t_pattern <- c(4.1, 8.0, 11.0, 23.3, 15.3, 8.1, 7.0, 7.0, 5.1, 6.1, 3.1, 1.9)
- t_step <- 0.5
- t_step_count <- 1:length(t_pattern)
- t_hour <- t_step_count * t_step
- hyeto <- t_pattern * rain_design/100
- hyeto_df <- data_frame(t_step_count, t_hour, rain = hyeto)
- hyeto_df
- # Chart the hyetograph
- hyeto_df %>%
- ggplot(aes( x = t_step_count )) +
- geom_bar(aes(y = rain), stat = 'identity', fill = 'blue', width = 1, colour = 'black') +
- ylab("Rainfall (mm)") +
- scale_x_continuous(name = "Time (hours)",
- breaks = c(t_step, t_step_count + t_step),
- labels = c(0, t_hour)) + # some fiddling here to get the x-axis labelling correct
- theme_gray(base_size = 14)
- # Function to remove initial loss from a hyetograph
- #
- # Arguments
- # Hyeto = hyetograph (mm of rain for 1 or more time steps)
- # IL = Initial Loss (mm)
- Remove_IL <- function(Hyeto, IL){
- # IL <- 1
- # Hyeto <- c(30, 7, 4, 3, 7, 20, 15, 3)
- # #
- if(sum(Hyeto) < IL) {
- warning('Initial loss is not satisfied')
- Hyeto[] <- 0
- return(Hyeto)
- }
- rain_cum <- cumsum(Hyeto)
- tzero <- min(which(rain_cum > IL)) - 1
- if(tzero == 0) {
- Hyeto[1] <- Hyeto[1] - IL
- return(Hyeto)
- }
- Hyeto[tzero + 1] <- Hyeto[tzero + 1] - (IL - rain_cum[tzero])
- Hyeto[1:tzero] <- 0
- Hyeto
- }
- # Function to remove continuing loss from a hyetograph
- # Arguments
- # Hyetograph - rainfall for each time step
- # CL - continuing loss in mm/h
- # t_step - length of each time step in hours
- Remove_CL <- function(Hyeto, CL, t_step) {
- CL_t_step = CL * t_step # Continuing loss per time step
- Hyeto <- Hyeto - CL_t_step
- Hyeto[Hyeto < 0] <- 0
- Hyeto
- }
- # Plot rainfall excess hyetograph
- hyeto_df %>%
- mutate(rain_IL = Remove_IL(rain, 10),
- rain_ILCL = Remove_CL(rain_IL, 2, 0.5)) %>%
- ggplot(aes( x = t_step_count)) +
- geom_bar(aes(y = rain), stat = 'identity', fill = 'white', width = 1, colour = 'black') +
- geom_bar(aes(y = rain_IL), stat = 'identity', fill = 'yellow', width = 1, colour = 'black') +
- geom_bar(aes(y = rain_ILCL), stat = 'identity', fill = 'blue', width = 1, colour = 'black') +
- ylab("Rainfall (mm)") +
- scale_x_continuous(name = "Time (hours)",
- breaks = c(t_step, t_step_count + t_step),
- labels = c(0, t_hour)) + # some fiddling here to get the x-axis labelling correct
- annotate(geom = 'rect', xmin = 8, xmax = 8.5, ymin = 20, ymax = 21, fill = 'white', colour = 'black') +
- annotate(geom = 'text', x = 8.6, y = 20.5, hjust = 0, label = 'Initial Loss' ) +
- annotate(geom = 'rect', xmin = 8, xmax = 8.5, ymin = 18.5, ymax = 19.5, fill = 'yellow', colour = 'black') +
- annotate(geom = 'text', x = 8.6, y = 19, hjust = 0, label = 'Continuing loss' ) +
- theme_gray(base_size = 14)
- # Convert to a hydrograph
- sub_area <- 78.7
- t_step = 0.5 # time step in hours
- # Example calculation
- hyeto_df %>%
- mutate(rain_IL = Remove_IL(rain, 10),
- rain_ILCL = Remove_CL(rain_IL, 2, 0.5)) %>%
- select(rain_ILCL) %>%
- unlist(use.names = FALSE) # pull out the rainfall excess
- # [1] 0.0000 0.0000 8.9990 20.1797 12.9077 6.3629 5.3630 5.3630 3.6359 4.5449 1.8179 0.7271
- # Rainfall in the 3rd time step is 8.9 mm
- (8.999 / 0.5) * 1e-3 * 78.7 * 1e+6 / 3600
- # Calculate the rainfall excess hydrograph
- # Add an extra row to the bottom of the data frame so the flow returns to zero
- # and plot hydrograph
- hyeto_df %>%
- mutate(rain_excess = Remove_CL(Remove_IL(rain, 10), 2, 0.5)) %>%
- add_row(t_step_count = 13,
- t_hour = 6.5,
- rain = 0,
- rain_excess = 0 ) %>%
- mutate(flow = (rain_excess/t_step) * sub_area * 1e3/3600) %>%
- ggplot(aes(x = t_hour, y = flow)) +
- geom_line(alpha = 0.5) +
- geom_point(shape = 21, color = 'grey', fill = 'blue', size = 3) +
- ylab(bquote('Flow ('*m^-3~s^-1*')')) +
- xlab('Time (hour)') +
- theme_gray(base_size = 14)
- hyeto_df %>%
- mutate(rain_excess = Remove_CL(Remove_IL(rain, 10), 2, 0.5)) %>%
- add_row(t_step_count = 13,
- t_hour = 6.5,
- rain = 0,
- rain_excess = 0 ) %>%
- mutate(flow = (rain_excess/t_step) * sub_area * 1e3/3600)
- # # A tibble: 13 × 5
- # t_step_count t_hour rain rain_excess flow
- # <dbl> <dbl> <dbl> <dbl> <dbl>
- # 1 1 0.5 3.7269 0.0000 0.00000
- # 2 2 1.0 7.2720 0.0000 0.00000
- # 3 3 1.5 9.9990 8.9990 393.45628
- # 4 4 2.0 21.1797 20.1797 882.30133
- # 5 5 2.5 13.9077 12.9077 564.35333
- # 6 6 3.0 7.3629 6.3629 278.20013
- # 7 7 3.5 6.3630 5.3630 234.48228
- # 8 8 4.0 6.3630 5.3630 234.48228
- # 9 9 4.5 4.6359 3.6359 158.96963
- # 10 10 5.0 5.5449 4.5449 198.71313
- # 11 11 5.5 2.8179 1.8179 79.48263
- # 12 12 6.0 1.7271 0.7271 31.79043
- # 13 13 6.5 0.0000 0.0000 0.00000
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement