The QuantiSlakeTest, dynamic weighting of soil under water to measure soil structural stability
Openscience & Reproducible Research of Vanwin & Hardy (2022)

Table of Contents

1. Introduction

This is the public repository linked to the paper "The QuantiSlakeTest, dynamic weighting of soil under water to measure soil structural stability" submitted to the SOIL journal. You will find open data, meta data and commented R-scripts.

[to be amended]

This README file is done in a reproducible research spirit based on the powerful combo of emacs, org and babel. All the code blocks of this file have been extracted in one file, init.R in the root of the repository. The raw data are stored in the data-raw folder of the repository.

2. Initialisation

2.1. Cleaning session and workspace

Geek zonerm(list=ls())
load(file= "./data-output/qst-openscience.rda")

2.2. Loading packages

Geek zonelibrary(dplyr)
library(ggplot2) 
require(ggpubr) ## for one graph with ggarrange
require(lme4) ## for models
require(gt) # for beautiful table
require(car)
require(emmeans)
library(conflicted)

2.3. Conflicts

Geek zone## ggplot2
conflict_prefer("annotate", "ggplot2")

2.4. Read data

2.4.1. QST data (1 line = 1 qst experiment)

Geek zoneqst.opendata <- read.csv(
  file = c(
    "./data-raw/qst-opendata.csv")
)
row.names(qst.opendata) <- qst.opendata$slake

qst.opendata$modalite_uk_abb <- factor(
  qst.opendata$modalite_uk_abb,
  levels = c(
    "P1 K0","P1 K1","P1 K2","P2 K0","P2 K1","P2 K2",
    "P","RT",
    "RE","FYM","RR"
  )
)

head(qst.opendata, n = 10) %>%
  mutate(across(where(is.double),
                \(x) signif(x, digits = 3)))
  slake parcelle.x sample.id series.id starting.date starting.time duration tmax Wt0 Wend Wmax.Wt0 Wt0.Wend Warchi Sl.0.max Sl.max.30 Sl.max.60 Sl.max.300 Sl.max.600 local.slope.max.30 local.slope.30.60 local.slope.60.300 local.slope.300.600 t25 t50 t75 t90 t95 t99 delta.t.max.25 delta.t.25.50 delta.t.50.75 delta.t.75.90 delta.t.90.95 delta.t.95.99 oldAUCfalse AUC AUCWendcurve AUCWendrect AUCWt0curve AUCWt0rect AUCWendsum AUCWt0sum AUC.100 AUC.200 AUC.300 AUC.400 AUC.500 AUC.600 AUC.700 AUC.800 AUC.900 AUC.1000 campagne essai.x culture subzone X modalite essaib parcellea parcelleb traitement.x bloc terre plateforme ligne colonne traitment2 traitement3 traitement4 newbloc essai.b parcelle.b parcelle.c sampleid corg humus pHeau sampleidgranulo clay siltf siltc silt sandf sandc sand ib nlisting essai.y parcelle.y traitement.y oc.clay tare tarepfrais tarepsec tarep105 pfrais psec p105 watercontent pctms watercontentcorr dapp rhob Wroot sample parcelle mode profondeur dryingtemp dryingj eau essaiuk modaliteuk modaliteukabb
382 382 39 ErnageA-39-PK-39-P2K0-2019-D-C Ernage A-39 2019-06-20 11:37:06 971 4.69 0.94 0.061 0.06 0.88 0.42 0.013 -0.017 -0.011 -0.0029 -0.0015 -0.017 -0.0052 -0.00083 -0.00013 21.4 32.5 67.2 146 360 785 21.4 11.1 34.7 78.7 214 425 126 126 72.9 53.5 0.463 126 126 126 47.9 63 75.7 86.7 96.4 106 114 121 126   Ernage A Ernage A Froment 39 5 P2 K0 PK PK39 E-39-2019 P2K0 2 Ernage A PK     P2 K0     PK 39 PK39 T193414 11.5 2.3 6.95 G190839 15 28 50.3 78.3 6.14 0.583 6.73 2.1 31 PK E-39-2019 P2 K0 0.077 8.96 164 136 137 155 127 128 18.2 0.986 18.5 1.27 1.25   PK-39-P2K0-2019-D-C 39 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P2 K0 P2 K0
383 383 39 ErnageA-39-PK-39-P2K0-2019-D-B Ernage A-39 2019-06-20 11:50:47 719 2.11 0.93 0.21 0.07 0.72 0.44 0.033 -0.014 -0.0091 -0.0026 -0.0013 -0.014 -0.0039 -0.00093 -8e-05 19.9 29.4 96.5 265 298 331 19.9 9.44 67.1 169 33 33 211 249 61.5 187 0.489 248 249 249 58.1 94.7 124 145 165 186 205       Ernage A Ernage A Froment 39 5 P2 K0 PK PK39 E-39-2019 P2K0 2 Ernage A PK     P2 K0     PK 39 PK39 T193414 11.5 2.3 6.95 G190839 15 28 50.3 78.3 6.14 0.583 6.73 2.1 31 PK E-39-2019 P2 K0 0.077 8.96 164 136 137 155 127 128 18.2 0.986 18.5 1.27 1.25   PK-39-P2K0-2019-D-B 39 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P2 K0 P2 K0
384 384 40 ErnageA-40-PK-40-P2K2-2019-D-B Ernage A-40 2019-06-20 12:04:49 972 25.2 0.9 0.25 0.1 0.65 0.44 0.004 -0.0005 -0.0008 -0.0015 -0.0012 -0.0005 -0.0011 -0.0016 -0.00087 162 296 383 558 725 911 162 134 87.2 174 167 186 449 449 229 220 7.56 442 449 449 96.7 180 247 298 334 367 398 424 449   Ernage A Ernage A Froment 40 6 P2 K2 PK PK40 E-40-2019 P2K2 2 Ernage A PK     P2 K2     PK 40 PK40 T193415 11.7 2.35 6.8 G190840 16.6 30.3 47.3 77.5 5.16 0.651 5.81 2.02 32 PK E-40-2019 P2 K2 0.0705 8.57 184 150 149 175 141 140 19.5 0.986 19.8 1.41 1.39   PK-40-P2K2-2019-D-B 40 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P2 K2 P2 K2
386 386 16 ErnageA-16-PK-16-P1K0-2019-D-C Ernage A-16 2019-06-20 13:24:07 971 12.4 0.93 0.64 0.07 0.29 0.43 0.0057 -0.0013 -0.0013 -0.0009 -0.00054 -0.0013 -0.0013 -0.0008 -0.00019 85.5 155 328 600 755 940 85.5 69.7 173 273 154 185 648 648 85.8 562 2.56 645 648 648 94.5 176 251 322 391 457 527 588 648   Ernage A Ernage A Froment 16 2 P1 K0 PK PK16 E-16-2019 P1K0 1 Ernage A PK     P1 K0     PK 16 PK16 T193411 11.9 2.38 6.86 G190836 19.9 27.5 47.4 74.9 5.06 0.214 5.27 1.76 28 PK E-16-2019 P1 K0 0.06 8.96 169 140 133 160 131 124 18.3 0.982 18.7 1.31 1.28   PK-16-P1K0-2019-D-C 16 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P1 K0 P1 K0
387 387 16 ErnageA-16-PK-16-P1K0-2019-D-B Ernage A-16 2019-06-20 13:29:54 1130 32.5 0.93 0.62 0.07 0.31 0.44 0.0022 -0.002 -0.0017 -0.00093 -0.00059 -0.002 -0.0014 -0.00074 -0.00024 86 195 349 569 726 881 86 109 154 220 157 155 642 642 96 546 3.16 639 642 642 95.1 176 252 323 392 457 525 584 642 700 Ernage A Ernage A Froment 16 2 P1 K0 PK PK16 E-16-2019 P1K0 1 Ernage A PK     P1 K0     PK 16 PK16 T193411 11.9 2.38 6.86 G190836 19.9 27.5 47.4 74.9 5.06 0.214 5.27 1.76 28 PK E-16-2019 P1 K0 0.06 8.96 169 140 133 160 131 124 18.3 0.982 18.7 1.31 1.28   PK-16-P1K0-2019-D-B 16 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P1 K0 P1 K0
388 388 41 ErnageA-41-PK-41-P1K2-2019-D-C Ernage A-41 2019-06-20 13:45:11 726 11.8 0.85 0.25 0.15 0.6 0.43 0.013 -0.0022 -0.0045 -0.0024 -0.0012 -0.0022 -0.0068 -0.0019 -8.3e-05 65.8 77.7 177 286 319 471 65.8 12 99.2 109 32.8 153 275 319 93.5 225 6.22 312 319 319 83.3 129 163 190 216 241 267       Ernage A Ernage A Froment 41 7 P1 K2 PK PK41 E-41-2019 P1K2 2 Ernage A PK     P1 K2     PK 41 PK41 T193416 12.5 2.5 6.85 G190841 16.9 31.9 46.2 78.1 4.38 0.661 5.04 1.97 33 PK E-41-2019 P1 K2 0.074 8.03 176 143 143 168 135 135 19.5 0.986 19.8 1.35 1.33   PK-41-P1K2-2019-D-C 41 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P1 K2 P1 K2
391 391 40 ErnageA-40-PK-40-P2K2-2019-D-D Ernage A-40 2019-06-24 11:46:49 1560 23.5 0.9 0.46 0.1 0.44 0.44 0.0043 -0.0024 -0.0033 -0.0015 -0.00082 -0.0024 -0.0042 -0.0011 -0.00013 72.6 91.5 252 535 941 1530 72.6 18.9 161 283 405 587 512 512 107 404 4.01 508 512 512 90.5 155 213 266 317 368 420 466 512 557 Ernage A Ernage A Froment 40 6 P2 K2 PK PK40 E-40-2019 P2K2 2 Ernage A PK     P2 K2     PK 40 PK40 T193415 11.7 2.35 6.8 G190840 16.6 30.3 47.3 77.5 5.16 0.651 5.81 2.02 32 PK E-40-2019 P2 K2 0.0705 8.57 184 150 149 175 141 140 19.5 0.986 19.8 1.41 1.39   PK-40-P2K2-2019-D-D 40 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P2 K2 P2 K2
392 392 40 ErnageA-40-PK-40-P2K2-2019-D-E Ernage A-40 2019-06-24 11:50:40 626 50.3 0.83 0.74 0.17 0.095 0.41 0.0034 -0.00091 -0.0012 -0.00082   -0.00091       110 150 275 297 440 473 110 39.8 125 22 143 33 508 710 46.2 664 19.4 691 710 710 97.2 184 264 339 413 486         Ernage A Ernage A Froment 40 6 P2 K2 PK PK40 E-40-2019 P2K2 2 Ernage A PK     P2 K2     PK 40 PK40 T193415 11.7 2.35 6.8 G190840 16.6 30.3 47.3 77.5 5.16 0.651 5.81 2.02 32 PK E-40-2019 P2 K2 0.0705 8.57 184 150 149 175 141 140 19.5 0.986 19.8 1.41 1.39   PK-40-P2K2-2019-D-E 40 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P2 K2 P2 K2
401 401 40 ErnageA-40-PK-40-P2K2-2019-D-F Ernage A-40 2019-06-24 13:16:25 1500 25.3 0.87 0.49 0.13 0.38 0.43 0.0051 -0.0006 -0.00096 -0.0011 -0.00076 -0.0006 -0.0013 -0.0011 -0.00044 185 284 447 661 816 1190 185 98.2 163 214 154 371 597 597 167 430 12.8 584 597 597 96.6 183 262 327 388 450 500 549 597 645 Ernage A Ernage A Froment 40 6 P2 K2 PK PK40 E-40-2019 P2K2 2 Ernage A PK     P2 K2     PK 40 PK40 T193415 11.7 2.35 6.8 G190840 16.6 30.3 47.3 77.5 5.16 0.651 5.81 2.02 32 PK E-40-2019 P2 K2 0.0705 8.57 184 150 149 175 141 140 19.5 0.986 19.8 1.41 1.39   PK-40-P2K2-2019-D-F 40 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P2 K2 P2 K2
403 403 39 ErnageA-39-PK-39-P2K0-2019-D-D Ernage A-39 2019-06-24 13:40:09 1180 14.7 0.89 0.47 0.11 0.42 0.43 0.0075 -0.0024 -0.0038 -0.0017 -0.00088 -0.0024 -0.0052 -0.0012 -7.8e-05 57.7 84.9 133 198 308 429 57.7 27.3 48.2 65.1 110 121 468 468 57.7 410 3.92 464 468 468 86.9 144 195 243 290 337 381 424 468 526 Ernage A Ernage A Froment 39 5 P2 K0 PK PK39 E-39-2019 P2K0 2 Ernage A PK     P2 K0     PK 39 PK39 T193414 11.5 2.3 6.95 G190839 15 28 50.3 78.3 6.14 0.583 6.73 2.1 31 PK E-39-2019 P2 K0 0.077 8.96 164 136 137 155 127 128 18.2 0.986 18.5 1.27 1.25   PK-39-P2K0-2019-D-D 39 K 3-8 Air libre 30 Eau déminéralisée P-K mineral fertiliser P2 K0 P2 K0

2.4.2. QST aggregated data (mean and sd by plot)

Geek zoneqst.mean.sd.opendata <- read.csv(
  file = c(
    "./data-raw/qst-mean-sd-opendata.csv"
  )
)

head(qst.mean.sd.opendata, n = 10) %>%
  mutate(across(where(is.double), \(x) signif(x, digits = 3)))
  series.id parcelle.c bloc campagne modaliteukabb essaiuk slakeMean slakeSD durationMean durationSD tmaxMean tmaxSD Wt0Mean Wt0SD WendMean WendSD Wmax.Wt0Mean Wmax.Wt0SD Wt0.WendMean Wt0.WendSD WarchiMean WarchiSD Sl.0.maxMean Sl.0.maxSD Sl.max.30Mean Sl.max.30SD Sl.max.60Mean Sl.max.60SD Sl.max.300Mean Sl.max.300SD Sl.max.600Mean Sl.max.600SD local.slope.max.30Mean local.slope.max.30SD local.slope.30.60Mean local.slope.30.60SD local.slope.60.300Mean local.slope.60.300SD local.slope.300.600Mean local.slope.300.600SD t25Mean t25SD t50Mean t50SD t75Mean t75SD t90Mean t90SD t95Mean t95SD t99Mean t99SD delta.t.max.25Mean delta.t.max.25SD delta.t.25.50Mean delta.t.25.50SD delta.t.50.75Mean delta.t.50.75SD delta.t.75.90Mean delta.t.75.90SD delta.t.90.95Mean delta.t.90.95SD delta.t.95.99Mean delta.t.95.99SD oldAUCfalseMean oldAUCfalseSD AUCMean AUCSD AUCWendcurveMean AUCWendcurveSD AUCWendrectMean AUCWendrectSD AUCWt0curveMean AUCWt0curveSD AUCWt0rectMean AUCWt0rectSD AUCWendsumMean AUCWendsumSD AUCWt0sumMean AUCWt0sumSD AUC.100Mean AUC.100SD AUC.200Mean AUC.200SD AUC.300Mean AUC.300SD AUC.400Mean AUC.400SD AUC.500Mean AUC.500SD AUC.600Mean AUC.600SD AUC.700Mean AUC.700SD AUC.800Mean AUC.800SD AUC.900Mean AUC.900SD AUC.1000Mean AUC.1000SD corgMean corgSD humusMean humusSD pHeauMean pHeauSD clayMean claySD siltfMean siltfSD siltcMean siltcSD siltMean siltSD sandfMean sandfSD sandcMean sandcSD sandMean sandSD ibMean ibSD nlistingMean nlistingSD oc.clayMean oc.claySD tareMean tareSD tarepfraisMean tarepfraisSD tarepsecMean tarepsecSD tarep105Mean tarep105SD pfraisMean pfraisSD psecMean psecSD p105Mean p105SD watercontentMean watercontentSD pctmsMean pctmsSD watercontentcorrMean watercontentcorrSD dappMean dappSD rhobMean rhobSD WrootMean WrootSD dryingjMean dryingjSD X essai traitement MWDT1 MWDT2 MWDT3 MAT1 MAT2 MAT3 MWDmean MWD2.MWD1 MWD3.MWD1 MWD3.MWD2 DiffMWD1MDW2
1 Ernage A-16 PK16 1 Ernage A P1 K0 P-K mineral fertiliser 413 23.8 1080 95.4 10.3 13.2 0.926 0.0472 0.342 0.279 0.074 0.0472 0.584 0.268 0.45 0.0158 0.0194 0.0158 -0.00666 0.00588 -0.00536 0.00402 -0.00197 0.00102 -0.00107 0.000474 -0.00666 0.00588 -0.0041 0.00281 -0.00112 0.000522 -0.000186 7.05e-05 49 35.3 103 68.5 194 134 378 199 539 216 763 177 49 35.3 54.3 35.3 90.1 67.2 184 83.8 161 67.4 224 154 386 250 386 250 85.2 7.57 300 246 1.85 1.61 384 249 386 250 386 250 75 20.5 125 49.3 168 80 207 110 247 139 284 167 321 197 353 224 386 250 345 257 11.9 0 2.38 0 6.86 0 19.9 0 27.5 0 47.4 0 74.9 0 5.06 0 0.214 0 5.27 0 1.76 0 28 0 0.06 0 8.96 0 169 0 140 0 133 0 160 0 131 0 124 0 18.3 0 0.982 0 18.7 0 1.31 0 1.28 0     30 0 2 PK P1 K0 0.264 1.22 0.768 0.443 0.802 0.816 0.75 4.62 2.91 0.631 0.504
2 Ernage A-2 PK02 1 Ernage A P1 K1 P-K mineral fertiliser 424 19.9 2000 2210 2.49 0.541 0.942 0.0386 0.214 0.227 0.0575 0.0386 0.73 0.198 0.435 0.0252 0.0236 0.0154 -0.0136 0.00427 -0.00872 0.00293 -0.0025 0.000707 -0.00143 0.000208 -0.0136 0.00427 -0.00467 0.000929 -0.000987 0.00028 -0.000121 7.9e-05 19.4 2.18 33.9 4.21 81.7 19.6 234 128 339 218 625 298 19.4 2.18 14.4 4.46 47.8 16.7 152 117 106 102 285 277 209 115 254 188 63.3 28.1 191 204 0.282 0.226 254 188 254 188 254 188 58.5 13.5 89.6 33 116 53.8 141 75.6 164 95.4 137 74 149 84.9 160 94.8 170 104 183 117 9.75 0 1.95 0 6.91 0 20 0 28.3 0 46 0 74.2 0 5.52 0 0.286 0 5.8 0 1.95 0 27 0 0.0488 0 8.94 0 171 0 142 0 151 0 162 0 133 0 142 0 17.7 0 0.981 0 18 0 1.33 0 1.31 0     30 0 1 PK P1 K1 0.225 0.997 0.537 0.358 0.762 0.751 0.586 4.44 2.39 0.538 0.312
3 Ernage A-20 PK20 1 Ernage A P1 K2 P-K mineral fertiliser 423 18.9 1200 239 5.95 6.94 0.93 0.0339 0.248 0.108 0.07 0.0339 0.682 0.125 0.448 0.0148 0.0199 0.0115 -0.0066 0.00297 -0.00566 0.00168 -0.00208 0.000327 -0.00113 0.000207 -0.0066 0.00297 -0.00478 0.00165 -0.00117 0.000551 -0.000162 0.000169 37.9 17 77.8 25.2 258 225 526 238 724 263 1070 161 37.9 17 39.9 17.4 180 218 269 226 198 188 341 334 352 81.5 352 81.5 135 54.6 217 94 1.42 2.04 351 80.5 352 81.5 352 81.5 75.1 9.22 121 16.1 161 24.3 198 33.6 233 44.4 268 55.4 298 65.3 326 73.8 352 81.5 372 103 11.4 0 2.28 0 7 0 19.9 0 26.8 0 47.4 0 74.2 0 5.57 0 0.398 0 5.96 0 1.77 0 29 0 0.0574 0 8.96 0 171 0 141 0 140 0 162 0 132 0 131 0 18.1 0 0.983 0 18.5 0 1.32 0 1.3 0     30 0 3 PK P1 K2 0.279 1.17 0.762 0.407 0.803 0.828 0.737 4.2 2.73 0.651 0.483
4 Ernage A-35 PK35 2 Ernage A P2 K1 P-K mineral fertiliser 425 13.7 3170 2200 17.2 11.6 0.896 0.027 0.278 0.122 0.104 0.027 0.618 0.0986 0.436 0.0152 0.0122 0.0147 -0.00266 0.000699 -0.00358 0.0015 -0.00224 0.000404 -0.00117 0.000189 -0.00266 0.000699 -0.00454 0.00244 -0.00194 0.000404 -0.000108 0.000101 72.4 24.9 115 29.6 210 72.6 264 75.8 352 58.8 1620 1630 72.4 24.9 42.5 12.8 94.8 51.9 54.3 25.3 88.2 64.8 1260 1630 351 101 351 101 108 24.7 243 106 3.89 2.33 347 101 351 101 351 101 87.1 7.52 139 23 178 37.8 210 48.8 240 59.1 269 70 297 79.8 324 90.2 351 101 389 134 10.8 0 2.17 0 6.61 0 16.2 0 27.9 0 49.8 0 77.7 0 5.61 0 0.513 0 6.12 0 2.09 0 30 0 0.0668 0 9 0 178 0 146 0 146 0 169 0 137 0 137 0 19 0 0.986 0 19.2 0 1.37 0 1.35 0     30 0 4 PK P2 K1 0.295 0.868 0.53 0.492 0.734 0.758 0.565 2.94 1.79 0.611 0.235
5 Ernage A-39 PK39 2 Ernage A P2 K0 P-K mineral fertiliser 395 11.8 1090 266 8.49 5.5 0.898 0.0432 0.276 0.222 0.102 0.0432 0.622 0.226 0.43 0.01 0.0168 0.0111 -0.00936 0.00627 -0.0075 0.00328 -0.00226 0.000666 -0.00116 0.000347 -0.00936 0.00627 -0.0056 0.00214 -0.000922 0.000204 -9.8e-05 4.67e-05 35.5 15.6 52.1 23.6 101 26.5 253 102 481 223 777 396 35.5 15.6 16.6 8.8 48.6 12.7 152 107 228 193 295 209 308 183 315 179 73.6 24.1 242 193 2.19 1.65 313 178 315 179 315 179 67.7 17.2 107 36.7 141 56.4 173 76.9 203 97.6 233 118 261 138 304 179 332 202 447 214 11.5 0 2.3 0 6.95 0 15 0 28 0 50.3 0 78.3 0 6.14 0 0.583 0 6.73 0 2.1 0 31 0 0.077 0 8.96 0 164 0 136 0 137 0 155 0 127 0 128 0 18.2 0 0.986 0 18.5 0 1.27 0 1.25 0     30 0 5 PK P2 K0 0.319 1.24 0.628 0.508 0.763 0.833 0.728 3.88 1.97 0.508 0.309
6 Ernage A-40 PK40 2 Ernage A P2 K2 P-K mineral fertiliser 392 6.98 1160 444 31.1 12.8 0.875 0.0332 0.485 0.201 0.125 0.0332 0.391 0.229 0.43 0.0141 0.0042 0.000707 -0.0011 0.000882 -0.00156 0.00117 -0.00123 0.000332 -0.000927 0.000239 -0.0011 0.000882 -0.0022 0.00173 -0.00127 0.000289 -0.00048 0.000372 133 50.9 205 101 339 91.6 513 154 730 213 1020 446 133 50.9 72.8 52.9 134 35.8 173 111 218 126 294 239 517 60.9 567 113 138 78.6 430 182 10.9 6.69 556 107 567 113 567 113 95.2 3.19 175 13.8 246 23.4 307 32.5 363 44.8 417 59.9 439 53.9 480 63.9 519 74.3 601 62 11.7 0 2.35 0 6.8 0 16.6 0 30.3 0 47.3 0 77.5 0 5.16 0 0.651 0 5.81 0 2.02 0 32 0 0.0705 0 8.57 0 184 0 150 0 149 0 175 0 141 0 140 0 19.5 0 0.986 0 19.8 0 1.41 0 1.39 0     30 0 6 PK P2 K2 0.312 1.33 0.779 0.482 0.801 0.79 0.807 4.26 2.5 0.587 0.467
7 Ernage A-41 PK41 2 Ernage A P1 K2 P-K mineral fertiliser 414 20.7 1290 418 10.1 1.3 0.864 0.0472 0.302 0.1 0.136 0.0472 0.562 0.076 0.458 0.0192 0.0139 0.00565 -0.00678 0.00268 -0.00646 0.00184 -0.00208 0.000327 -0.00108 0.000187 -0.00678 0.00268 -0.0062 0.00291 -0.000968 0.000576 -0.000106 6.63e-05 40 14.8 61.7 17.2 141 55.9 418 239 651 228 1060 384 40 14.8 21.7 13.5 79.5 39.2 277 199 233 211 408 175 354 84.1 363 75.6 98.2 29.6 265 86 3.41 2.05 360 77.3 363 75.6 363 75.6 74.6 6.92 119 13.8 159 20.6 196 28.7 233 38.5 269 48.1 302 57.4 342 72.2 374 82.4 414 99.1 12.5 0 2.5 0 6.85 0 16.9 0 31.9 0 46.2 0 78.1 0 4.38 0 0.661 0 5.04 0 1.97 0 33 0 0.074 0 8.03 0 176 0 143 0 143 0 168 0 135 0 135 0 19.5 0 0.986 0 19.8 0 1.35 0 1.33 0     30 0 7 PK P1 K2 0.297 1.64 0.642 0.526 0.834 0.844 0.859 5.51 2.16 0.392 0.345
8 Ernage A-47 PK47 2 Ernage A P1 K0 P-K mineral fertiliser 427 1.41 1140 157 31.1 8.34 0.82 0.0707 0.63 0.141 0.18 0.0707 0.19 0.0707 0.44 0 0.0057 0.000707 -0.00175 0.000354 -0.00155 0.000636 -0.000985 0.000163 -0.000595 0.000205 -0.00175 0.000354 -0.00137 0.00103 -0.00084 0.000368 -0.000211 0.000239 80.6 76.9 127 143 277 193 444 304 510 256 655 182 80.6 76.9 46.6 66 150 49.9 168 112 65 48.1 145 74.6 634 46.1 634 46.1 84.5 74.5 550 121 18 16.1 616 62.3 634 46.1 634 46.1 94.4 1.74 175 6.32 249 5 318 0.595 386 7.98 456 11.3 517 20.8 576 33.4 634 46.1 712 63.2 11 0 2.21 0 6.69 0 12.8 0 28.7 0 52.4 0 81 0 5.76 0 0.413 0 6.18 0 2.36 0 34 0 0.0862 0 8.06 0 172 0 140 0 142 0 164 0 132 0 134 0 19.6 0 0.973 0 20.1 0 1.32 0 1.28 0     30 0 8 PK P1 K0 0.272 1.09 0.574 0.472 0.761 0.778 0.647 4.02 2.11 0.525 0.302
9 Ernage A-51 PK51 2 Ernage A P1 K1 P-K mineral fertiliser 434 14.5 967 165 15.5 5.03 0.86 0.0216 0.393 0.247 0.14 0.0216 0.468 0.265 0.442 0.025 0.00948 0.00256 -0.006 0.00428 -0.00595 0.00297 -0.0019 0.000821 -0.000995 0.000438 -0.006 0.00428 -0.00575 0.00315 -0.000882 0.000646 -8.68e-05 4.53e-05 46.6 13.4 69.5 24.2 116 42.6 266 141 374 123 733 239 46.6 13.4 23 11.1 46.8 18.9 150 127 108 76.2 359 120 405 202 407 197 64.2 27 343 216 4.72 3.14 403 195 407 197 407 197 78.7 11.5 127 31.5 172 54 213 77.6 254 102 296 128 333 151 444 114 491 127 554 144 9.69 0 1.94 0 6.69 0 13 0 31.3 0 50.3 0 81.6 0 5.07 0 0.313 0 5.38 0 2.61 0 35 0 0.0743 0 8.09 0 166 0 137 0 139 0 158 0 128 0 131 0 18.5 0 0.987 0 18.7 0 1.28 0 1.27 0     30 0 9 PK P1 K1 0.271 1.15 0.519 0.509 0.781 0.772 0.647 4.24 1.91 0.451 0.247
10 Les Fonds-13NT LF13 5 Les Fonds RT Tillage 579 7.54 3150 2930 46.8 35.3 0.8 0.0548 0.76 0.0898 0.2 0.0548 0.0395 0.107 0.47 0.00816 0.00628 0.00407 -0.00128 0.000483 -0.00121 0.000416 -0.000535 0.000222 -0.000332 0.000135 -0.00128 0.000483 -0.00112 0.000508 -0.000368 0.000221 -0.00013 8.97e-05 97.7 55.7 206 112 589 159 846 154 1160 450 1650 949 97.7 55.7 108 59.8 383 206 257 131 317 317 392 545 740 51.1 740 51.1 74 30.3 666 76.4 61.8 50.9 678 28.3 740 51.1 740 51.1 94.9 1.74 182 5.5 267 10.3 351 16.5 432 23.3 517 30.2 594 35.4 667 43.4 740 51.1 814 85.5 13 0 2.6 0 7.17 0 13.1 0 30.4 0 49.2 0 79.6 0 5.86 0 1.49 0 7.35 0 2.08 0 3 0 0.0994 0 8.08 0 159 0 131 0 129 0 150 0 123 0 121 0 18.4 0 0.988 0 18.6 0 1.23 0 1.21 0 55.8 10.2 30 0 30 Tvl du sol TCS 0.28 1.21 1 0.473 0.729 0.794 0.832 4.34 3.57 0.824 0.72

2.4.3. Raw data from one QST experiment

Geek zonedata.slake.full.example <- read.csv(
   file = c(
     "./data-raw/data-slake-full-example.csv"
   )
)

data.slake.example <- read.csv(
   file = c(
     "./data-raw/data-slake-example.csv"
   )
)

head(data.slake.example, n = 10) %>%
   mutate(across(where(is.double), \(x) signif(x, digits = 3)))
  X.1 sample.id slake i campaign serie sample.x time mass mass.2 time.2 derivmass derivsecmass diffmass.full pdiffmass.full difftime.full difftime difftime.num diffmass pdiffmass sample.y campagne parcelle mode profondeur dryingtemp dryingj eau operator starttime device device2 essai culture series.id subzone X modalite essaib parcellea parcelleb corg a cn hum xxx1 xxx2 xxx3 xxx4 xxx5 xxx6 xxx7 traitement bloc terre plateforme ligne colonne traitment2 traitement3 traitement4 newbloc slake.sample
1 1 ErnageA-20-PK-20-P1K2-2019-D-F 410 11 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:22.739357 54.3 54.3 2019-07-01 10:46:22.739357 -3.88 49.9 -78.7 0.408 8.62 0 0 -7.55 0.878 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       P1K2 1 Ernage A PK     P1 K2     410 - PK-20-P1K2-2019-D-F
2 2 ErnageA-20-PK-20-P1K2-2019-D-F 410 12 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:23.496483 58.7 58.7 2019-07-01 10:46:23.496483 4.42 8.3 -74.3 0.442 9.37 0.757 0.757 -3.13 0.949 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       P1K2 1 Ernage A PK     P1 K2     410 - PK-20-P1K2-2019-D-F
3 3 ErnageA-20-PK-20-P1K2-2019-D-F 410 13 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:24.333913 60 60 2019-07-01 10:46:24.333913 1.31 -3.11 -73 0.451 10.2 1.59 1.59 -1.82 0.971 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       P1K2 1 Ernage A PK     P1 K2     410 - PK-20-P1K2-2019-D-F
4 4 ErnageA-20-PK-20-P1K2-2019-D-F 410 14 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:25.300481 59.9 59.9 2019-07-01 10:46:25.300481 -0.09 -1.4 -73 0.451 11.2 2.56 2.56 -1.91 0.969 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       P1K2 1 Ernage A PK     P1 K2     410 - PK-20-P1K2-2019-D-F
5 5 ErnageA-20-PK-20-P1K2-2019-D-F 410 15 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:26.260727 60.3 60.3 2019-07-01 10:46:26.260727 0.33 0.42 -72.7 0.453 12.1 3.52 3.52 -1.58 0.974 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       P1K2 1 Ernage A PK     P1 K2     410 - PK-20-P1K2-2019-D-F
6 6 ErnageA-20-PK-20-P1K2-2019-D-F 410 16 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:27.217738 60.5 60.5 2019-07-01 10:46:27.217738 0.25 -0.08 -72.5 0.455 13.1 4.48 4.48 -1.33 0.978 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       P1K2 1 Ernage A PK     P1 K2     410 - PK-20-P1K2-2019-D-F
7 7 ErnageA-20-PK-20-P1K2-2019-D-F 410 17 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:28.187114 60.8 60.8 2019-07-01 10:46:28.187114 0.3 0.05 -72.2 0.457 14.1 5.45 5.45 -1.03 0.983 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       P1K2 1 Ernage A PK     P1 K2     410 - PK-20-P1K2-2019-D-F
8 8 ErnageA-20-PK-20-P1K2-2019-D-F 410 18 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:29.140590 60.9 60.9 2019-07-01 10:46:29.14059 0.08 -0.22 -72.1 0.458 15 6.4 6.4 -0.95 0.985 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       P1K2 1 Ernage A PK     P1 K2     410 - PK-20-P1K2-2019-D-F
9 9 ErnageA-20-PK-20-P1K2-2019-D-F 410 19 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:30.096646 61.1 61.1 2019-07-01 10:46:30.096646 0.21 0.13 -71.9 0.46 16 7.36 7.36 -0.74 0.988 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       P1K2 1 Ernage A PK     P1 K2     410 - PK-20-P1K2-2019-D-F
10 10 ErnageA-20-PK-20-P1K2-2019-D-F 410 20 Ernage A 20 PK-20-P1K2-2019-D-F 2019-07-01 10:46:30.901920 61.2 61.2 2019-07-01 10:46:30.90192 0.15 -0.06 -71.7 0.461 16.8 8.16 8.16 -0.59 0.99 PK-20-P1K2-2019-D-F Ernage A 20 K 3-8 Air libre 30 Eau déminéralisée Etudiant 2 20200000000000.0     Ernage A Froment Ernage A-20 20 3 P1 K2 PK PK20 E-20-2019                       P1K2 1 Ernage A PK     P1 K2     410 - PK-20-P1K2-2019-D-F

