library(rsofun)

## Warning: replacing previous import 'dplyr::intersect' by
## 'lubridate::intersect' when loading 'rsofun'

## Warning: replacing previous import 'dplyr::union' by 'lubridate::union'
## when loading 'rsofun'

## Warning: replacing previous import 'dplyr::setdiff' by 'lubridate::setdiff'
## when loading 'rsofun'

## Warning: replacing previous import 'purrr::invoke' by 'rlang::invoke' when
## loading 'rsofun'

## Warning: replacing previous import 'purrr::flatten_int' by
## 'rlang::flatten_int' when loading 'rsofun'

## Warning: replacing previous import 'purrr::flatten_dbl' by
## 'rlang::flatten_dbl' when loading 'rsofun'

## Warning: replacing previous import 'purrr::list_along' by
## 'rlang::list_along' when loading 'rsofun'

## Warning: replacing previous import 'purrr::%||%' by 'rlang::%||%' when
## loading 'rsofun'

## Warning: replacing previous import 'purrr::flatten_chr' by
## 'rlang::flatten_chr' when loading 'rsofun'

## Warning: replacing previous import 'purrr::prepend' by 'rlang::prepend'
## when loading 'rsofun'

## Warning: replacing previous import 'purrr::flatten' by 'rlang::flatten'
## when loading 'rsofun'

## Warning: replacing previous import 'Metrics::ll' by 'rlang::ll' when
## loading 'rsofun'

## Warning: replacing previous import 'purrr::rep_along' by 'rlang::rep_along'
## when loading 'rsofun'

## Warning: replacing previous import 'purrr::splice' by 'rlang::splice' when
## loading 'rsofun'

## Warning: replacing previous import 'purrr::modify' by 'rlang::modify' when
## loading 'rsofun'

## Warning: replacing previous import 'purrr::flatten_lgl' by
## 'rlang::flatten_lgl' when loading 'rsofun'

## Warning: replacing previous import 'purrr::as_function' by
## 'rlang::as_function' when loading 'rsofun'

## Warning: replacing previous import 'purrr::%@%' by 'rlang::%@%' when
## loading 'rsofun'

library(dplyr)

## Warning: package 'dplyr' was built under R version 3.3.2

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(readr)
library(lubridate)

## Warning: package 'lubridate' was built under R version 3.3.2

## 
## Attaching package: 'lubridate'

## The following object is masked from 'package:base':
## 
##     date

source("get_obs_bysite_fluxnet2015.R")
source("train_predict_fvar.R")
source("predict_nn.R")
source("prepare_trainingdata_fvar.R")
source("remove_outliers.R")
source("get_consecutive.R")
source("add_alpha.R")
source("get_droughts_fvar.R")
source("align_events.R")

Get FLUXNET 2015 data

Read observational data including soil moisture for one site ("AU-How") from a FLUXNET 2015 data file (standard).

getvars <- c( 
  "GPP_NT_VUT_REF", 
  "GPP_DT_VUT_REF",                
  "LE_F_MDS", 
  "LE_F_MDS_QC", 
  "NEE_VUT_REF_NIGHT_QC", 
  "NEE_VUT_REF_DAY_QC", # quality flag that goes with GPP 
  "TA_F",
  "VPD_F",
  "P_F",
  "SW_IN_F",
  "NETRAD"
  )
df <- get_obs_bysite_fluxnet2015( sitename="AU-How", path_fluxnet2015="./", timescale="d", getvars=getvars )

## Warning: package 'bindrcpp' was built under R version 3.3.2

## Warning in as.POSIXlt.POSIXct(x, tz): unknown timezone 'zone/tz/2018g.1.0/
## zoneinfo/Europe/Madrid'

## [1] "Converting: soilm_obs_mean = mean across different soil depths (SWC_F_MDS), with na.rm = TRUE"
## [1] "Converting: temp = TA_F"
## [1] "Converting: vpd = VPD_F * 1e2 (given in hPa, required in Pa)"
## [1] "Converting: prec = P_F"
## [1] "Converting: transp_obs = LE_F_MDS (given in W m-2, converting to MW m-2)"
## [1] "Converting: swin = SW_IN_F * 60 * 60 * 24 (given in W m-2, required in J m-2 d-1)"
## [1] "Converting: swin = SW_IN_F * 60 * 60 * 24 (given in W m-2, required in J m-2 d-1)"
## [1] "Converting: ppfd = swin * kfFEC * 1.0e-6 (convert from J/m2/d to mol/m2/d; kfFEC = 2.04 is the flux-to-energy conversion, micro-mol/J (Meek et al., 1984))"
## [1] "Converting: netrad = NETRAD * 60 * 60 * 24 (given in W m-2 (avg.), required in J m-2 (daily total))"
## [1] "Converting: gpp_obs = mean( GPP_NT_VUT_REF, GPP_DT_VUT_REF )"

## get some mo' (modelled soil moisture)
load("/alphadata01/bstocker/nn_fluxnet2015/data/modobs_fluxnet2015_s11_s12_s13_with_SWC_v3.Rdata")
df_wcont_splash  <- fluxnet[[ "AU-How" ]]$ddf$s11 %>% 
                      as_tibble() %>% 
                      mutate( date=ymd( paste0( as.character(year), "-01-01") ) + days(doy-1), wcont_splash=wcont )
df_wcont_splash %>% write_csv( path="inst/ext/df_wcont_splash.csv")

df_wcont_swbm    <- fluxnet[[ "AU-How" ]]$ddf$s12 %>% 
                      as_tibble() %>% 
                      mutate( date=ymd( paste0( as.character(year), "-01-01") ) + days(doy-1), wcont_swbm=wcont )
df_wcont_swbm %>% write_csv( path="inst/ext/df_wcont_swbm.csv")

df_wcont_et      <- fluxnet[[ "AU-How" ]]$ddf$swc_by_etobs %>% 
                      as_tibble() %>% 
                      mutate( date=ymd( paste0(as.character(year), "-01-01") ) + months(moy-1) + days(dom-1), wcont_et=soilm_from_et )
df_wcont_et %>% write_csv( path="inst/ext/df_wcont_et.csv")

df_wcont_et_orth <- fluxnet[[ "AU-How" ]]$ddf$swc_by_etobs %>% 
                      as_tibble() %>% 
                      mutate( date=ymd( paste0(as.character(year), "-01-01") ) + months(moy-1) + days(dom-1), wcont_et_orth=soilm_from_et_orthbucket )
df_wcont_et_orth %>% write_csv( path="inst/ext/df_wcont_et_orth.csv")

Prepare data for training

Prepare training data, removing NAs and outliers.

## use all observational soil moisture data
varnams_soilm <- df %>% 
  dplyr::select( starts_with("SWC_") ) %>% 
  dplyr::select( -ends_with("QC") ) %>% 
  names()

## define settings used by multiple functions as a list
settings <- list( 
  target        = "transp_obs", 
  predictors    = c("temp","vpd","swin"), 
  varnams_soilm = varnams_soilm 
  )

df_train <- prepare_trainingdata_fvar( df, settings )

Get soil moisture threshold

profile_soilmthreshold_fvar( df_train, settings )

Train models and predict fvar

Train models for one set of soil moisture input data. In this case it’s observational data.

df_nn_soilm_obs <- train_predict_fvar( 
  df_train,
  settings,
  hidden_good        = 10, 
  hidden_all         = 10, 
  soilm_threshold    = 0.6, 
  nrep               = 3, 
  weights            = NA, 
  package            = "nnet" 
  )

## [1] "NN repetition 1"

## Loading required package: nnet

## Loading required package: caret

## Warning: package 'caret' was built under R version 3.3.2

