Formula to calculate the equal area slope. See https://tonyladson.wordpress.com/2017/03/04/on-the-calculation-of-equal-area-slope/

Calc_slope_ea.R

# Description
# Calc_slope_ea(x,y) calculates teh equal area slope 


# Arguments
#
# x = distance along the stream in km
# y = elevation along the stream in m

# Details
# See references

# References
# Appendix A of
# Technical Memorandum No 61
# A method for estimating design peak discharge
# Planning and Technical Services
# Water and Soil Division
# Ministry of Works and Development
# Wellington 1980
# http://web.archive.org/web/20170228021947/http://docs.niwa.co.nz/library/public/TM61_1980.pdf



Calc_slope_ea <- function(x,y){
  
  y = y - min(y)
  
  # Function to calcuate area under a curve using the trapazoid rule
  # See http://svitsrv25.epfl.ch/R-doc/library/caTools/html/trapz.html
  .Calc_area_trapazoid <- function(x,y){
    idx = 2:length(x)
    return (as.double( (x[idx] - x[idx-1]) %*% (y[idx] + y[idx-1])) / 2)
    
  }
  
  A <- .Calc_area_trapazoid(x,y)
  2*A/(1000 * (tail(x,1) - head(x,1))^2)
}


# Test and playing around...

library(pracma)
library(tidyverse)


# Generate a randome stream profile
set.seed(17)

# The random profile doesn't always work so well but its sufficient for a test

profile_df <- data_frame(x = cumsum(runif(10, min = 0.1, max = 0.4)),
                         y = cumsum(c(runif(5, min = 0.1, max = 0.3), 
                                      runif(3, min = 0.3, max = 1),
                                      runif(2, min = 2, max = 4))))


profile_df <- profile_df %>% 
  mutate(y = y - min(y)) # Calculations are easier if we set the elevation of the outlet at zero
                           
plot_ea <- profile_df %>% 
  ggplot(aes(x,y)) + 
  geom_point(shape = 1, size = 2) + 
  geom_line() +
  scale_x_continuous(name = 'Distance (km)') +
  scale_y_continuous(name = 'Elevation (m)') +
  theme_gray(base_size = 14)

plot_ea

# Area under the curve

# Function to calcuate area under a curve using the trapazoid rule
# See http://svitsrv25.epfl.ch/R-doc/library/caTools/html/trapz.html

Calc_area_trapazoid <- function(x,y){
  idx = 2:length(x)
  return (as.double( (x[idx] - x[idx-1]) %*% (y[idx] + y[idx-1])) / 2)
  
}


A <- Calc_area_trapazoid(profile_df$x, profile_df$y)
L <- last(profile_df$x) - first(profile_df$x)
h = 2 * A/L

# C = equal area slope ordinate

C_x = max(profile_df$x) 
C_y = min(profile_df$y) + h

# Add to plot


plot_ea +
  annotate(geom = "point", x = C_x, y = C_y, size = 2, colour = 'blue') +
  annotate(geom = "segment", x = min(x), xend = max(x), y = min(y), yend = C_y, linetype = 'dashed')

# Next, define the polygons for the two parts
# plot and calculate the area 

# We need to find the intersection. 

# Expression for the y-ordinate of the profile

Profile_y <- approxfun(x = profile_df$x, y = profile_df$y) # expression for the 

# Expression for the y-ordinate of the equal area slope line
Slope_ea_y <- approxfun(x = c(min(profile_df$x), max(profile_df$x)),
                        y = c(min(profile_df$y), min(profile_df$y) + h ))

# Find where these lines intersect


# function that defines the different in y values for a given x
Calc_ydiff <- function(x){
  Profile_y(x) - Slope_ea_y(x)
  
}

x_int <- uniroot(Calc_ydiff, interval = c(profile_df$x[2], max(profile_df$x)))$root
y_int <- Slope_ea_y(x_int)

# Can now define the two polygons
# A1 and A2 in Figure 1
# or X and Y in Figure 2


x_1 = c(profile_df$x[profile_df$x < x_int], x_int, min(profile_df$x))
y_1 = c(profile_df$y[profile_df$x < x_int], y_int, min(profile_df$y))

poly_1_df = data_frame(x = x_1, y = y_1)


x_2 = c(x_int, profile_df$x[profile_df$x > x_int], max(profile_df$x), x_int)
y_2 = c(y_int, profile_df$y[profile_df$x > x_int], min(profile_df$y) + h, y_int)

poly_2_df = data_frame(x = x_2, y = y_2)

plot_ea +
  annotate(geom = "point", x = C_x, y = C_y, size = 2, colour = 'blue') +
  annotate(geom = "segment", x = min(profile_df$x), xend = max(profile_df$x), y = min(profile_df$y), yend = C_y, linetype = 'dashed') +
  geom_polygon(data = poly_1_df, aes(x,y), colour = 'blue', alpha = 0.2, fill = 'blue') +
  geom_polygon(data = poly_2_df, aes(x,y), colour = 'green', alpha = 0.2, fill = 'green') 


# Calculate areas of the two polygons

A1 <- pracma::polyarea(x_1[-length(x_1)], y_1[-length(y_1)])
A2 <- pracma::polyarea(rev(x_2[-length(x_2)]), rev(y_2[-length(y_2)]))

all.equal(A1, A2)


test_slope_ea = h/L
test_slope_ea

Calc_slope_ea(profile_df$x, profile_df$y)

all.equal(test_slope_ea/1000, Calc_slope_ea(profile_df$x, profile_df$y))