3. Personal functions

3.1. Illustrative plot of early phase of a QST curve (bw, no data)

Geek zoneillustr.qst.curve.early.nb.nd <- function(my.data.slake, my.qst.opendata) {
  ## data for testing the function
  ## my.data.slake <- data.slake.full.example
  ## my.qst.opendata <- qst.opendata
  ## max.time = my.qst.opendata[as.character(id.slakes), 'tmax']
  id.slakes <- my.data.slake$slake[1]
  max.time = my.qst.opendata[as.character(id.slakes), 'tmax']

  graph <- (my.data.slake %>%
            ggplot2::ggplot(aes(x = difftime.num, y = pdiffmass))
    + geom_ribbon(aes(ymax = pdiffmass),ymin=0, xmin=0,
                  fill="gray60",colour=NA)
    + geom_ribbon(data = my.data.slake[8:11,],aes(ymax = pdiffmass),ymin=0, xmin=my.data.slake$difftime.num[8],
                  fill="gray80",colour=NA)
    + geom_ribbon(data = my.data.slake[1:8,],aes(ymax = pdiffmass),ymin=0,
                  fill="gray90",colour=NA)
    + geom_hline(yintercept = 1, col = "gray15", linewidth = 0.3)
    + geom_hline(yintercept = my.qst.opendata[as.character(id.slakes), 'Wt0'], col = "gray15", linewidth = 0.3)
    + geom_vline(xintercept = my.data.slake$difftime.num[1:100], col = "gray5", linewidth = 0.3)
    + geom_vline(xintercept = my.data.slake$difftime.num[c(1,5,10,15,20,25,30,35)],
                 col = "gray5",
                 linewidth = 0.6)
    + geom_ribbon(aes(ymin=pdiffmass), ymax = 10, fill="white", colour=NA, alpha=0.8)
    + geom_line(colour = 'gray20', linewidth = 3)
    + geom_point(colour = 'gray80', size = 0.5)
    ## t0
    + geom_vline(xintercept = 0, col = "gray15", linewidth = 1)
    + geom_point(x = 0 , y = my.qst.opendata[as.character(id.slakes), 'Wt0'], col = "gray89", size = 3)
    + annotate(geom = "label", x = 0,
               y = my.qst.opendata[as.character(id.slakes), 'Wt0'],
               label = "W[t0]",
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), label.size = 1,
               hjust = 1.05,
               vjust = 1.2
               )
    ## Archimedes buoyancy
    + geom_vline(xintercept = my.data.slake$difftime.num[1],
                 col = "gray15", linewidth = 1)
    + geom_vline(xintercept = my.data.slake$difftime.num[8],
                 col = "gray15", linewidth = 1)
    + geom_segment(x = my.data.slake$difftime.num[3],
                   xend = my.data.slake$difftime.num[3],
                   y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']),
                   yend = my.qst.opendata[as.character(id.slakes), 'Wt0'],
                   col = "gray15",
                   linewidth = 1,
                   arrow = arrow(ends="both", length = unit(0.008, "npc"))
                   )
    + geom_point(x = my.data.slake$difftime.num[3],
                 y = (my.qst.opendata[as.character(id.slakes), 'Wt0']+(1/my.qst.opendata[as.character(id.slakes), 'Warchi']))/2,
                 col = "gray15",
                 size = 2
                 )
    + annotate(geom = "label",
               x = my.data.slake$difftime.num[3],
               y = (my.qst.opendata[as.character(id.slakes), 'Wt0']+(1/my.qst.opendata[as.character(id.slakes), 'Warchi']))/2,
               label = "Archimedes \n buoyancy",
               ## parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), label.size = 1,
               hjust = -0.05
               )
    ## tmax
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 'tmax'], col = "gray15", linewidth = 1)
    + geom_point(x = my.qst.opendata[as.character(id.slakes), 'tmax'], y = 1, col = "gray89", size = 3)
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 'tmax'],
               y = 1.1,
               ## label = paste("t[max]==", round(my.qst.opendata[as.character(id.slakes), 'tmax'],0),"~sec"),
               label = "t[max]",
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), label.size = 1,
               hjust = -0.05
               )
    ## data acquisition (i)
    + annotate(geom = "label", x = my.data.slake$difftime.num[1],
               y = 0.2,
               ## y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi'] + 1.3),
               label = "Acquired data (i):",
               ## size=5,
               label.size=0,
               hjust = 0.1
               )
    + annotate(geom = "text", x = my.data.slake$difftime.num[c(1,5,10,15,20,25,30,35)],
               y = 0,
               ## y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi'] + 1.2),
               label = c(1,5,10,15,20,25,30,35),
               ## size=5,
               angle = 90, hjust = 0, vjust = 0
               )
    ## Under water info
    + annotate(geom = "label",
               x = my.data.slake$difftime.num[1],
               y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']+0.6),
               label = "1: Prepartion phase, \nsoil sample is in the air",
               ## size=3,
               hjust = 0, vjust = 0
               )
    + geom_segment(x = my.data.slake$difftime.num[1],
                   xend = my.data.slake$difftime.num[8],
                   y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi'] + 0.5),
                   yend = (1/my.qst.opendata[as.character(id.slakes), 'Warchi'] + 0.5),
                   col = "gray15",
                   linewidth = 0.3,
                   arrow = arrow(ends="both", length = unit(0.008, "npc"))
                   )
    + annotate(geom = "label",
               x = my.data.slake$difftime.num[8],
               y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']+0.2),
               label = "2: Immersion phase",
               ## size=3,
               hjust = 0, vjust = 0
               )
    + geom_segment(x = my.data.slake$difftime.num[8],
                   xend = my.data.slake$difftime.num[11],
                   y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']+0.1),
                   yend = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']+0.1),
                   col = "gray15",
                   linewidth = 0.3,
                   arrow = arrow(ends="both", length = unit(0.008, "npc"))
                   )
    + annotate(geom = "label",
               x = my.data.slake$difftime.num[11],
               y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']-0.3),
               label = "3: QuantiSlakeTest, \nsoil sample is under the water",
               ## size=3,
               hjust = 0, vjust = 0
               )
    + geom_segment(x = my.data.slake$difftime.num[11],
                   xend = my.data.slake$difftime.num[50],
                   y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']-0.4),
                   yend = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']-0.4),
                   col = "gray15",
                   linewidth = 0.3,
                   arrow = arrow(ends="both", length = unit(0.008, "npc"))
                   )
    + annotate(geom = "label",
               x = my.data.slake$difftime.num[12],
               y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']-0.7),
               label = "Early increase in soil mass",
               ## size=3,
               hjust = 0
               )
    + geom_segment(x = my.data.slake$difftime.num[11],
                   xend = my.qst.opendata[as.character(id.slakes), 'tmax'],
                   y = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']-0.8),
                   yend = (1/my.qst.opendata[as.character(id.slakes), 'Warchi']-0.8),
                   col = "gray15",
                   linewidth = 0.3,
                   arrow = arrow(ends="both", length = unit(0.008, "npc"))
                   )
    ## Slope ?
    + geom_segment(x = 0,
                   xend = my.qst.opendata[as.character(id.slakes), 'tmax'],
                   y = my.qst.opendata[as.character(id.slakes), 'Wt0'],
                   yend = 1,
                   col = "gray89",
                   linewidth = 1
                   )
    + geom_point(x = my.qst.opendata[as.character(id.slakes), 'tmax']/2,
                 y = (my.qst.opendata[as.character(id.slakes), 'Wt0']+1)/2,
                 col = "gray89",
                 size = 2
                 )
    + annotate(geom = "label",
               x = my.qst.opendata[as.character(id.slakes), 'tmax']/2,
               y = (my.qst.opendata[as.character(id.slakes), 'Wt0']+1)/2,
               label = "Slope[0-max]",
               parse = TRUE,
               col = "gray15", fill = "gray88",  label.r = unit(0, "lines"),
               label.size = 1, hjust = 0.5, vjust = 1.2
               )
    ## Legend
    + annotate("rect",
               xmin=240, xmax =1000,
               ymin=0.80,ymax =1.15,
               fill="white", colour="gray16"
               )
    + annotate(geom = "text", x = 250, y = 1.12,
               label = "bold('QuantiSlakeTest curve')",
               parse = TRUE,
               size = 5, hjust = 0
               )
    + annotate(geom = "text", x = 250, y = 1.07,
               label = "bold('and main indicators')",
               size=5, parse = TRUE,
               hjust = 0
               )
    + annotate(geom = "label", x = 280, y = 1.01,  
               label = "'(i) Indicators related to the early increase in soil mass'",
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), size = 1,
               label.size = 3, hjust = 0
               )
    + annotate(geom = "label", x = 280, y = 0.95,
               label = "'(ii) Slopes in the decreasing part of the curve'",
               parse = TRUE,
               col = "gray16", size = 1,
               label.size = 3, hjust = 0
               )
    + annotate(geom = "label", x = 280, y = 0.89,
               label = "(iii) Threshold values of mass loss", size = 3,
               fill = "white",
               hjust = 0
               )
    + annotate(geom = "label", x = 280, y = 0.83,
               label = "bold('(iv) Global indicators')",
               parse = TRUE,
               col = "gray90",
               fill = "gray10",
               label.size = 3, hjust = 0
               )
    + labs(y = "Relative soil mass [-]",
           x = "Time [sec]"# , 
           )
    + coord_cartesian(xlim=c(NA,max.time+2),ylim=c(0,3))
    + scale_x_continuous(breaks=c(0,
                                  round(my.data.slake$difftime.num[1],1),
                                  round(my.data.slake$difftime.num[8],1),
                                  round(my.data.slake$difftime.num[20],1),
                                  ## round(my.qst.opendata[as.character(id.slakes), 'tmax'],0),
                                  round(my.qst.opendata[as.character(id.slakes), 'tmax'],1)
                                  ))
    + scale_y_continuous(breaks=c(0,
                                  my.qst.opendata[as.character(id.slakes), 'Wt0'],
                                  round(1/my.qst.opendata[as.character(id.slakes), 'Warchi'],2),
                                  1))
    + theme_classic(base_size=16)
  )
  print(graph)
}
## IDEA with DERIV'
## graph <- (my.data.slake %>%
##           dplyr::select(difftime.num,pdiffmass,derivmass,derivsecmass) %>%
##           tidyr::pivot_longer(cols = -1) %>%          
##           ggplot2::ggplot(aes(x = difftime.num, y = value))
##   + geom_line(colour = 'gray20', linewidth = 3)
##   + geom_line(colour = 'gray20', linewidth = 3)
##   + geom_vline(xintercept = 0, col = "gray15", linewidth = 0.3)
##   + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 'tmax'], col = "gray15", linewidth = 0.3)
##   ## Legend
##   + annotate("rect",
##              xmin=240, xmax =1000,
##              ymin=0.80,ymax =1.15,
##              fill="white", colour="gray16"
##              )
##   + annotate(geom = "text", x = 250, y = 1.12,
##              label = "bold('QuantiSlakeTest curve')",
##              parse = TRUE,
##              size = 5, hjust = 0
##              )
##   + annotate(geom = "text", x = 250, y = 1.07,
##              label = "bold('and main indicators')",
##              size=5, parse = TRUE,
##              hjust = 0
##              )
##   + annotate(geom = "label", x = 280, y = 1.01,  
##              label = "'(i) Indicators related to the early increase in soil mass'",
##              parse = TRUE,
##              col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), size = 1,
##              label.size = 3, hjust = 0
##              )
##   + annotate(geom = "label", x = 280, y = 0.95,
##              label = "'(ii) Slopes in the decreasing part of the curve'",
##              parse = TRUE,
##              col = "gray16", size = 1,
##              label.size = 3, hjust = 0
##              )
##   + annotate(geom = "label", x = 280, y = 0.89,
##              label = "(iii) Threshold values of mass loss", size = 3,
##              fill = "white",
##              hjust = 0
##              )
##   + annotate(geom = "label", x = 280, y = 0.83,
##              label = "bold('(iv) Global indicators')",
##              parse = TRUE,
##              col = "gray90",
##              fill = "gray10",
##              label.size = 3, hjust = 0
##              )
##   + labs(y = "Relative soil mass [-]",
##          x = "Time [sec]"# , 
##          )
##   + coord_cartesian(xlim=c(NA,max.time+2))
##   + scale_x_continuous(breaks=c(0,
##                                 ## round(my.qst.opendata[as.character(id.slakes), 'tmax'],0),
##                                 round(my.qst.opendata[as.character(id.slakes), 'tmax'],0)
##                                 ))
##   + scale_y_continuous(breaks=c(0,
##                                 my.qst.opendata[as.character(id.slakes), 'Wt0'],
##                                 my.qst.opendata[as.character(id.slakes), 'Wend'],
##                                 1)
##   + facet_grid(name~., scales = "free_y")
##   + theme_classic(base_size=16)

3.2. Illustrative plot of one QST curve (bw, no data)

Geek zone
illustr.qst.curve.indic.nb.nd <- function(my.data.slake, my.qst.opendata,  max.time = 1000) {

  id.slakes <- my.data.slake$slake[1]

  graph <- (ggplot2::ggplot(my.data.slake, aes(x = difftime.num, y = pdiffmass))
    + geom_ribbon(aes(ymax =pdiffmass),ymin=0, fill="gray60",colour=NA)
    + geom_vline(xintercept = 0, col = "gray15", linewidth = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 'tmax'], col = "gray15", linewidth = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't25'], col = "gray15", linewidth = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't50'], col = "gray15", linewidth = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't75'], col = "gray15", linewidth = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't90'], col = "gray15", linewidth = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't95'], col = "gray15", linewidth = 0.3)
    + geom_vline(xintercept = my.qst.opendata[as.character(id.slakes), 't99'], col = "gray15", linewidth = 0.3)
    + geom_hline(yintercept = 1, col = "gray15", linewidth = 0.3)
    + geom_hline(yintercept = my.qst.opendata[as.character(id.slakes), 'Wend'], col = "black", linewidth = 0.3)
    + geom_hline(yintercept = my.qst.opendata[as.character(id.slakes), 'Wt0'], col = "gray15", linewidth = 0.3)
    + geom_ribbon(aes(ymin=pdiffmass),ymax =2, fill="white",colour=NA,alpha=0.8)
    + annotate("rect", xmin=1000,xmax =1200,ymin=0.5,ymax =2, fill="white",colour=NA,alpha=0.8)
    + geom_line(colour = 'gray20', linewidth = 3)
    + geom_point(x = 0 , y = my.qst.opendata[as.character(id.slakes), 'Wt0'], col = "gray89", size = 3)
    + geom_point(x = my.qst.opendata[as.character(id.slakes), 'tmax'], y = 1, col = "gray89", size = 3)
    + geom_segment(x = my.qst.opendata[as.character(id.slakes), 'tmax'],
                   xend = my.data.slake$difftime.num[
                                          my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'
                                                                                        ] + 60)][1],
                   y = 1,
                   yend = my.data.slake$pdiffmass[
                                          my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'
                                                                                        ] + 60)][1],
                   col = "gray89",
                   linewidth = 0.6
                   )
    + geom_segment(x = my.qst.opendata[as.character(id.slakes), 'tmax'],
                   xend = my.data.slake$difftime.num[
                                          my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'
                                                                                        ] + 300)][1],
                   y = 1,
                   yend = my.data.slake$pdiffmass[
                                          my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'
                                                                                        ] + 300)][1],
                   col = "gray66",
                   linewidth = 0.6
                   )
    + geom_segment(x = my.data.slake$difftime.num[
                                          my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'
                                                                                        ] + 60)][1],
                   xend = my.data.slake$difftime.num[
                                          my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'
                                                                                        ] + 300)][1],
                   y = my.data.slake$pdiffmass[
                                          my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'
                                                                                        ] + 60)][1],
                   yend = my.data.slake$pdiffmass[
                                          my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'
                                                                                        ] + 300)][1],
                   col = "gray66",
                   linewidth = 0.6
                   )
    + geom_segment(x = my.qst.opendata[as.character(id.slakes), 'tmax'],
                   xend = my.data.slake$difftime.num[
                                          my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'
                                                                                        ] + 600)][1],
                   y = 1,
                   yend = my.data.slake$pdiffmass[
                                          my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'
                                                                                        ] + 600)][1],
                   col = "gray66",
                   linewidth = 0.6
                   )
    + geom_segment(x = my.data.slake$difftime.num[
                                       my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'
                                                                                     ] + 300)][1],
                   xend = my.data.slake$difftime.num[
                                          my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'
                                                                                        ] + 600)][1],
                   y = my.data.slake$pdiffmass[
                                       my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'
                                                                                     ] + 300)][1],
                   yend = my.data.slake$pdiffmass[
                                          my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'
                                                                                        ] + 600)][1],
                   col = "gray66",
                   linewidth = 0.6
                   )
    + geom_point(x = my.qst.opendata[as.character(id.slakes), 'tmax']/2,
                 y = (my.qst.opendata[as.character(id.slakes), 'Wt0']+1)/2,
                 col = "gray89",
                 size = 1
                 )
    + geom_segment(x = 0,
                   xend = my.qst.opendata[as.character(id.slakes), 'tmax'],
                   y = my.qst.opendata[as.character(id.slakes), 'Wt0'],
                   yend = 1,
                   col = "gray89",
                   linewidth = 0.3
                   )
    + geom_segment(x = my.qst.opendata[as.character(id.slakes), 'tmax']/2,
                   y = (my.qst.opendata[as.character(id.slakes), 'Wt0']+1)/2,
                   xend = -10,
                   yend = 1.1,
                   col = "gray15",
                   linewidth = 0.3
                   )
    + annotate(geom = "label", -10,
               y = 1.1,
               label = "Slope[0-max]",
               parse = TRUE,
               col = "gray15", fill = "gray88",  label.r = unit(0, "lines"),## , size = 1,
               label.size = 1, hjust = 1
               )
    + geom_segment(x = -40,
                   xend = -40,
                   y = 1,
                   yend = my.qst.opendata[as.character(id.slakes), 'Wt0'],
                   col = "gray15",
                   linewidth = 0.3,
                   arrow = arrow(ends="both", length = unit(0.008, "npc"))
                   )
    + geom_segment(x = my.qst.opendata[as.character(id.slakes), 't50'],
                   xend = my.qst.opendata[as.character(id.slakes), 't75'],
                   y = 0.1,
                   yend = 0.1,
                   col = "gray15",
                   linewidth = 0.3,
                   arrow = arrow(ends="both", length = unit(0.008, "npc"))
                   )
    + geom_segment(x = my.qst.opendata[as.character(id.slakes), 't75'],
                   xend = my.qst.opendata[as.character(id.slakes), 't90'],
                   y = 0.05,
                   yend = 0.05,
                   col = "gray15",
                   linewidth = 0.3,
                   arrow = arrow(ends="both", length = unit(0.008, "npc"))
                   )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't50'],
               y = 0.1,
               label = "dt[50-75]", 
               parse=TRUE,
               label.size = 1,
               fill = "white",
               hjust = 1.1
               )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't90'],
               y = 0.05,
               label = "dt[75-90]",
               parse=TRUE,
               label.size = 1,
               fill = "white",
               hjust = -0.05
               )
    + geom_point(x = -40,
                 y = mean(c(1,my.qst.opendata[as.character(id.slakes), 'Wt0'])),
                 col = "gray15",
                 size = 1
                 )
    + geom_point(x = my.qst.opendata[as.character(id.slakes), 'tmax'],
                 y = 1.1,
                 col = "gray15",
                 size = 1
                 )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 'tmax'],
               y = 1.1,
               label = "t[max]",
               label.size = 1,
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), ## size = 1,
               hjust = -0.05
               )
    + geom_point(x = my.qst.opendata[as.character(id.slakes), 't25'],
                 y = 0.30,
                 size = 1
                 )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't25'],
               y = 0.30,
               label = "t25", 
               label.size = 1,
               fill = "white",
               hjust = -0.05
               )
    + geom_point(x = my.qst.opendata[as.character(id.slakes), 't50'],
                 y = 0.25,
                 size = 1
                 )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't50'],
               y = 0.25,
               label = "t50",
               label.size = 1,
               fill = "white",             
               hjust = -0.05
               )
    + geom_point(x = my.qst.opendata[as.character(id.slakes), 't75'],
                 y = 0.2,
                 size = 1
                 )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't75'],
               y = 0.2,
               label = "t75", 
               label.size = 1,
               fill = "white",
               hjust = -0.05
               )
    + geom_point(x = my.qst.opendata[as.character(id.slakes), 't90'],
                 y = 0.19,
                 size = 1
                 )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't90'],
               y = 0.19,
               label = "t90",
               label.size = 1,
               fill = "white",
               hjust = -0.05
               )
    + geom_point(x = my.qst.opendata[as.character(id.slakes), 't95'],
                 y = 0.19,
                 size = 1
                 )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't95'],
               y = 0.19,
               label = "t95", 
               label.size = 1,
               fill = "white",
               hjust = -0.05
               )
    + geom_point(x = my.qst.opendata[as.character(id.slakes), 't99'],
                 y = 0.19,
                 size = 1
                 )
    + annotate(geom = "label", x = my.qst.opendata[as.character(id.slakes), 't99'],
               y = 0.19,
               label = "t99",
               label.size = 1,
               fill = "white",
               hjust = 1.05
               )
    + geom_point(x = -100,
                 y = 1,
                 col = "gray5",
                 size = 1
                 )
    + annotate(geom = "label", x = -100,
               y = 1,
               label = "W[max]",
               label.size = 1,
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), # size = 1,
               hjust = 1.05,
               vjust = -0.2
               )
    + geom_point(x = -100,
                 y = my.qst.opendata[as.character(id.slakes), 'Wt0'],
                 col = "gray5",
                 size = 1
                 )
    + annotate(geom = "label", x = -100,
               y = my.qst.opendata[as.character(id.slakes), 'Wt0'],
               label = "W[t0]",
               label.size = 1,
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), # size = 1,
               hjust = 1.05,
               vjust = 1.2
               )
    + geom_point(x = my.data.slake$difftime.num[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)][1],
                 y = my.data.slake$pdiffmass[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)][1],
                 col = "gray89", size = 3)
    + geom_point(x = my.data.slake$difftime.num[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 600)][1],
                 y = my.data.slake$pdiffmass[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 600)][1],
                 col = "gray89", size = 3)
    + geom_point(x = my.data.slake$difftime.num[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 60)][1],
                 y = my.data.slake$pdiffmass[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 60)][1],
                 col = "gray89", size = 3)
    + annotate(geom = "label", x = -50,
               y = mean(c(1,my.qst.opendata[as.character(id.slakes), 'Wt0'])),
               label = "W[max]-W[t0]",
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), ## size = 1,
               label.size = 1, hjust = 1.05,
               )
    + geom_point(x = -50,
                 y = my.qst.opendata[as.character(id.slakes), 'Wend'],
                 col = "black",
                 size = 1
                 )
    + annotate(geom = "label", x = -50,
               y = my.qst.opendata[as.character(id.slakes), 'Wend'],
               label = "bold('W')[bold('end')]",
               label.size = 1,
               parse = TRUE,
               col = "gray90",
               ## fill = "lightgrey",
               fill = "gray10",
               hjust = 1.05,
               vjust = -0.05
               )
    + geom_point(
        x = mean(c(as.numeric(
                     my.data.slake$difftime.num[
                                     my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes),
                                                                                   'tmax'] + 60)
                                   ][1]),
                   as.numeric(
                     my.data.slake$difftime.num[
                                     my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes),
                                                                                   'tmax'] + 300)
                                   ][1]))
                 ),
        y = mean(
          c(my.data.slake$pdiffmass[
                               my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 60)
                             ][1],
            my.data.slake$pdiffmass[
                               my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)
                             ][1])
        ),
                 col = "gray66",
                 size = 1
      )
    + annotate(geom = "label",
               x = mean(c(as.numeric(
                 my.data.slake$difftime.num[
                                 my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes),
                                                                               'tmax'] + 60)
                               ][1]),
                 as.numeric(
                   my.data.slake$difftime.num[
                                   my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes),
                                                                                 'tmax'] + 300)
                                 ][1]))
                 ),
               y = mean(
                 c(my.data.slake$pdiffmass[
                                   my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 60)
                                 ][1],
                   my.data.slake$pdiffmass[
                                   my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)
                                 ][1])
               ),
               label = "Slope[60-300]", 
               parse = TRUE,
               col = "gray16", ## size = 1,
               label.size = 0.3, hjust = 1.1, vjust = 1.1
               )
    + geom_point(
        x = mean(c(as.numeric(
          my.data.slake$difftime.num[
                                     my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes),
                                                                                   'tmax'] + 600)
                                   ][1]),
                   as.numeric(
                     my.data.slake$difftime.num[
                                     my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes),
                                                                                   'tmax'] + 300)
                                   ][1]))
                 ),
        y = mean(
          c(my.data.slake$pdiffmass[
                               my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 600)
                             ][1],
            my.data.slake$pdiffmass[
                               my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)
                             ][1])
        ),
                 col = "gray66",
                 size = 1
      )
    + annotate(geom = "label",
               x = mean(c(as.numeric(
                 my.data.slake$difftime.num[
                                 my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes),
                                                                               'tmax'] + 600)
                               ][1]),
                 as.numeric(
                   my.data.slake$difftime.num[
                                   my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes),
                                                                                 'tmax'] + 300)
                                 ][1]))
                 ),
               y = mean(
                 c(my.data.slake$pdiffmass[
                                   my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 600)
                                 ][1],
                   my.data.slake$pdiffmass[
                                   my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)
                                 ][1])
               ),
               label = "Slope[300-600]", 
               parse = TRUE,
               col = "gray16", ## size = 1,
               label.size = 0.3, hjust = 0.5, vjust = 1.1
               )
    + geom_point(
        x = mean(c(my.qst.opendata[as.character(id.slakes), 'tmax'],
                   as.numeric(
                     my.data.slake$difftime.num[
                                     my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes),
                                                                                   'tmax'] + 60)
                                   ][1]))
                 ),
        y = mean(
          c(1, my.data.slake$pdiffmass[
                               my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 60)
                             ][1])
        ),
                 col = "gray89",
                 size = 1
                 )
    + annotate(geom = "label",
               x = mean(c(my.qst.opendata[as.character(id.slakes), 'tmax'], as.numeric(my.data.slake$difftime.num[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 60)][1]))),
               y = mean(c(1, my.data.slake$pdiffmass[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 60)][1])),
               label = "Slope[max-60]", 
               parse = TRUE,
               col = "gray16", ## size = 1,
               label.size = 0.3, hjust = -0.05, vjust = -0.05
               )
    + geom_point(x = quantile(c(my.qst.opendata[as.character(id.slakes), 'tmax'], as.numeric(my.data.slake$difftime.num[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)][1])),0.5),
                 y = quantile(c(1, my.data.slake$pdiffmass[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)][1]),0.5),
                 col = "gray66",
                 size = 1
                 )
    + annotate(geom = "label",
               x = quantile(c(my.qst.opendata[as.character(id.slakes), 'tmax'], as.numeric(my.data.slake$difftime.num[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)][1])),0.5),
               y = quantile(c(1, my.data.slake$pdiffmass[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 300)][1]),0.5),
               label = "Slope[max-300]", 
               parse = TRUE,
               col = "gray16", ## size = 1,
               label.size = 0.3, hjust = -0.05, vjust = -0.05
               )
    + geom_point(x = quantile(c(my.qst.opendata[as.character(id.slakes), 'tmax'], as.numeric(my.data.slake$difftime.num[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 600)][1])),0.75),
                 y = quantile(c(1, my.data.slake$pdiffmass[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 600)][1]),0.25),
                 col = "gray66",
                 size = 1
                 )
    + annotate(geom = "label",
               x = quantile(c(my.qst.opendata[as.character(id.slakes), 'tmax'],
                             as.numeric(my.data.slake$difftime.num[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 600)][1])), 0.75),
               y = quantile(c(1, my.data.slake$pdiffmass[my.data.slake$difftime.num > (my.qst.opendata[as.character(id.slakes), 'tmax'] + 600)][1]),0.25),
               label = "Slope[max-600]", 
               parse = TRUE,
               col = "gray16", ## size = 1,
               label.size = 0.3, hjust = -0.05, vjust = -0.05
               )
    + annotate(geom = "label", x = 500, y = 0.25,
               label = "bold('Area under curve (AUC)')",
               parse = TRUE,
               col = "gray90",
               fill = "gray10",
               label.size = 1, hjust = 0
               )
    ## Legend
    + annotate("rect",
               xmin=260, xmax =1000,
               ymin=0.80,ymax =1.15,
               fill="white", colour="gray16"
               )
    + annotate(geom = "text", x = 270, y = 1.12,
               label = "bold('QuantiSlakeTest curve')",
               parse = TRUE,
               ## label.size = 5,
               hjust = 0
               )
    + annotate(geom = "text", x = 270, y = 1.08,
               label = "bold('and main indicators')",
               parse = TRUE,
               ## label.size = 5,
               hjust = 0
               )
    + annotate(geom = "label", x = 285, y = 1.01,  
               label = "'(i) Indicators of the early increase in soil mass'",
               parse = TRUE,
               col = "gray5", fill = "gray88",  label.r = unit(0, "lines"), # size = 1,
               label.size = 1, hjust = 0
               )
    + annotate(geom = "label", x = 285, y = 0.95,
               label = "'(ii) Slopes in the decreasing part of the curve'",
               parse = TRUE,
               col = "gray16", # size = 1,
               label.size = 0.3, hjust = 0
               )
    + annotate(geom = "label", x = 285, y = 0.89,
               label = "(iii) Threshold values of mass loss", label.size = 1,
               fill = "white",
               hjust = 0
               )
    + annotate(geom = "label", x = 285, y = 0.83,
               label = "bold('(iv) Global indicators')",
               parse = TRUE,
               col = "gray90",
               fill = "gray10",
               label.size = 1, hjust = 0
               )
    + labs(y = "Relative soil mass [-]",
           x = "Time [sec]"# , 
           )
    + coord_cartesian(xlim=c(-200,max.time), ylim = c(0, 1.1))
    + scale_x_continuous(breaks=c(0,
                                  ## round(my.qst.opendata[as.character(id.slakes), 'tmax'],0),
                                  round(my.qst.opendata[as.character(id.slakes), 't50'],0),
                                  round(my.qst.opendata[as.character(id.slakes), 't90'],0),
                                  round(my.qst.opendata[as.character(id.slakes), 't95'],0),
                                  round(my.qst.opendata[as.character(id.slakes), 't99'],0)
                                  ))
    + scale_y_continuous(breaks=c(0,
                                  my.qst.opendata[as.character(id.slakes), 'Wt0'],
                                  my.qst.opendata[as.character(id.slakes), 'Wend'],
                                  1))
    + theme_classic(base_size=16)
  )

