Example application

The goal of our approach helps users to reproduce elements of a research workflow. As a demonstrator, this app allows users to reproduce elements (e.g. figures) of our paper ‘Rainfall and weather type effects on atmospheric deposition chemistry in the British Isles (1993 – 2015) and the potential influence of climate change’, which is currently under review.

Setup

I have preloaded the following datasets, packages, and scripts. Note that only the setup chunk in the Rmd document can be called by other code chunks by default. Future releases of learnr should allow better linking of outputs between chunks.

This code chunk reads all the required datasets and define functions to be used, which includes:

• Reading ECN rainfall chemistry data
• Reading ECN weather data
• Reading ERA5 weather reanalysis data
• Reading daily Lamb weather type data
• Reading daily Atlantic multidecadal oscillation (AMO) data
• Various functions to read, aggregate, and transform data in a usable form

# https://beta.rstudioconnect.com/content/2671/Combining-Shiny-R-Markdown.html
library(learnr)
tutorial_options(exercise.timelimit = 60, 
                 exercise.completion=TRUE,
                 exercise.eval=TRUE)
knitr::opts_chunk$set(
  tidy = FALSE
)

library(ncdf4)
library(lubridate)
library(plotly)
library(ggplot2)
library(lubridate)
library(stringr)
library(cowplot)
library(dplyr)
library(changepoint)
library(feather)
library(stringi)
library(stringr)
library(readr)
library(tidyr)
library(maps)
library(here)
library(ggrepel)

