Skip to content

Instantly share code, notes, and snippets.

@h-a-graham
Last active October 18, 2021 11:58
Show Gist options
  • Save h-a-graham/9076242027dfcab76d70868a372fab97 to your computer and use it in GitHub Desktop.
Save h-a-graham/9076242027dfcab76d70868a372fab97 to your computer and use it in GitHub Desktop.
# 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'])
@kent37
Copy link

kent37 commented Oct 18, 2021

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment