Skip to content

Instantly share code, notes, and snippets.

@TonyLadson
Last active November 14, 2016 00:43
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save TonyLadson/fad3552c0ebe7bc87810a2f6a231effb to your computer and use it in GitHub Desktop.
Save TonyLadson/fad3552c0ebe7bc87810a2f6a231effb to your computer and use it in GitHub Desktop.
---
title: "RORB: comparing modelled and measured hydrographs."
output: html_notebook
---
# Introduction
This is the [R Markdown](http://rmarkdown.rstudio.com) Notebook to accompany the blog: [RORB: comparing modelled and measured hydrographys](https://tonyladson.wordpress.com/2016/11/14/rorb-comparing-modelled-and-measured-hydrographs/)
```{r setup}
library(tidyverse)
```
# Data
Example calculations are made using data from the Werribee River worked example (Section 10.2 of the [RORB manual](https://dl.dropboxusercontent.com/u/10963448/rorb_manual.pdf)). The actual and calculated hydrographs are listed on page 111 and 112.
```{r example_data}
calculated <- c(0,0,8.007,37.806,92.598,163.077,236.298,302.075,
337.999,327.515,292.287,254.062,216.021,178.231,
142.87,117.351,100.052,84.542,72.073,64,58.707,54.95,
52.309,50.392,47.78,43.812,40.403,38.49,37.407)
actual <- c(0,0,8,34,64,147,245,310,356,330,290,245,216,185,150,122,104,
96,90,76,68,62,59,57,55,53,50,42,36)
```
# Peak discharge
The peak discharge is just the maximum discharge that occurs within the hydrograph.
```{r peak_discharge}
# Hydrograph peaks
max(actual)
max(calculated)
# Absolute difference in peaks
max(calculated) - max(actual)
# Percentage difference in peaks
100 * (max(calculated) - max(actual))/max(actual)
```
# Time to peak
For the Werribee case study, hydrographs start at time step zero, end at time step 28 and the time increment is 2 hours.
```{r time_to_peak}
time_start <- 0
time_end <- 28
time_inc <- 2
total_inc <- time_end - time_start + 1
tt <- time_inc*(time_start:time_end)
# time of peak for actual hydrograph
peak_time_actual <- tt[which.max(actual)]
peak_time_actual
# time of peak for calculated hydrograph
peak_time_calculated <- tt[which.max(calculated)]
peak_time_calculated
# Absolute difference
peak_time_calculated - peak_time_actual
# Percentage difference
100*(peak_time_calculated - peak_time_actual)/peak_time_actual
```
# Volume
RORB bases the volume on the average flow times the length of the hydrograph
```{r volume}
# Mean flow, actual hydrograph
mean(actual)
# Mean flow, calculated hydrograph
mean(calculated)
# Volume actual
vol_actual <- mean(actual) * total_inc * time_inc * 60 * 60
vol_actual
format(vol_actual, scientific = TRUE)
# Volume calculated
vol_calculated <- mean(calculated) * total_inc * time_inc * 60 * 60
vol_calculated
format(vol_calculated, scientific = TRUE)
# Absolute error
vol_error <- vol_calculated - vol_actual
format(vol_error, scientific = TRUE)
# Percentage error
100 * vol_error/vol_actual
```
# Average absolute coordinate error
The average absolute coordinate error is the mean of the absolute values of the differences between the actual and calculated hydrographs.
```{r ave_abs_coord_error}
# Absolute error
abs_error <- mean(abs(actual-calculated))
abs_error
# Percentage error
100*abs_error/mean(actual)
```
# Time to centroid
These calculations match the results produced by RORB.
```{r time_to_centroid}
# Time to centroid actual
centroid_actual <- sum(tt*actual)/sum(actual)
centroid_actual
# Time to centroid calculated
centroid_calculated <- sum(tt*calculated)/sum(calculated)
centroid_calculated
# Absoluate error
sum(tt*calculated)/sum(calculated) - sum(tt*actual)/sum(actual)
# Percentage error
100*(sum(tt*calculated)/sum(calculated) - sum(tt*actual)/sum(actual) )/(sum(tt*actual)/sum(actual))
```
# Lag (centre of mass to centre of mass)
```{r lag_cm}
# In the Werribee case, the concentrated input is the hydrograph at Melton Reservoir
h_melton <- c(0,0,66,150,253,325,391,420,309,247,211,166,139,88,
86,82,63,55,54,52,50,49,48,47,37,36,36,36,36)
# Time to centroid actual
sum(tt*actual)/sum(actual)
# Time to centroid calculated
sum(tt*calculated)/sum(calculated)
# Time to centroid input
sum(tt*h_melton)/sum(h_melton)
# difference: input to actual
sum(tt*actual)/sum(actual) - sum(tt*h_melton)/sum(h_melton)
# difference: input to calculated
sum(tt*calculated)/sum(calculated) - sum(tt*h_melton)/sum(h_melton)
# Absolute difference
abs_err <- (sum(tt*calculated)/sum(calculated) - sum(tt*h_melton)/sum(h_melton)) -
(sum(tt*actual)/sum(actual) - sum(tt*h_melton)/sum(h_melton))
abs_err
# Pecentage difference
100*abs_err/(sum(tt*actual)/sum(actual) - sum(tt*h_melton)/sum(h_melton) )
```
# Lag to peak
```{r lag_to_peak}
# Time to centroid inputs
sum(tt*h_melton)/sum(h_melton)
# Time of peak (actual)
tt[which.max(actual)]
# Lag to peak
# Actual
tt[which.max(actual)] - sum(tt*h_melton)/sum(h_melton)
# Calculated
tt[which.max(calculated)] - sum(tt*h_melton)/sum(h_melton)
# Absolute error
abs.err <- (tt[which.max(calculated)] - sum(tt*h_melton)/sum(h_melton)) -
(tt[which.max(actual)] - sum(tt*h_melton)/sum(h_melton))
# Percentage error
abs.err/(tt[which.max(actual)] - sum(tt*h_melton)/sum(h_melton))
```
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment