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) # xylem kmax por un factor de 0.3 sp_params$xylem_kmax <- sp_params$xylem_kmax*0.3 } # 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)
data.frame( Pkg_version = as.character(packageVersion('medfate')), Site = params$code, Date = format(Sys.Date()) )
remarks
# 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 # 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 = 'PlantAbsorbedSWRLeaf') # plant transpiration by leaf plot(res_simple, type = 'PlantTranspirationLeaf') # reset layout par(mfrow = c(1,1)) }
# check if complex mode has to be performed if (params$transpMode %in% c('complex', 'both')) { # control object control_obj <- medfate::defaultControl() control_obj$verbose <- FALSE control_obj$transpirationMode <- 'Complex' 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') }
# 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) } }
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))
print(plot_swc_layers_gg(models_dfs))
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))
print(plot_res_gg('Eplanttot', models_dfs, soil_object, measuredData, params$transpMode, leaf_norm = TRUE))
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) )
if (params$transpMode %in% c('simple', 'both')) { print(soil_object) simple_input[['below']] } else { print(soil_object) complex_input[['below']] }
if (params$transpMode == 'simple') { print(simple_input) } else { if (params$transpMode == 'complex') { print(complex_input) } else { print(simple_input) print(complex_input) } }
# saving objects file_name <- file.path('Output', packageVersion('medfate')[[1]], params$code, paste0(format(Sys.time(), "%Y%m%d_%H%M"), '_', params$code, '_', 'global_kmax_factor_report_objects.RData')) save(list = ls(all.names = TRUE), file = file_name) # printing session info devtools::session_info()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.