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)
}

General Info

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

Remarks

Site Aspect: r terrainData$aspect

Complex Model

# check if complex mode has to be performed
if (params$transpMode %in% c('complex', 'both')) {

  # get the mean wind value for correct dates
  wind_start <- which(meteoData$DOY == 1)
  wind_end <- length(meteoData$DOY)

  wind_mean <- mean(meteoData$WindSpeed[wind_start:wind_end], na.rm = TRUE)

  # control object
  control_obj <- medfate::defaultControl()
  control_obj$verbose <- FALSE
  control_obj$transpirationMode <- 'Complex'
  control_obj$defaultWindSpeed <- wind_mean
  control_obj$taper <- params$tapering

  # 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')
}

Model SWB plots

# show the plot only if simple mode must be performed
if (params$transpMode %in% c('complex', 'both')) {
  # config layout
  par(mfrow = c(3,2))
  # P/PET plot
  plot(res_complex, yearAxis = TRUE)
  # ET plot
  plot(res_complex, type = 'ET', yearAxis = TRUE)
  # SWC plot
  plot(x = 1:nrow(meteoData),
       y = res_complex$SoilWaterBalance$W.1 * soil_object$Theta_FC[[1]],
       col = 'black', ylab = 'SWC', xlab = 'Days (absolute)')
  # Stress plot
  plot(res_complex, type = 'PlantStress', yearAxis = TRUE)
  # radiation plot
  plot(meteoData$DOY, meteoData$Radiation)
  # wind plot
  barplot(meteoData$WindSpeed)
  # reset layout
  par(mfrow = c(1,1))
}

Model TB plots

if (params$transpMode %in% c('complex', 'both')) {
  # config layout
  par(mfrow = c(2,1))
  # SoilEnergyBalance
  plot(res_complex, 'SoilEnergyBalance')
  # CanopyEnergyBalance
  plot(res_complex, type = 'CanopyEnergyBalance')
  # reset layout
  par(mfrow = c(1,1))
}

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

Temperature Validation

Temp stats

IMPORTANT! Temperature statistics are truncated to start in January to avoid meteoland wind prediction problems in the first months

trunc_st <- which(meteoData$DOY == 1)
trunc_end <- length(meteoData$DOY)

Temp_stats <- statistics_summary('Temperature', models_dfs,
                                 measuredData, meteo_data = meteoData,
                                 trunc = trunc_st:trunc_end)
knitr::kable(Temp_stats, row.names = FALSE)

Temp plots

plots_list <- plot_temp_complex_gg(models_dfs, measuredData, meteoData)

temp_plots <- list(
  complex_ecoh = cowplot::plot_grid(
    plotlist = plots_list[['temp']],
    ncol = 1,
    align = 'h', axis = "tblr"
  ),
  complex_cor = cowplot::plot_grid(
    plotlist = plots_list[['cor']],
    ncol = 1,
    align = 'h', axis = "tblr"
  )
)

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

Soil Info

if (params$transpMode %in% c('simple', 'both')) {
  print(soil_object)
  simple_input[['below']]
} else {
  print(soil_object)
  complex_input[['below']]
}

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, '_',
                              'temperature_report_objects.RData'))
save(list = ls(all.names = TRUE), file = file_name)
# printing session info
devtools::session_info()