print(graph)
}

3.3. Correlation tables

Geek zonef.df.corr <- function(df, var.l, var.c,
                      my.round = 3
                      ){

  ## keep only numeric data
  df.num <- df[, sapply(df, is.numeric)] 

  ## na.omit(df.num)
  mcor <- cor(df.num)

  df.corr <- as.data.frame(round(mcor,
                                 my.round))[var.l, var.c]
  return(df.corr)
}

f.gt.corr.color <- function(df,
                            my.pal = RColorBrewer::brewer.pal(n=8, name="RdYlBu")){
  tab <- df %>%
    dplyr::mutate(" " = "    ") %>%
    tibble::rownames_to_column("rownames") %>%
    gt(rowname_col = "rownames") %>%
    data_color(
      columns = names(df),
      fn = scales::col_numeric(
                         palette = my.pal,
                         domain = c(-1,1)
                       )) 
  return(tab)
}

3.4. Matrix of boxplots

Geek zoneplot.compare.trt <- function(trial = "P-K mineral fertiliser",
                               my.color="black"){ 
    df <- qst.opendata %>%
      dplyr::filter(essai_uk == trial) %>%
      dplyr::select(modalite_uk_abb, parcelle, bloc,
                    Wmax.Wt0,
                    tmax,
                    Sl.max.30,
                    Sl.max.60,
                    Sl.max.300,
                    ## Sl.max.600,              
                    t75,
                    t95,
                    Wend,
                    AUC) %>%
      dplyr::rename(
               "(a) Wmax-Wt0 [-]" = Wmax.Wt0,
               "(b) tmax [sec]" = tmax,
               "(c) Slope max-30 [sec-1]" = Sl.max.30,
               "(d) Slope max-60 [sec-1]" = Sl.max.60,
               "(e) Slope max-300 [sec-1]" = Sl.max.300,
               ## "(e) Slope max-600 [sec-1]" = Sl.max.600,
               "(f) t75 [sec]" = t75,
               "(g) t95 [sec]" = t95,
               "(h) Wend [-]" = Wend,
               "(i) AUC [sec]" = AUC) %>% 
      tidyr::pivot_longer(!c(modalite_uk_abb,parcelle,bloc),
                          names_to = "indicator",
                          values_to = "value"
                          )

    df$indicator <- factor(df$indicator,
                           levels=c(
                             "(a) Wmax-Wt0 [-]",
                             "(b) tmax [sec]",
                             "(c) Slope max-30 [sec-1]",
                             "(d) Slope max-60 [sec-1]",
                             "(e) Slope max-300 [sec-1]",
                             ## "(e) Slope max-600 [sec-1]" ,
                             "(f) t75 [sec]",
                             "(g) t95 [sec]",
                             "(h) Wend [-]",
                             "(i) AUC [sec]"
                           )
                           )

    levels(df$indicator) <- c(
        expression("(a)"~"W"["max"]~"-W"["t0"]~"[-]"),
        expression("(b)"~"t"["max"]~"[sec]"),
        expression("(c)"~"Slope"["max-30"]~"[sec"^{-1}~"]"),
        expression("(d)"~"Slope"["max-60"]~"[sec"^{-1}~"]"),
        expression("(e)"~"Slope"["max-300"]~"[sec"^{-1}~"]"),
        expression("(f)"~"t75"~"[sec]"),
        expression("(g)"~"t95"~"[sec]"),
        expression("(h)"~"W"["end"]~"[-]"),
        expression("(i)"~"AUC"~"[sec]")
        )

    g <- (ggplot(data = df, aes(x=modalite_uk_abb,
                                y=value))
      + geom_boxplot(colour="black",fill="gray80")
      + facet_wrap(.~indicator,
                   nrow=3, scales="free_y",
                   labeller = label_parsed
                   )
      + labs(title = paste0("Trial - ",trial),
             x = "Treatment",
             y = "Value of the QuantiSlakeTest indicator")
      + theme_bw()
    )
}

