Influence of catchment shape factor on regional flood estimation https://tonyladson.wordpress.com/2016/04/13/shape-factor-in-regional-flood-estimation/

ShapeFactor.R

####################################################################
#
# Influence of shape factor in regional flood estimation
# See https://tonyladson.wordpress.com/2016/04/13/shape-factor-in-regional-flood-estimation/
#
####################################################################


# Plot of shape factor distribution

my.p <- c(1, 10, 50, 90, 99)
my.p <- my.p/100
my.z <- qnorm(my.p, lower.tail = FALSE)
Sf <- c(0.32, 0.51, 0.76, 1.06, 1.51)


par(oma = c(1,2,0,0))
plot(my.z, Sf, log = 'y', xaxt = 'n', xlab = '', type = 'b', las = 1, ylab = 'Shape factor',
     pch = 21,
     bg = 'grey')
axis(side = 1, at = my.z, label = my.p*100)
title(xlab = '% small than')


### Ellipses

library(plotrix)

plot(c(0,2), c(-0.2,2), type="n", xaxt = 'n', yaxt = 'n', bty = 'n', xlab = '', ylab = '')
draw.ellipse(x = 1, y = 1, a = 1, b = 1)
points(1, 1, pch = 19)
points(2, 1, pch = 19)

segments(1,1, 1,2, col = 'grey')
segments(1,1, 2,1, col = 'grey')
text(x = 1, y = 1.5, labels = 'b', pos = 4 )
text(x = 1.5, y = 1, labels = 'a', pos = 1 )
text(x = 1, y = -0.2, labels = expression(paste('Area = ', pi, ' a b')))
text(1, 1, labels = 'Centroid', pos = 1)
text(1.8, 0.8, labels = 'Outlet', pos = 1)
segments(1.8, 0.8, 2,1, col = 'black')

# a direction is side to side
# b direction is up and down

# We need a function that will plot an ellipse that has the correct aspect
# ratio based on a and b.  Scaling factors can be calculated from the
# plot dimensions par(pin) and the limiting x and y coordinates.

# x.c, y.c coordinates of the center 
# a and b ellipse size
# a direction is side to side
# b direction is up and down
# scale.factor makes the elipse larger or smaller while retaining
# the aspect ratio


MyEllipse <- function(x.c = 0.5, y.c = 0.5, a, b, scale.factor = 1){
  
  my.pin <- par('pin')
  my.usr <- par('usr')
  
  h.scale = my.pin[1]/(my.usr[2] - my.usr[1])

  v.scale = my.pin[2]/(my.usr[4] - my.usr[3])
  
  a.s <- scale.factor * a / h.scale
  b.s <- scale.factor * b / v.scale
  
  
  plot(c(0,1), c(0,1), type="n", xaxt = 'n', yaxt = 'n', bty = 'n', xlab = '', ylab = '', xaxs = 'i', yaxs = 'i')
  points(x = x.c, y = y.c, pch = 19, cex = 0.2)
  draw.ellipse(x = x.c, y = y.c, a = a.s, b = b.s)
  
  segments(x.c, y.c, x.c + a.s, y.c, col = 'grey')
  #segments(x.c, y.c, y.c, y.c - b.s, col = 'grey')
  text(x = 0.5*(x.c + x.c + a.s), y = y.c, labels = 'a', pos = 1 )
  #text(x = x.c, y = 0.5*(y.c + y.c - b.s), labels = 'b', pos = 4 )
 
  
}


# plot ellipses with various scale factors

my.a <- c(0.32, 0.51, 0.76, 1.06, 1.51)
my.b <- 1/(pi * my.a)


# remove all the margins and plot as 1 row and 5 columns

op <- par(oma = c(0,0,0,0), 
          mfrow = c(5,1), 
          omi = c(0,0,0,0), 
          mar = c(0,0,0,0),
          mai = c(0,0,0,0))

for(a in my.a) {
  b <- 1/(pi * a)
  MyEllipse(0.5, 0.5, a, b, 1.5)
  text(0.2, 0.6, cex = 1.5, bquote(paste('S'['f']*' = ', .(a))))
  
}

par(op)

### Sensitivity of flood estimates to the shape factor.

sf <-
  structure(list(ShapeFactor = c(0.04, 0.07, 0.23, 0.29, 0.35, 0.44, 0.62, 0.73, 0.83, 1.72), 
                      flood = c(614, 457, 204, 174, 155, 133, 107, 95.4, 88, 54.6)), 
                     .Names = c("ShapeFactor", "flood" ), 
                      class = c("data.frame"), row.names = c(NA, -10L))

dev.off() # reset graph settings

par(oma = c(1, 2, 0, 0))
plot(flood ~ ShapeFactor, data = sf, type = 'b', las = 1, ylab = '1% flood estimate',
     pch = 21, 
     bg = 'grey', 
     ylim = c(0,650),
     xlab = 'Shape factor')

## log-log
par(oma = c(1, 2, 0, 0))
plot(flood ~ ShapeFactor, data = sf, type = 'b', las = 1, ylab = '1% flood estimate',
     pch = 21, 
     bg = 'grey', 
     ylim = c(50,650),
     xlab = 'Shape factor',
     log = 'xy')


# Comparison at 226209
# Moe River at Darnum

# Estimates from At site Frequency Analysis
# These are taken from RFFE. The 'output_nearby' file

# outlet
# lat = -38.21
# lon = 146
# Centroid
# lat = -38.1861
# lon = 145.966
# Flow record length = 40 years
# Mean annual rainfall = 1067.75
# Shape factor = 0.281
# Catchment area = 214

flood.atSite <- c(28.19, 41.61, 49.65, 56.69, 64.92, 70.51)
flood.atSite.LCL <- c(24.27, 36.93, 44.74, 50.12, 55.5, 58.31)
flood.atSite.UCL <- c(32.7, 47.4, 56.3, 64.7, 76.4, 86.5)

flood.rffe <- c(38.3, 66.2, 88.7, 113, 150, 181)
flood.rffe.LCL <- c(19.6, 35.1, 45.9, 56.2, 69.2, 79)
flood.rffe.UCL <- c(74.6, 125, 173, 230, 320, 411)


aep <- c(0.5, 0.2, 0.1, 0.05, 0.02, 0.01)

aep.labels <- 1/aep
aep.labels

my.z <- qnorm(aep, lower.tail = FALSE )


par(oma = c(1,2,0,0))
plot(my.z, flood.atSite, 
     xaxt = 'n', 
     log = 'y',
     las = 1, 
     xlab = 'AEP: 1 in Y year',
     ylab = 'Flood peak (cumec)',
     pch = 21,
     bg = 'grey',
     ylim = c(20, 450), 
     type = 'b')
axis(side = 1, at = my.z, labels = aep.labels)
lines(my.z, flood.atSite.LCL, lty = 2)
lines(my.z, flood.atSite.UCL, lty = 2)
points(my.z, flood.rffe, type = 'b', pch = 21, bg = 'blue', col = 'blue')
lines(my.z, flood.rffe.LCL, lty = 2, col = 'blue', lwd = 2)
lines(my.z, flood.rffe.UCL, lty = 2, col = 'blue', lwd = 2)

legend('topleft', legend = c('At-site', 'Regional'), pch = 19, col = c('grey','blue'), lty = 1, bty = 'n')