In lvaudor/woody: Use wood occurrence data and discharge data for wood flux prediction

library(woody)
library(tidyverse)
library(DiagrammeR)

Global process

grViz("
digraph boxes_and_circles {
  graph [overlap = true, fontsize = 10]
  node [shape = box, style= filled]
  wpath [fillcolor = burlywood1]
  qpath  [fillcolor = cadetblue1]
  W [fillcolor = burlywood1]
  Q  [fillcolor = cadetblue1]
  Wc [fillcolor = darkkhaki]
  Qc  [fillcolor = cadetblue1]
  Q [fillcolor = cadetblue1]
  Wwt [fillcolor = darkkhaki]
  Wwt_test [fillcolor = darkkhaki]
  myrf [fillcolor = gold]
  Wp [fillcolor = darkolivegreen2]

  node [shape = ellipse,fixedsize = false,color = grey, style= filled, fillcolor=grey90]
  'import_Wdata()'
  'import_Qdata()'
  'complete_Qdata()' 
  'complete_Wdata_with_Qdata()'
  'Wdata_as_waiting_times()'
  'run_rf()'
  'predict_rf()'

  # several 'edge' statements
  wpath -> 'import_Wdata()' [arrowhead = none]
  'import_Wdata()'-> W 
  qpath -> 'import_Qdata()' [arrowhead = none]
  'import_Qdata()'-> Q 
  Q -> 'complete_Qdata()' [arrowhead=none]
  'complete_Qdata()' -> Qc 
  W -> 'complete_Wdata_with_Qdata()' [arrowhead=none]
  Qc-> 'complete_Wdata_with_Qdata()' [arrowhead=none]
  'complete_Wdata_with_Qdata()' -> Wc
  Wc -> 'Wdata_as_waiting_times()' [arrowhead=none]
  'Wdata_as_waiting_times()' -> Wwt
  Wwt -> 'run_rf()' [arrowhead = none]
  'run_rf()' -> myrf 
  myrf -> 'predict_rf()' [arrowhead = none]
  'Wwt_test' -> 'predict_rf()' [arrowhead = none]
  'predict_rf()' -> 'Wp' 
}
")

Prepare data

Get W data about wood occurrences

grViz("
digraph boxes_and_circles {
  graph [overlap = true, fontsize = 10]
  node [shape = box, style= filled]
  wpath [fillcolor = burlywood]
  W [fillcolor = burlywood]

  node [shape = ellipse,fixedsize = false,color = grey, style= filled, fillcolor=grey90]
  'import_Wdata()'


  # several 'edge' statements
  wpath -> 'import_Wdata()' [arrowhead = none]
  'import_Wdata()'-> W 
}
")

wpath="../data-raw/mem_data/woody_data_train/"
mystation="V2942010"
fs::dir_tree(wpath)

library(woody)

Import all wood log files inside a file (1 directory per event)

W=import_Wdata(wpath,site="training") 

Get Q data about discharge

grViz("
digraph boxes_and_circles {
  graph [overlap = true, fontsize = 10]
  graph [overlap = true, fontsize = 10]
  node [shape = box, style= filled]
  qpath  [fillcolor = cadetblue1]
  Q  [fillcolor = cadetblue1]


  node [shape = ellipse,fixedsize = false,color = grey, style= filled, fillcolor=grey90]
  'import_Qdata()'

  qpath -> 'import_Qdata()' [arrowhead = none]
  'import_Qdata()'-> Q 

}
")

Import Q:

qpath="../data-raw/Qdata/Qdatc_Ain.csv"
Q=import_Qdata(path=qpath, site="Ain")

The first lines of Q look like this

knitr::kable(head(Q))

Complete Q into Qc

grViz("
digraph boxes_and_circles {
   graph [overlap = true, fontsize = 10]
  node [shape = box, style= filled]
  Qc  [fillcolor = cadetblue1]
  Q [fillcolor = cadetblue1]


  node [shape = ellipse,fixedsize = false,color = grey, style= filled, fillcolor=grey90]
  'complete_Qdata()' 


  # several 'edge' statements
  Q -> 'complete_Qdata()' [arrowhead=none]
  'complete_Qdata()' -> Qc 

}
")

Now complete Q with variables T_Q and S

result_file="../data-raw/results/Q.RDS"
if(!file.exists(result_file)){
  Qc=complete_Qdata(qtvar=Q,site="Ain")
  saveRDS(Qc,result_file)
}
Qc=readRDS(result_file)

The first lines of Qc look like this:

knitr::kable(head(Qc))

Complete W with Q and get Wc

grViz("
digraph boxes_and_circles {
 graph [overlap = true, fontsize = 10]
  node [shape = box, style= filled]
  W [fillcolor = burlywood]
  Wc [fillcolor = darkkhaki]
  Qc  [fillcolor = cadetblue1]

  node [shape = ellipse,fixedsize = false,color = grey, style= filled, fillcolor=grey90]
  'complete_Wdata_with_Qdata()'


  # several 'edge' statements
  W -> 'complete_Wdata_with_Qdata()' [arrowhead=none]
  Qc-> 'complete_Wdata_with_Qdata()' [arrowhead=none]
  'complete_Wdata_with_Qdata()' -> Wc
}
")

Complete W with discharge data (interpolating Q)

Wc=complete_Wdata_with_Qdata(Wdata=W,
                             Qdata=Qc)

The first lines of Wc look like this:

knitr::kable(head(Wc))

Transform Wc into Wwt

grViz("
digraph boxes_and_circles {
  graph [overlap = true, fontsize = 10]
  node [shape = box, style= filled]
  Wc [fillcolor = darkkhaki]
  Wwt [fillcolor = darkkhaki]

  node [shape = ellipse,fixedsize = false,color = grey, style= filled, fillcolor=grey90]
  'Wdata_as_waiting_times()'

  # several 'edge' statements

  Wc -> 'Wdata_as_waiting_times()' [arrowhead=none]
  'Wdata_as_waiting_times()' -> Wwt
}
")

Then calculate Wwt, where 1 row= 1 waiting time between two wood occurrences.

Wwt=Wdata_as_waiting_times(Wc)

The first lines of Wwt look like this:

knitr::kable(head(Wwt))

Run random forest

grViz("
digraph boxes_and_circles {
   graph [overlap = true, fontsize = 10]
  node [shape = box, style= filled]
  Wwt [fillcolor = darkkhaki]
  myrf [fillcolor = gold]

  node [shape = ellipse,fixedsize = false,color = grey, style= filled, fillcolor=grey90]

  'run_rf()'

  # several 'edge' statements
  Wwt -> 'run_rf()' [arrowhead = none]
  'run_rf()' -> myrf 
}
")

myrf <- run_rf(Wwt)
myrf

Get predictions and assess random forest performance

grViz("
digraph boxes_and_circles {
   graph [overlap = true, fontsize = 10]
  node [shape = box, style= filled]
  Wwt_test [fillcolor = darkkhaki]
  myrf [fillcolor = gold]
  Wp [fillcolor = darkolivegreen2]

  node [shape = ellipse,fixedsize = false,color = grey, style= filled, fillcolor=grey90]
  'predict_rf()'

  # several 'edge' statements

  Wwt_test -> 'predict_rf()' [arrowhead = none] 
  myrf -> 'predict_rf()' [arrowhead = none]
  'predict_rf()' -> 'Wp' 
}
")

On the training dataset itself

Predict on training data itself

Wp=predict_rf(obj_rf=myrf,
              newdata=Wwt)

R2:

calc_rf_R2(Wp)

On another dataset

wpathpred="../data-raw/mem_data/woody_data_predict/"

Wwt_test=import_Wdata(wpathpred) %>%
  complete_Wdata_with_Qdata(Qc) %>% 
  Wdata_as_waiting_times() 

Wp_test=predict_rf(Wwt_test,myrf)
calc_rf_R2(Wp_test)

lvaudor/woody documentation built on March 22, 2022, 9:53 a.m.