PosNegDifferenceLinePlot.R

# Find area difference of two lines and group by whether one line is higher or
# lower than the other. Possiblt a long way round but it works!

#--------- Load packages -------------
library(wk)
library(sf)
library(tidyverse)

#------------- Create some fake data in long form ---------------
fake_df <- tibble(scenario = c(rep('high',6),rep('low',6)),
                  time = c(c(1:6), c(1:6)),
                  Q = c(0.9, 5, 6, 4, 1, 1, 1, 4, 5, 4.1, 2, 2)) 

# This just repeats the fake data over and over. test on longer dataset...
extend <- function(x, n){
  do.call("rbind", replicate(n, x, simplify = FALSE)) %>%
    mutate(time= row_number())
}

fake_df <- fake_df %>% group_by(scenario) %>% group_split() %>%
  lapply(., extend,n=6) %>% # how many times to repeat
  bind_rows()

# #------------- Convert data to wide form -----------------------
t <- fake_df %>% pivot_wider(values_from=Q, names_from =scenario) %>%
  mutate(.min = ifelse(high<low, high, low),
         .max = ifelse(high>low, high, low),
         .lab = ifelse(high>low, 'attenuation', 'storage')) 

#------------- Create Polygon from the two lines ---------------------
poly <- (wk_polygon(xy(c(t$time,rev(t$time)), c(t$low, rev(t$high)))))%>%
  st_as_sfc() %>%
  st_make_valid()%>%
  st_cast('MULTIPOLYGON') %>%
  st_cast('POLYGON') %>%
  st_as_sf() %>%
  mutate(diffID = as.factor(row_number()),
         shp_area = st_area(.)) 
  
#------------- Get polygon ID number based on .lab value ------------------
polyID <- t  %>%
  mutate(diffID = as.factor(cumsum(.lab != lag(.lab, def = first(.lab)))+1)) %>%
  select(.lab, diffID) %>%
  distinct() 

#------------ Join polygon with id table ---------------------
polyJoin <- poly %>%
  left_join(.,polyID, by='diffID')

#------------ Plot ---------------------------------------
ggplot(polyJoin) +
  geom_sf( aes(fill=.lab), alpha=.5, colour=NA) +
  geom_line(data=fake_df, mapping=aes(x=time, y=Q, linetype=scenario), lwd=0.5 ) +
  scale_fill_brewer(palette='Set1') +
  scale_linetype_manual(values=c(1,2)) +
  theme_bw()

#-------------Get stats - ------------------------------
# for now just the sum but could get max mean attenuation/storage events 
# if you liked...

#TOTAL VOLUME of stored water
sum(polyJoin$shp_area[polyJoin$.lab=='storage'])
#TOTAL VOLUME of attenuated water
sum(polyJoin$shp_area[polyJoin$.lab=='attenuation'])

See also https://www.nsgrantham.com/fill-between-two-lines-ggplot2