3.5. Matrix of boxplots - Version 2

Geek zoneplot.compare.trt.v2 <- function(trial = "P-K mineral fertiliser",
                               my.color="black"){ 
    df <- qst.opendata %>%
      dplyr::filter(essai_uk == trial) %>%
      dplyr::select(modalite_uk_abb, parcelle, bloc,
                    tmax,
                    local.slope.30.60,
                    delta.t.50.75,
                    Wend) %>%
      dplyr::rename(
               "(a) tmax [sec]" = tmax,
               "(b) Slope 30-60 [sec-1]" = local.slope.30.60,
               "(c) dt 50-75 [sec]" = delta.t.50.75,
               "(d) Wend [-]" = Wend) %>% 
      tidyr::pivot_longer(!c(modalite_uk_abb,parcelle,bloc),
                          names_to = "indicator",
                          values_to = "value"
                          )

    df$indicator <- factor(df$indicator,
                           levels=c(
                             "(a) tmax [sec]",
                             "(b) Slope 30-60 [sec-1]",
                             "(c) dt 50-75 [sec]",
                             "(d) Wend [-]"
                           )
                           )

    levels(df$indicator) <- c(
      expression("(a)"~"t"["max"]~"[sec]"),
      expression("(b)"~"Slope"["30-60"]~"[sec"^{-1}~"]"),
      expression("(c)"~"dt"["50-75"]~"[sec]"),
      expression("(d)"~"W"["end"]~"[-]")
    )

    g <- (ggplot(data = df, aes(x=modalite_uk_abb,
                                y=value))
      + geom_boxplot(colour="black",fill="gray80")
      + facet_wrap(.~indicator,
                   nrow=3, scales="free_y",
                   labeller = label_parsed
                   )
      + labs(title = paste0("Trial - ",trial),
             x = "Treatment",
             y = "Value of the QuantiSlakeTest indicator")
      + theme_bw()
    )
}

3.6. Matrix of boxplots - Version 3

Geek zoneplot.compare.trt.v3 <- function(trial = "P-K mineral fertiliser",
                               my.color="black"){ 
    df <- qst.opendata %>%
      dplyr::filter(essai_uk == trial) %>%
      dplyr::select(modalite_uk_abb, parcelle, bloc,
                    tmax,
                    local.slope.60.300,
                    delta.t.50.75,
                    Wend) %>%
      dplyr::rename(
               "(a) tmax [sec]" = tmax,
               "(b) Slope 60-300 [sec-1]" = local.slope.60.300,
               "(c) dt 50-75 [sec]" = delta.t.50.75,
               "(d) Wend [-]" = Wend) %>% 
      tidyr::pivot_longer(!c(modalite_uk_abb,parcelle,bloc),
                          names_to = "indicator",
                          values_to = "value"
                          )

    df$indicator <- factor(df$indicator,
                           levels=c(
                             "(a) tmax [sec]",
                             "(b) Slope 60-300 [sec-1]",
                             "(c) dt 50-75 [sec]",
                             "(d) Wend [-]"
                           )
                           )

    levels(df$indicator) <- c(
      expression("(a)"~"t"["max"]~"[sec]"),
      expression("(b)"~"Slope"["60-300"]~"[sec"^{-1}~"]"),
      expression("(c)"~"dt"["50-75"]~"[sec]"),
      expression("(d)"~"W"["end"]~"[-]")
    )

    g <- (ggplot(data = df, aes(x=modalite_uk_abb,
                                y=value))
      + geom_boxplot(colour="black",fill="gray80")
      + facet_wrap(.~indicator,
                   nrow=3, scales="free_y",
                   labeller = label_parsed
                   )
      + labs(title = paste0("Trial - ",trial),
             x = "Treatment",
             y = "Value of the QuantiSlakeTest indicator")
      + theme_bw()
    )
}

3.7. Statistical analyses, model, tests, ANOVA and comparisons of means

Geek zonef.qst.anova  <- function(my.campagne, indic="Wend"){
  res <- list()

  df <- qst.opendata %>%
    dplyr::filter(campagne == my.campagne) %>%
    dplyr::mutate(trmt = as.factor(as.character(modalite_uk_abb))) %>%
    dplyr::mutate(bloc = as.factor(bloc)) %>%
    dplyr::mutate(parcelle = as.factor(parcelle)) %>%
    dplyr::select(indic,bloc,trmt,parcelle)

  ## Test normality
  test.norm <- shapiro.test(df[[indic]])
  res[["Shapiro"]] <-  test.norm$p.value

  test.var <- bartlett.test(df[[indic]] ~ trmt,
                            data=df)

  ## Test homoscedacity
  res[["Bartlett"]] <-  test.var$p.value

  ## ANOVA
  rmodel <- lmer(
    get(indic) ~ trmt + (1|bloc) + (1|parcelle),
    data=df
  )
  mod <- car::Anova(rmodel, test = "F")
  res[["Anova"]] <- mod

  return(res)
}

f.qst.means.comp <- function(my.campagne,indic="Wend"){
  df.lf <- qst.opendata %>%
    dplyr::filter(campagne == my.campagne) %>%
    dplyr::mutate(trmt = as.factor(as.character(modalite_uk_abb))) %>%
    dplyr::mutate(bloc = as.factor(bloc)) %>%
    dplyr::mutate(parcelle = as.factor(parcelle))

  rmodel <- lmer(
    get(indic) ~ trmt + (1|bloc) + (1|parcelle),
    data=df.lf
  )

  t <- emmeans::contrast(
                  emmeans::emmeans(rmodel, ~ trmt),
                  alpha = 0.05,
                  method = "pairwise"
                )  %>%
    as_tibble() %>%
    mutate(across(where(is.double), signif, digits = 3))
  return(t)
}

s.f.qst.means.comp <- function(my.campagne,my.indic="Wend"){
  t <- f.qst.means.comp(my.campagne,my.indic) %>%
    tidyr::separate(contrast,c("group1","group2"),sep = " - ") %>%
    dplyr::mutate(p.value = signif(p.value,2)) %>%
    dplyr::mutate(campagne = my.campagne, indic = my.indic) %>%
    dplyr::select(campagne,indic,group1,group2,p.value) %>%
    as.data.frame()
}

3.8. Scatter plot 2 variables

Geek zonef.scatter.fred <- function(var_x,var_y,database,databasemeansd) {
  g <- (ggplot(data=database,
               aes(x=.data[[var_x]],
                   y=.data[[var_y]])
               )
    + geom_smooth(data=databasemeansd,
                  aes(x=.data[[paste0(var_x,"_Mean")]],
                      y=.data[[paste0(var_y,"_Mean")]],
                      color = essai_uk),
                  method=lm,
                  formula = y~x,
                  ## color="black",
                  se=F)
    + geom_smooth(data=databasemeansd,
                  aes(x=.data[[paste0(var_x,"_Mean")]],
                      y=.data[[paste0(var_y,"_Mean")]],
                      ),
                  method=lm,
                  formula = y~x,
                  color="black",
                  se=F)
    + stat_summary(
        fun = 'mean',
        fun.min = function(x) 0,
        geom = 'point', size = 6,
        aes(color=essai_uk))
    + stat_summary(
        fun = mean,
        fun.min = function(y) mean(y) - sd(y), 
        fun.max = function(y) mean(y) + sd(y), 
        geom ="pointrange", aes(color=essai_uk),
        show.legend = FALSE,linewidth=1)
    + scale_color_manual(
        values = c(viridis::viridis(4)[1:3]),
        name = "Trials",
        drop = FALSE
      )
    ## + geom_jitter(aes(group=essai_uk, color=essai_uk),
    ##               alpha=0.5,size=1, width=0.0005)
    ## + annotate("text",label = "r=0.925",x=0.05,y=0.3,linewidth=6,hjust = 0)
    ## + labs(
    ##     y=expression("W"["max"] - "W"["t0"] ~ "[-]"),
    ##     x="SOC:Clay ratio [-]")
    + theme_classic(base_size = 16)
  )
  return(g)
}

4. Figures and Tables of the article

4.1. Material & Methods

4.1.1. Figure 1 - Material & Methods - Illustration of the material (see montage photo)

illustration_device.png

4.1.2. Figure 2a - Material & Methods - Illustration of a QST curve and main indicators

Geek zonefig2a <- illustr.qst.curve.indic.nb.nd(
  data.slake.example,
  qst.opendata
)
ggsave("./fig/fig02a.pdf",fig2a, width=2200, height=2000, units="px")

illustration_indicators.png

4.1.3. Figure 2b - Material & Methods - Illustration of the early phase of the QST and Archimedes Buoyancy

Geek zonefig2b <- illustr.qst.curve.early.nb.nd(
  data.slake.full.example,
  qst.opendata
)
ggsave("./fig/fig02b.pdf",
       fig2b, width=2200, height=2000, units="px")

illustration_early_phase.png

4.2. Results

4.2.1. Table 1 - Results - Soil properties

  1. Data cooking
    Geek zonedf.soil.prop <- qst.opendata %>%
  dplyr::select(
           campagne,
           essai_uk,
           modalite_uk,
           parcelle,
           clay,
           ## silt_f,
           ## silt_c,
           silt,
           ## sand_f,
           ## sand_c,
           sand,
           corg,
           oc.clay,
           pH_eau,
           rho_b,
           ) %>%
    dplyr::group_by(essai_uk,parcelle) %>%
  dplyr::slice_head() %>%
  dplyr::ungroup() %>%
  dplyr::mutate(plot = 1:length(parcelle)) %>%
  dplyr::rename(
           Treatment = modalite_uk,
           Plot = plot,
           "Clay (<2 μm)" = clay,
           ## "Silt (fine)" = silt_f,
           ## "Silt (coarse)" = silt_c,
           "Silt (2 μm - 50 μm)" = silt,
           ## "Sand (fine)" = sand_f,
           ## "Sand (coarse)" = sand_c,
           "Sand (50 – 2000 μm)" = sand,    
           "Organic C (g.kg^-1)" = corg,
           "SOC:Clay" = oc.clay,
           "pH [-]" = pH_eau,
           "Bulk density [-]" = rho_b
         )

write.csv(df.soil.prop,file="./data-output/check_up_table_soil_prop.csv")

df.soil.prop <- df.soil.prop %>%
  dplyr::select(-c(campagne,essai_uk,parcelle))
  2. Table as an image with gt
    Geek zone    tab1 <- df.soil.prop %>%
      mutate(empty = "") %>%
      gt() %>%
      tab_spanner(
        label = md("Texture (%)"),
        columns = names(df.soil.prop)[2:4]
      ) %>%
      cols_move(
        columns = "empty",
        after = "Treatment"
      ) %>%
      cols_move(
        columns = "Plot",
        after = "empty"
      ) %>%
      cols_move(
        columns = "Treatment",
        after = "Plot"
      ) %>%
      tab_row_group(
        label = "P-K mineral fertiliser",
        rows = 27:35
      ) %>%
      tab_row_group(
        label = "Tillage",
        rows = 19:26
      ) %>%
      tab_row_group(
        label = "Organic matter",
        rows = 1:18
      ) %>%
      fmt_number(
        columns = c("Clay (<2 μm)","Silt (2 μm - 50 μm)",
                     "Sand (50 – 2000 μm)"),
             decimals = 1
           ) %>%
      fmt_number(
             columns = "SOC:Clay",
       decimals = 3
      ) %>%
          fmt_number(
             columns = c("pH [-]","Bulk density [-]"),
             decimals = 2
           ) %>%
      tab_options(table.font.size = 10)

gtsave(tab1,"./fig/tab01.pdf", zoom = 1)

gtsave(tab1,"./fig/tab01.tex")