# color blind friendly palette (with grey):
cbf_1 <- c("#999999", "#E69F00", "#56B4E9", "#009E73", 
           "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# color_blind_friendly palette (with black): 
cbf_2 <- c("#000000", "#E69F00", "#56B4E9", "#009E73", 
           "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

era5 = function(ncfname){
  era5_single = function(ncfname){
    ncin <- nc_open(ncfname)
    #print(ncin)
    out = tibble(t  = make_datetime(1900,1,1,ncvar_get(ncin,"time"),0,0), 
                 cp = ncvar_get(ncin,"cp"),
                 lsp = ncvar_get(ncin,"lsp"),
                 lspf = ncvar_get(ncin,"lspf"),
                 ERA5_rain = ncvar_get(ncin,"tp"),)
    nc_close(ncin)
    return(out)
  }
  out = tibble()
  for (i in 1:length(ncfname)){
    out = rbind(out,era5_single(ncfname[i]))
  }
  
  return(out)
}

get_datechem <- function(site){
  lamb_type_lookup <-  read_csv(here('data',"lamb-type-lookup.csv"))  
  
  lamb_data = read_csv(here('data','20CR_1871-1947_ncep_1948-2020_12hrs_UK.csv')) %>% 
    dplyr::mutate(DATE = make_date(year,month,day)) %>% 
    dplyr::select(DATE,LWT) %>% 
    dplyr::mutate(LWT = as.character(LWT)) %>% 
    dplyr::mutate(LWT = setNames(lamb_type_lookup$Type,lamb_type_lookup$Num)[LWT])
  
  
  ### get ERA5
  
  # files = list.files(pattern = paste0('download_',site,"_"),path = '/data/rain_pattern/ERA5',full.names = T)
  # ERA5 = era5(files) %>% 
  #   mutate(cp = 1000*cp, lsp = 1000*lsp,ERA5_rain =cp+lsp, frac = cp/ERA5_rain) %>%  # m to mm
  #   group_by(date(t)) %>% 
  #   mutate(frac_daily = mean(frac)) %>% 
  #   ungroup() %>% 
  #   mutate(frac = ifelse(!between(frac,0,1),NA,frac) ) %>% 
  #   mutate(frac = ifelse(frac <1e-6,0,frac) ) %>% 
  #   tidyr::fill(frac,.direction="down")
  
  ### get rain flux
  rain_daily = read_feather("/data/rain_pattern/rain_hourly.feather") %>% 
    dplyr::filter(SITECODE == site) %>% 
    dplyr::mutate(SDATE = as.Date(SDATE,format="%d-%b-%y")) %>% 
    dplyr::rename(DATE = SDATE) %>%
    group_by(DATE) %>% 
    summarise(VALUE = sum(VALUE))
  
  ### end get rain flux
  
  # list of PC dates that passes QC
  QC_date = readr::read_csv('/data/rain_pattern/ECN_PC_chemistry-QC.csv') %>% 
    dplyr::filter(OVERALL_CHECK_OK == 1, SITECODE== site) %>% 
    dplyr::rename(DATE = SDATE) %>% 
    dplyr::filter(SITECODE == site) %>% 
    dplyr::mutate(DATE = as.Date(DATE,format="%d-%b-%y")) %>% 
    dplyr::select(DATE) %>% pull
  
  
  rain_chem = read_csv("/data/rain_pattern/ECN_PC1.csv")   %>% 
    dplyr::rename(DATE = SDATE) %>% 
    dplyr::filter(SITECODE == site) %>% 
    dplyr::filter(!FIELDNAME %in% c("VACUUM","PHAQCS","PHAQCU")) %>%  ## remove variables
    dplyr::filter(!str_detect(FIELDNAME, "Q"))    %>%
    #dplyr::mutate(FIELDNAME = paste(FIELDNAME,"(Rain)")) %>% 
    dplyr::mutate(DATE = as.Date(DATE,format="%d-%b-%y"))
  
  # outlier removal, note 1 mg/L to 1ppm
  # http://inside.mines.edu/~epoeter/_GW/resultsNOV03.pdf
  rain_chem = rain_chem %>% 
    tidyr::spread(FIELDNAME,VALUE) %>% 
    #dplyr::filter(!PO4P > 1.0) %>% 
    #dplyr::filter(between( (SODIUM/23)/(CHLORIDE/35.5), 0.8,1.2)) %>% # seasalt check
    dplyr::mutate(`non-marine sulphate` = (SO4S*3*20.83-28.21*0.104*CHLORIDE)/1000 ) %>% # in meq/L 
    dplyr::mutate(#ALKY = ALKY/61, # as for Hco3, everything to meq/L
      H = 10^(6-PH)/1000, # ph to meq/L
      ALUMINIUM=ALUMINIUM/26.98*3,
      CALCIUM=CALCIUM/40.08*2,
      CHLORIDE=CHLORIDE/35.45,
      IRON=IRON/55.845*2,
      MAGNESIUM=MAGNESIUM/24.3*2,
      POTASSIUM=POTASSIUM/39.1,
      SODIUM=SODIUM/22.99,
      NH4N=NH4N/14.01,
      NO3N=NO3N/14.01,
      PO4P=PO4P/15,
      SO4S=SO4S*3/96*2,
      EC=(H*349.8+SO4S*79.8+NO3N*71.4+NH4N*73.5+SODIUM*50.1+
            MAGNESIUM*53.05+CALCIUM*59.5+CHLORIDE*76.4+POTASSIUM*73.5), # tehorectical EC based on diffusion coefficient (Nernst-Einstein Equation),https://www.aqion.de/site/77
      EC_ratio = CONDY/EC,
      EC_balance = (CONDY-EC)/(CONDY+EC)
    ) %>%  # convert all ions to meq/L
    dplyr::mutate(sum_anion = CHLORIDE+NO3N+PO4P+SO4S, 
                  sum_cation= NH4N+CALCIUM+MAGNESIUM+POTASSIUM+SODIUM,
                  charge_balance = (sum_cation-sum_anion)/(sum_cation+sum_anion)) %>% 
    #dplyr::filter(!abs(charge_balance) > 0.05) %>% 
    #dplyr::filter(!abs(EC_balance) > 0.2) %>%  # measured TDS/ theorectical TDS should be between 1.0 and 1.2 
    tidyr::gather(FIELDNAME,VALUE,-SITECODE,-DATE) %>% 
    dplyr::filter(DATE %in% QC_date)
  
 
  # convert all ions to meq/L
  
  date_chem = distinct(rain_chem,DATE) %>% 
    dplyr::arrange(DATE) %>% 
    dplyr::mutate(SDATE = as.Date(DATE,format="%d-%b-%y")) %>% 
    dplyr::mutate(DIFF = as.numeric(difftime(DATE,lag(DATE,1)))) %>% 
    dplyr::mutate(STARTDATE = lag(DATE,1)+1, ENDDATE = DATE) %>% 
    dplyr::filter(!is.na(DIFF)) %>%
    dplyr::filter(between(DIFF,6,8)) #%>% 
  #dplyr::mutate(lamb_str2 = NA,mean_rain_per_day = NA) # maximum len of continuous rain > 0.1
  
  
  get_lamb_str2 <- function(lamb_data,STARTDATE,ENDDATE){
    lamb_data %>% 
      filter(between(DATE,STARTDATE,ENDDATE)) %>% 
      select(LWT) %>% 
      mutate(LWT = case_when(LWT == "C" | LWT == "W" |
                               LWT == "SW" |LWT == "CW" |
                               LWT == "CSW"  ~ "1" , TRUE ~ "0")) %>% #event generating types
      pull() %>% paste(collapse = '')
  }
  get_lamb_strW <- function(lamb_data,STARTDATE,ENDDATE){
    lamb_data %>% 
      filter(between(DATE,STARTDATE,ENDDATE)) %>% 
      select(LWT) %>% 
      mutate(LWT = case_when(LWT == "W" |LWT == "SW" |LWT == "CW" |LWT == "CSW"  ~ "1" ,
                             TRUE ~ "0")) %>% #event generating types
      pull() %>% paste(collapse = '')
  } 
  get_lamb_strC <- function(lamb_data,STARTDATE,ENDDATE){
    lamb_data %>% 
      filter(between(DATE,STARTDATE,ENDDATE)) %>% 
      select(LWT) %>% 
      mutate(LWT = case_when(LWT=="C"|LWT=="CSW" |LWT=="CW"|LWT=="CNW"|LWT=="CN"|
                               LWT=="CNE"|LWT=="CE"|LWT=="CSE"|LWT=="CS"  ~ "1" ,
                             TRUE ~ "0")) %>% #event generating types
      pull() %>% paste(collapse = '')
  }
  
  date_chem = date_chem %>% 
    rowwise() %>%
    mutate(lamb_str2= get_lamb_str2(get("lamb_data"),STARTDATE,ENDDATE),
           lamb_strW= get_lamb_strW(get("lamb_data"),STARTDATE,ENDDATE),
           lamb_strC= get_lamb_strC(get("lamb_data"),STARTDATE,ENDDATE),
           mean_rain_per_day = filter(rain_daily,between(DATE,STARTDATE,ENDDATE)) %>% 
             select(VALUE) %>% pull() %>% mean(na.rm=T),
           #cp_week=filter(ERA5,between(`date(t)`,STARTDATE,ENDDATE)) %>% select(cp) %>% sum(),
           #era5_week = filter(ERA5,between(`date(t)`,STARTDATE,ENDDATE)) %>% select(ERA5_rain) %>% sum() 
    ) %>% 
    ungroup()
  
  date_chem = date_chem %>% mutate(SITECODE=site) %>% select(-SDATE,-STARTDATE,-ENDDATE) %>% 
    left_join(rain_chem %>% select(-SITECODE) %>% tidyr::spread(FIELDNAME,VALUE), by="DATE") %>%
    gather(FIELDNAME,VALUE,-DATE,-DIFF,-lamb_str2,-lamb_strW,-lamb_strC,-SITECODE,-mean_rain_per_day,
           #-cp_week,-era5_week,
           -VOLUME)
  
  print(paste0('Complete ',site,'!'))
  
  
  return(date_chem)  
}


data0 = rbind(get_datechem('T02'),get_datechem('T04'),get_datechem('T06'),get_datechem('T11'))

## [1] "Complete T02!"
## [1] "Complete T04!"
## [1] "Complete T06!"
## [1] "Complete T11!"

date_chem = data0 %>%
   tidyr::spread(FIELDNAME,VALUE)  %>% 
  mutate(SITECODE = case_when(
    SITECODE == 'T02' ~ 'Glensaugh',
    SITECODE == 'T04' ~ 'Moor_House',
    SITECODE == 'T06' ~ 'Rothamsted',
    SITECODE == 'T01' ~ 'Drayton',
    SITECODE == 'T07' ~ 'Sourthrope',
    SITECODE == 'T11' ~ 'Snowdon',
    TRUE ~ SITECODE
  )) %>% 
mutate(lamb2count=sapply(strsplit(lamb_str2, ''), function(x) sum(as.numeric(x))) ,
       lambCcount=sapply(strsplit(lamb_strC, ''), function(x) sum(as.numeric(x))) ,
       lambWcount=sapply(strsplit(lamb_strW, ''), function(x) sum(as.numeric(x))) ,
       yr_grp = case_when(
         between(year(DATE),1992,1995) ~ '1992-1995',
         between(year(DATE),1996,1999) ~ '1996-1999',
         between(year(DATE),2000,2003) ~ '2000-2003',
         between(year(DATE),2004,2007) ~ '2004-2007',
         between(year(DATE),2008,2011) ~ '2008-2011',
         between(year(DATE),2012,2015) ~ '2012-2015',
         TRUE                      ~ NA_character_),
       lamb_grp2 = case_when(
         between(lamb2count,0,2) ~ '0-2',
         #between(lamb2count,2,3) ~ '2-3',
         #         between(lamb2count,3,3) ~ '3',
         #         between(lamb2count,4,4) ~ '4',
         between(lamb2count,3,5) ~ '3-5',
         between(lamb2count,6,8) ~ '6 or more',
         TRUE                      ~ NA_character_),
       lambCcountgrp = case_when(
         between(lambCcount,0,1) ~ '0-1',
         between(lambCcount,2,3) ~ '2-3',
         between(lambCcount,4,5) ~ '4-5',
         between(lambCcount,6,8) ~ '6 or more',
         TRUE                      ~ NA_character_),
       lambWcountgrp = case_when(
         between(lambWcount,0,1) ~ '0-1',
         between(lambWcount,2,3) ~ '2-3',
         between(lambWcount,4,5) ~ '4-5',
         between(lambWcount,6,8) ~ '6 or more',
         TRUE                      ~ NA_character_),
       # conv3grp = case_when(
       #   between(cp_week/era5_week,0,0.2) ~ '0-20%',
       #   between(cp_week/era5_week,0.2,0.5) ~ '20-50%',
       #   between(cp_week/era5_week,0.5,0.8) ~ '50-80%',
       #   between(cp_week/era5_week,0.8,1.0) ~ '80-100%',
       #   TRUE                      ~ NA_character_),
       season = case_when(
         month(DATE) %in% 3:5 ~ 'Spring (MAM)',
         month(DATE) %in% 6:8 ~ 'Summer (JJA)',
         month(DATE) %in% 9:11 ~ 'Autumn (SON)',
         month(DATE) %in% c(12,1,2) ~ 'Winter (DJF)',
         TRUE                      ~ NA_character_),
       season = factor(season,levels=c('Spring (MAM)','Summer (JJA)','Autumn (SON)','Winter (DJF)')),
       #wm = weighted.mean(VALUE,VOLUME)
) #%>% 
  #ungroup() %>% 
  # filter(!(FIELDNAME == 'CONDY' & VALUE > 100)) %>%
  # filter(!(FIELDNAME == 'NO3N' & VALUE > 5)) %>%
  # filter(!(FIELDNAME == 'NH4N' & VALUE > 5)) %>%
  # filter(!(FIELDNAME == 'CHLORIDE' & VALUE > 7.5)) %>%
  # filter(!(FIELDNAME == 'SODIUM' & VALUE > 5)) %>%
  # filter(!(FIELDNAME == 'non-marine sulphate' & VALUE > 50)) %>%
  #filter(!(VALUE < 0))

date_chem = date_chem %>% 
  mutate(SODIUM=SODIUM*1000, CHLORIDE=CHLORIDE*1000, NH4N=NH4N*1000,
         NO3N=NO3N*1000, `non-marine sulphate`=`non-marine sulphate`*1000)
# date_chem = date_chem %>% group_by(SITECODE) %>% 
#   mutate(conv3quant = rank(cp_week/era5_week)/length(cp_week/era5_week)) %>% 
#   mutate( conv3grp = case_when(
#     between(conv3quant,0,0.33) ~ '0-33%',
#     between(conv3quant,0.33,0.67) ~ '33-67%',
#     between(conv3quant,0.67,1.0) ~ '67-100%',
#     TRUE                      ~ NA_character_))