In MalditoBarbudo/MedfateValidation: Medfate model validation

knitr::opts_chunk$set(echo = FALSE)
knitr::opts_chunk$set(error = TRUE)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(warning = FALSE)
knitr::opts_knit$set(root.dir = params$wd)

# libraries
library(medfate)
library(MedfateValidation)
library(dplyr)

# site/plot code
code <- params$code

# raw data
treeData <- load_treeData(code)
shrubData <- load_shrubData(code)
seedData <- load_seedData(code)
customParams <- load_customParams(code)
measuredData <- load_measuredData(code)
meteoData <- load_meteoData(code)
miscData <- load_miscData(code)
soilData <- load_soilData(code)
terrainData <- load_terrainData(code)
remarks <- load_remarks(code)

# model objects
forest_object <- buildForest(treeData, shrubData, seedData, miscData)
soil_object <- buildSoil(soilData)

# SpParams
# check if new or old will be used
if (params$SPParams == 'old') {
  # load
  data('SpParamsMED', package = 'medfate')
  # modify
  sp_params <- SpParams_mod(SpParamsMED, customParams)
}

if (params$SPParams == 'new') {
  # load
  data('newParams', package = 'MedfateValidation')
  # modify
  sp_params <- SpParams_mod(newParams, customParams)
}

# Warning for custom FC
if (!all(is.na((soilData[['FC']])))) {
  warning('Custom FC values has been provided: ', soilData[['FC']][1], ' ...')
}

# dinamic fig size for soilwater
swc_vars <- c(
  'SWC',
  names(measuredData)[stringr::str_detect(names(measuredData), '^SWC_[0-9]$')]
)

constant_mode <- c('both' = 2, 'simple' = 1, 'complex' = 1)

height_chunk <- 3.5*length(swc_vars)*constant_mode[[params$transpMode]]
knitr::opts_chunk$set(fig.height = height_chunk)

General Info

data.frame(
    Pkg_version  = as.character(packageVersion('medfate')),
    Site = params$code,
    Date = format(Sys.Date())
  )

Remarks

remarks

Simple Model

# check if simple mode has to be performed
if (params$transpMode %in% c('simple', 'both')) {
  # control object
  control_obj <- medfate::defaultControl()
  control_obj$verbose <- FALSE
  control_obj$transpirationMode <- 'Simple'
  control_obj$taper <- params$tapering
  control_obj$averageFracRhizosphereResistance <- params$rhizosphere

  # input object
  simple_input <- medfate::forest2swbInput(forest_object, soil_object,
                                           sp_params, control_obj)

  # input modifications (if any)
  simple_input <- inputMod(simple_input, customParams)

  # model run
  res_simple <- medfate::swb(simple_input, soil_object, meteoData)
} else {
  # if simple mode is not performed, report about it
  cat('Simple model has not been selected to be performed')
}

# show the plot only if simple mode must be performed
if (params$transpMode %in% c('simple', 'both')) {
  # config layout
  par(mfrow = c(3,2))
  # P/PET plot
  plot(res_simple, yearAxis = TRUE)
  # ET plot
  plot(res_simple, type = 'ET', yearAxis = TRUE)
  # SWC plot
  plot(res_simple, 'Theta')
  # Stress plot
  plot(res_simple, type = 'PlantStress', yearAxis = TRUE)
  # Short wave radiation leaf
  plot(res_simple, type = 'PlantPhotosynthesisLeaf')
  # plant transpiration by leaf
  plot(res_simple, type = 'PlantTranspirationLeaf')
  # reset layout
  par(mfrow = c(1,1))
}

Complex Model

# check if complex mode has to be performed
if (params$transpMode %in% c('complex', 'both')) {
  # control object
  control_obj <- medfate::defaultControl()
  control_obj$verbose <- TRUE
  control_obj$transpirationMode <- 'Complex'
  control_obj$taper <- params$tapering
  control_obj$averageFracRhizosphereResistance <- params$rhizosphere

  # input object
  # we also rebuild the soil_object to avoid using W data from the previous
  # runs
  soil_object <- buildSoil(soilData)
  complex_input <- medfate::forest2swbInput(forest_object, soil_object,
                                            sp_params, control_obj)

  # input modifications (if any)
  complex_input <- inputMod(complex_input, customParams)

  # model run
  res_complex <- medfate::swb(complex_input, soil_object, meteoData,
                             terrainData$latitude, terrainData$elevation,
                             terrainData$slope, terrainData$aspect)
} else {
  # if simple mode is not performed, report about it
  cat('Complex model has not been selected to be performed')
}

# show the plot only if simple mode must be performed
if (params$transpMode %in% c('complex', 'both')) {
  # config layout
  par(mfrow = c(4,2))
  # P/PET plot
  plot(res_complex, yearAxis = TRUE)
  # ET plot
  plot(res_complex, type = 'ET', yearAxis = TRUE)
  # SWC layers plot
  plot(res_complex, 'Theta')
  # Stress plot
  plot(res_complex, type = 'PlantStress', yearAxis = TRUE)
  # Canopy Energy Balance plot
  plot(res_complex, type = 'CanopyEnergyBalance')
  # Canopy Temp plot
  plot(res_complex, type = 'CanopyTemperature')
  plot(res_complex, type = 'AirTemperature', lty = 2, add = TRUE)
  # Short wave radiation leaf
  plot(res_complex, type = 'PlantAbsorbedSWRLeaf')
  # plant transpiration by leaf
  plot(res_complex, type = 'PlantTranspirationLeaf')
  # reset layout
  par(mfrow = c(1,1))
}

e_meas <- measuredData %>%
  dplyr::summarize_all(funs(all(is.na(.)))) %>%
  as.logical()

e_measured_coh <- names(measuredData)[!e_meas]

if (params$transpMode %in% c('both')) {
  sp_names <- paste0('E_sp', simple_input$above$SP)
  models_dfs <- saveRes(res_simple, res_complex, e_measured_coh,
                        spParams = sp_params,
                        site_code = params$code)
} else {
  if (params$transpMode %in% c('simple')) {
    sp_names <- paste0('E_sp', simple_input$above$SP)
    models_dfs <- saveRes(res_simple, e_measured_coh,
                          spParams = sp_params,
                          site_code = params$code)
  } else {
    sp_names <- paste0('E_sp', complex_input$above$SP)
    models_dfs <- saveRes(complex_res = res_complex, e_measured_coh,
                          spParams = sp_params, site_code = params$code)
  }
}

Models comparision

SWC in the measured layers

SWC_stats <- statistics_summary('SWC', models_dfs, measuredData,
                                soil_object)
knitr::kable(SWC_stats, row.names = FALSE)

print(plot_res_gg('SWC', models_dfs, soil_object, measuredData, params$transpMode))

SWC in all layers

print(plot_swc_layers_gg(models_dfs))

Total Tranpiration

Eplanttot_stats <- statistics_summary('Eplanttot', models_dfs, measuredData,
                                      soil_object)

Comparision | MAE | r² | Bias ------------|-----|----|------- Simple vs. Measured | r Eplanttot_stats$Eplanttot_MAE_simple | r Eplanttot_stats$Eplanttot_r_sq_simple | r Eplanttot_stats$Eplanttot_bias_simple Complex vs. Measured | r Eplanttot_stats$Eplanttot_MAE_complex | r Eplanttot_stats$Eplanttot_r_sq_complex | r Eplanttot_stats$Eplanttot_bias_complex Simple vs. Complex | r Eplanttot_stats$Eplanttot_MAE_both | r Eplanttot_stats$Eplanttot_r_sq_both | r Eplanttot_stats$Eplanttot_bias_both

# dinamic fig size for total transpiration

if (params$transpMode == 'both') {
  height_chunk <- 9
} else {
  height_chunk <- 3
}

knitr::opts_chunk$set(fig.height = height_chunk)

print(plot_res_gg('Eplanttot', models_dfs, soil_object, measuredData, params$transpMode))

Total transpiration normalized to leaf area

print(plot_res_gg('Eplanttot', models_dfs, soil_object,
                  measuredData, params$transpMode, leaf_norm = TRUE))

Transpiration by cohorts

Ecohorts_stats <- statistics_summary('E_by_Cohort', models_dfs, measuredData,
                                     soil_object)

rbind(
  round(Ecohorts_stats[["Esp_MAE_simple"]], 5),
  round(Ecohorts_stats[["Esp_r_sq_simple"]], 5),
  round(Ecohorts_stats[["Esp_bias_simple"]], 5),
  round(Ecohorts_stats[["Esp_MAE_complex"]], 5),
  round(Ecohorts_stats[["Esp_r_sq_complex"]], 5),
  round(Ecohorts_stats[["Esp_bias_complex"]], 5),
  round(Ecohorts_stats[["Esp_MAE_both"]], 5),
  round(Ecohorts_stats[["Esp_r_sq_both"]], 5),
  round(Ecohorts_stats[["Esp_bias_both"]], 5)
) %>%
  cbind(
    c('MAE Simple vs. Measured',
      'rsq Simple vs. Measured',
      'Bias Simple vs. Measured',
      'MAE Complex vs. Measured',  
      'rsq Complex vs. Measured',  
      'Bias Complex vs. Measured',
      'MAE Simple vs. Complex',
      'rsq Simple vs. Complex',
      'Bias Simple vs. Complex'),
    .
  ) %>%
  as.data.frame() %>%
  knitr::kable(col.names = c('Statistics', Ecohorts_stats[['Cohort_name']]))

# dinamic fig size for total transpiration

if (params$transpMode %in% c('both', 'simple')) {
  ncoh <- ncol(res_simple[["PlantTranspiration"]])
} else {
  ncoh <- ncol(res_complex[["PlantTranspiration"]])
}

if (ncoh < 4) {
  height_chunk = 3
} else {
  height_chunk <- 3*ceiling(ncoh/3)
}

knitr::opts_chunk$set(fig.height = height_chunk)

plots_cohorts <- plot_res_gg('E_by_Cohort', models_dfs, soil_object,
                             measuredData, params$transpMode)

purrr::walk(
  plots_cohorts,
  ~ print(.x)
)

Transpiration by cohorts normalized to leaf area

plots_cohorts_norm <- plot_res_gg('E_by_Cohort', models_dfs, soil_object,
                                  measuredData, params$transpMode,
                                  leaf_norm = TRUE)

purrr::walk(
  plots_cohorts_norm,
  ~ print(.x)
)

Soil Info

print(soil_object)

SoilWaterBalance Input

if (params$transpMode == 'simple') {
  print(simple_input)
} else {
  if (params$transpMode == 'complex') {
    print(complex_input)
  } else {
    print(simple_input)
    print(complex_input)
  }
}

Session Info

# saving objects
file_name <- file.path('Output', packageVersion('medfate')[[1]],
                       params$code,
                       paste0(format(Sys.time(), "%Y%m%d_%H%M"),
                              '_', params$code, '_',
                              'global_report_objects.RData'))
save(list = ls(all.names = TRUE), file = file_name)
# printing session info
devtools::session_info()

MalditoBarbudo/MedfateValidation documentation built on May 7, 2019, 1:22 p.m.