gtsave(tab1,"./fig/dataframe_soil_properties.png",
             expand=30)

    dataframe_soil_properties.png

4.2.2. Table 2 - Results - Correlation Matrix QST v.s. Le Bissonnais & Soil Properties

  1. Data cooking
    Geek zonedf.comp <- qst.mean.sd.opendata %>%
  dplyr::select(
           MWD_T1,
           MWD_T2,
           MWD_T3,
           MA_T1,
           MA_T2,
           MA_T3,
           ## Diff_MWD1_MDW2,
           ## MWD2.MWD1,
           ## MWD3.MWD1,
           ## MWD3.MWD2,
           corg_Mean,
           clay_Mean,
           oc.clay_Mean,
           pH_eau_Mean,
           rho_b_Mean,
           ## Wroot_Mean,
           Wmax.Wt0_Mean,
           Wt0.Wend_Mean,
           tmax_Mean,
           Sl.0.max_Mean,
           local.slope.max.30_Mean,
           local.slope.30.60_Mean,
           local.slope.60.300_Mean,
           local.slope.300.600_Mean,
           ## Sl.max.60_Mean,
           ## Sl.max.300_Mean,
           ## Sl.max.600_Mean,
           delta.t.max.25_Mean,
           delta.t.25.50_Mean, 
           delta.t.50.75_Mean,
           delta.t.75.90_Mean, 
           t50_Mean,
           t75_Mean,
           t90_Mean,
           t95_Mean,
           Wend_Mean,
           AUC_Mean
           ## AUC_Wend_curve_Mean,
           ## AUC_Wend_rect_Mean,
         ) %>%
  dplyr::rename(
           "MWD 1" = MWD_T1,
           "MWD 3" = MWD_T2, ## NB INVERSION 2 et 3 !!!
           "MWD 2" = MWD_T3, ## NB INVERSION 2 et 3 !!!
           "MA 1" = MA_T1,
           "MA 3" = MA_T2, ## NB INVERSION 2 et 3 !!!
           "MA 2" = MA_T3, ## NB INVERSION 2 et 3 !!!
           ## "MWD 1 - MWD 2" = Diff_MWD1_MDW2,
           ## "MWD3:MWD1" = MWD2.MWD1,
           ## "MWD2:MWD1" = MWD3.MWD1,
           ## "MWD2:MWD3" = MWD3.MWD2,
           SOC = corg_Mean,
           Clay = clay_Mean,
           "SOC:Clay" = oc.clay_Mean,
           pH = pH_eau_Mean,
           "Bulk density" = rho_b_Mean,
           ## "Root biomass" = Wroot_Mean, ## only 1 LTE
           "Slope 0-max" = Sl.0.max_Mean,
           tmax = tmax_Mean,
           "Wmax-Wt0" = Wmax.Wt0_Mean,
           "Wt0-Wend" = Wt0.Wend_Mean,
           "Slope max-30" = local.slope.max.30_Mean,
           "Slope 30-60" = local.slope.30.60_Mean,
           "Slope 60-300" = local.slope.60.300_Mean,
           "Slope 300-600" = local.slope.300.600_Mean,
           ## "Slope max-60" = Sl.max.60_Mean,
           ## "Slope max-300" = Sl.max.300_Mean,
           ## "Slope max-600" = Sl.max.600_Mean,
           "dt max-25" = delta.t.max.25_Mean,
           "dt 25-50" = delta.t.25.50_Mean, 
           "dt 50-75" = delta.t.50.75_Mean,
           "dt 75-90" = delta.t.75.90_Mean,
           t50 = t50_Mean,
           t75 = t75_Mean,
           t90 = t90_Mean,
           t95 = t95_Mean,
           "Wend" = Wend_Mean,
           AUC = AUC_Mean,
           ## "pAUC curve Wend" = AUC_Wend_curve_Mean,
           ## "pAUC rect. Wend" = AUC_Wend_rect_Mean
         )

df.corr.t2 <- f.df.corr(
  df.comp,
  c(
    "Slope 0-max",
    "tmax",
    "Wmax-Wt0",
    "Wt0-Wend",
    "Slope max-30",
    "Slope 30-60",
    "Slope 60-30",
    "Slope 300-600",
    ## "Slope max-60",
    ## "Slope max-300",
    ## "Slope max-600",
    "t50",
    "t75",
    "t90",
    "t95",
    "dt max-25",
    "dt 25-50",
    "dt 50-75",
    "dt 75-90",
    "Wend",
    "AUC"
    ## "pAUC curve Wend",
    ## "pAUC rect. Wend"
  ),
  c(
    "MWD 1",
    "MWD 2",
    "MWD 3",
    ## "MWD 1 - MWD 2",
    "MA 1",
    "MA 2",
    "MA 3",
    "SOC",
    "Clay",
    "SOC:Clay",
    "pH",
    "Bulk density"
  )
)
  2. Table as an image with gt
    Geek zonetab2 <- f.gt.corr.color(df.corr.t2) %>%
  tab_spanner(
    label = md("Le Bissonnais *et al.* (1996)"),
    columns = names(df.corr.t2)[1:6]
  ) %>%
  tab_spanner(
    label = md("Soil properties"),
    columns = names(df.corr.t2)[7:11]
  ) %>%
  cols_move(
    columns = " ",
    after = "MA 3"
  ) %>%
  tab_row_group(
    label = "(iv) QST global indicators",
    rows = row.names(df.corr.t2)[17:18]
  ) %>%
  tab_row_group(
    label = "(iii) QST threshold values of mass loss",
    rows = row.names(df.corr.t2)[9:16]
  ) %>%
  tab_row_group(
    label = "(ii) QST slopes in the decreasing part of the curve",
    rows = row.names(df.corr.t2)[5:8]
  ) %>%
  tab_row_group(
    label = "(i) QST indicators of the early increase in soil mass",
    rows = row.names(df.corr.t2)[1:4]
  )

gtsave(tab2,"./fig/tab02.pdf", zoom = 1)

gtsave(tab2,"./fig/dataframe_cor_qst_biss_soil.png",
       expand=30)

gtsave(tab2,"./fig/dataframe_cor_qst_biss_soil.html")

    dataframe_cor_qst_biss_soil.png

4.2.3. Figure 3 - Results - Scatterplot correlation SOC:CLay X Wmax.Wt0

Geek zoneqst.opendata$essai_uk <- factor(
  qst.opendata$essai_uk,
  levels = c("Organic matter","Tillage","P-K mineral fertiliser")
)

fig3 <- (ggplot(data=qst.opendata,
             aes(x=oc.clay,
                 y=Wmax.Wt0)
             )
  + geom_smooth(data=qst.mean.sd.opendata,
                aes(x=oc.clay_Mean,
                    y=Wmax.Wt0_Mean),
                method=lm,
                formula = y~x,
                color="black",
                se=F)
  + stat_summary(
      fun = 'mean',
      fun.min = function(x) 0,
      geom = 'point', size = 6,
      aes(color=essai_uk))
  + stat_summary(
      fun = mean,
      fun.min = function(y) mean(y) - sd(y), 
      fun.max = function(y) mean(y) + sd(y), 
      geom ="pointrange", aes(color=essai_uk),
      show.legend = FALSE,linewidth=1)
  + scale_color_manual(
      values = c(viridis::viridis(4)[1:3]),
      name = "Trials",
      drop = FALSE
    )
  + geom_jitter(aes(group=essai_uk, color=essai_uk),
                alpha=0.5,size=1, width=0.0005)
  ## + annotate("text",label = "r=0.925",x=0.05,y=0.3,linewidth=6,hjust = 0)
  + labs(
      y=expression("W"["max"] - "W"["t0"] ~ "[-]"),
      x="SOC:Clay ratio [-]")
  + theme_classic(base_size = 16)
)
ggsave("./fig/fig03.pdf",fig3, width=3000, height=2000, units="px")
print(fig3)

xy_wmax_wt0_oc_clay_cropland.png

4.2.4. Figure 4 - Results - Comparison treatments : Organic Matter

Geek zoneg.bp.om <- plot.compare.trt.v3("Organic matter", my.color = viridis::viridis(4)[2])
ggsave("./fig/fig04.pdf",g.bp.om, width=2300, height=1700, units="px")
print(g.bp.om)

boxplot_practices_OM.png

4.2.5. Figure 5 - Results - Curves and boxplot Tillage and Wend

Illustrative figure done using function from the "slaker" R-package, under devlopment : https://frdvnw.gitlab.io/slaker/dev/. Real data from the QST curves of the 3 long term field trials.

illustr_curves_wend_tillage.png

4.2.6. Figure 6 - Results - Comparison treatments : Tillage

Geek zoneg.bp.t <- plot.compare.trt.v2("Tillage", my.color = viridis::viridis(4)[3])
ggsave("./fig/fig06.pdf",g.bp.t, width=2300, height=1700, units="px")
print(g.bp.t)

boxplot_practices_tillage.png

4.3. Discussion

4.3.1. Figure 7 - Discussion - Root Weights

  1. Data cooking
    Geek zoneqst.opendata.roots <- qst.opendata %>%
  dplyr::filter(essai_uk == "Tillage") %>%
  dplyr::select(modalite_uk, Wroot, parcelle,
                Wmax.Wt0,
                Sl.max.30,
                Sl.max.60,
                Sl.max.300,
                Sl.max.600,              
                Wend
                ) 

df.roots.long <- qst.opendata.roots %>%
  tidyr::pivot_longer(
           !c(
              modalite_uk, Wroot, parcelle
            ),
           names_to = "indicator",
           values_to = "value"
         )

df.roots.long$indicator <- factor(
  df.roots.long$indicator,
  levels=c(
    "Wmax.Wt0",
    "Sl.max.30",
    "Sl.max.60",
    "Sl.max.300",
    "Sl.max.600",
    "Wend"
  )
)

levels(df.roots.long$indicator) <- c(
  expression("(a) W"["max"]~"-W"["t0"]~"[-]"),
  expression("(b) Slope"["max-30"]~"[sec"^{-1}~"]"),
  expression("(c) Slope"["max-60"]~"[sec"^{-1}~"]"),
  expression("(d) Slope"["max-300"]~"[sec"^{-1}~"]"),
  expression("(e) Slope"["max-600"]~"[sec"^{-1}~"]"),
  expression("(f) W"["end"]~"[-]")
)
  2. Graph
    Geek zonefig7 <- (ggplot(data=df.roots.long,
             aes(x=Wroot,
                 y=value,
                 color=modalite_uk))
  + geom_point()
  + geom_smooth(method="lm")
  + ggpmisc::stat_poly_eq(
               aes(label = paste(after_stat(rr.label), sep = "**\", \"**")),
               label.x = "right",
               label.y = "bottom",
               vstep=0.07, size = 3
             )
  + facet_wrap(.~indicator,
               nrow=2, scales="free_y",
               ## labeller = labeller(indicator = as_labeller(indic_names,  label_parsed)),
               labeller = label_parsed
               )
  + scale_color_manual(values = c(viridis::turbo(8)[c(2,7)]), name="Treatment")
  + labs(x = "Root biomass (mg)",
         y = "Value of the QuantiSlakeTest indicators")
  + theme_classic()
)
ggsave("./fig/fig07.pdf",fig7, width=3000, height=1500, units="px")
print(fig7)

    scatter_root_weight_indics.png

5. Statistical analyses - MOVED in a seperated file

./report_stat.html

6. Supplementary figures not shown in the article - MOVED in a seperated file

./report_supp_stig.html

7. Session info

The scripts have been run the Wed Jun 7 09:45:36 2023 in R version 4.3.0 (2023-04-21) on a platform x86_64-pc-linux-gnu (64-bit).