## Loading required package: lattice

## Loading required package: ggplot2

## Warning: package 'ggplot2' was built under R version 3.3.2

## [1] "NN repetition 2"
## [1] "NN repetition 3"

Do the training again using another set of soil moisture data from a model.

df <- df %>% left_join( select(df_wcont_swbm, date, wcont_swbm), by="date" )
settings$varnams_soilm <- "wcont_swbm"

df_train <- prepare_trainingdata_fvar( df, settings )

df_nn_soilm_swbm <- train_predict_fvar( 
  df_train,
  settings,
  hidden_good        = 10, 
  hidden_all         = 10, 
  soilm_threshold    = 0.6, 
  nrep               = 3, 
  weights            = NA, 
  package            = "nnet" 
  )

## [1] "NN repetition 1"
## [1] "NN repetition 2"
## [1] "NN repetition 3"

Identify soil moisture droughts

## first aggregate fvar derived from different soil moisture datasets
df_nn <- bind_rows( df_nn_soilm_obs, df_nn_soilm_swbm ) %>% 
  group_by( date ) %>%
  summarise( 
    fvar       = mean( fvar, na.rm=TRUE ), 
    fvar_min   = min(  fvar, na.rm=TRUE ),
    fvar_max   = max(  fvar, na.rm=TRUE ), 
    fvar_med   = median(  fvar, na.rm=TRUE ), 
    var_nn_pot = mean( var_nn_pot, na.rm=TRUE ), 
    var_nn_act = mean( var_nn_act, na.rm=TRUE )
    ) %>%
  left_join( select_( df, "date", settings$target ), by = "date" )

out_droughts <- get_droughts_fvar( 
  df_nn, 
  nam_target     = settings$target, 
  leng_threshold = 10, 
  df_soilm       = select(df_wcont_swbm, date, soilm=wcont_swbm),
  df_par         = select(df, date, par=ppfd),
  par_runmed     = 10
  )

## Warning in runmed(df_nn$fvar_filled[idxs], par_runmed): 'k' must be odd!
## Changing 'k' to 11

## [1] "get drought events ..."
## [1] "number of events before pruning: 20"
## [1] "number of events after trimming: 20"
## [1] "number of events after pruning by soil moisture: 20"

## Loading required package: hydroGOF

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 3.3.2

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## [1] "number of events after pruning by mod vs obs: 16"
## [1] "number of events after pruning by PAR level: 16"
## [1] "number of events after pruning: 16"

Plot what we got.

with( out_droughts$df_nn, plot(date, fvar_smooth_filled, type="l" ) )
rect( 
  out_droughts$df_nn$date[out_droughts$droughts$idx_start], rep( -99, nrow(out_droughts$droughts) ), 
  out_droughts$df_nn$date[(out_droughts$droughts$idx_start+out_droughts$droughts$len-1)], rep( 99, nrow(out_droughts$droughts) ), 
  col=rgb(0,0,0,0.2), border=NA )

Align data by droughts

dovars <- c("fvar", "soilm")
df_nn <- out_droughts$df_nn %>% mutate(site="AU-How")
df_alg <- align_events( 
  select( df_nn, site, date, one_of(dovars)), 
  select( df_nn, site, date, isevent=is_drought_byvar_recalc ),
  dovars,
  leng_threshold = 10,
  before=20,
  after=79,
  nbins=10,
  do_norm=FALSE,
  normbin=2 )

## Loading required package: tidyr

## Warning: package 'tidyr' was built under R version 3.3.2

Plot what we got

with( df_alg$df_dday_aggbydday, plot(    dday, fvar_median, type="l", ylim=c(0,1.2)) )
with( df_alg$df_dday_aggbydday, polygon( c(dday, rev(dday)), c(fvar_q33, rev(fvar_q66)), col=rgb(0,0,0,0.2), border = NA ) )
abline(h=1, lty=3)
abline(v=0, col='red')