MedfateValidation: Medfate model validation

#' Internal function to plot results for SWC in the simple model
#'
#' To use inside \code{\link{plot_res}} function
#'
#' @param models List with the data frames with the results of the models as
#'   generated by \code{\link{saveRes}}
#'
#' @param soil Soil object used in the model
#'
#' @param measured_data Measured data dataframe
#'
#' @export

plot_swc_simple <- function(models, soil, measured_data) {
  dates <- models[['simple']][['Dates']]
  SWC_vals_simple <- models[['simple']][['W.1']] * soil$Theta_FC[[1]]
  SWC_vals_meas <- measured_data[['SWC']]
  y_limits <- c(
    min(min(SWC_vals_simple, na.rm = TRUE),
        min(SWC_vals_meas, na.rm = TRUE)) - (min(min(SWC_vals_simple, na.rm = TRUE),
                                                 min(SWC_vals_meas, na.rm = TRUE)))*0.05,
    max(max(SWC_vals_simple, na.rm = TRUE),
        max(SWC_vals_meas, na.rm = TRUE)) + (max(max(SWC_vals_simple, na.rm = TRUE),
                                                 max(SWC_vals_meas, na.rm = TRUE)))*0.05
  )

  # check if there are measured values, if not, return an empty plot
  if (all(is.na(SWC_vals_meas))) {
    plot(dates,
         SWC_vals_simple,
         ylab = 'SWC', xlab = 'Date', type = 'n',
         main = 'Simple vs. Measured',
         ylim = y_limits)

    plot(SWC_vals_simple,
         SWC_vals_meas, pch = 20, xlab = 'Simple', ylab = 'Measured',
         type = 'n',
         ylim = y_limits, xlim = y_limits)
  } else {
    plot(dates,
         SWC_vals_simple,
         ylab = 'SWC', xlab = 'Date', type = 'n',
         main = 'Simple vs. Measured',
         ylim = y_limits)
    lines(dates,
          SWC_vals_simple,
          col = 'blue')
    lines(dates,
          SWC_vals_meas, col = 'red', lwd = 2)

    plot(SWC_vals_simple,
         SWC_vals_meas, pch = 20, xlab = 'Simple', ylab = 'Measured',
         ylim = y_limits, xlim = y_limits)
  }
}

#' Internal function to plot results for SWC in the complex model
#'
#' To use inside \code{\link{plot_res}} function
#'
#' @param models List with the data frames with the results of the models as
#'   generated by \code{\link{saveRes}}
#'
#' @param soil Soil object used in the model
#'
#' @param measured_data Measured data dataframe
#'
#' @export

plot_swc_complex <- function(models, soil, measured_data) {
  dates <- models[['complex']][['Dates']]
  SWC_vals_complex <- models[['complex']][['W.1']] * soil$Theta_FC[[1]]
  SWC_vals_meas <- measured_data[['SWC']]
  y_limits <- c(
    min(min(SWC_vals_complex, na.rm = TRUE), min(SWC_vals_meas, na.rm = TRUE)) - (min(min(SWC_vals_complex, na.rm = TRUE), min(SWC_vals_meas, na.rm = TRUE)))*0.05,
    max(max(SWC_vals_complex, na.rm = TRUE), max(SWC_vals_meas, na.rm = TRUE)) + (max(max(SWC_vals_complex, na.rm = TRUE), max(SWC_vals_meas, na.rm = TRUE)))*0.05
  )

  # check if there are measured values, if not, return an empty plot
  if (all(is.na(SWC_vals_meas))) {
    plot(dates,
         SWC_vals_complex,
         ylab = 'SWC', xlab = 'Date', type = 'n',
         main = 'Complex vs. Measured',
         ylim = y_limits)

    plot(SWC_vals_complex,
         SWC_vals_meas, pch = 20, xlab = 'Complex', ylab = 'Measured',
         type = 'n',
         ylim = y_limits, xlim = y_limits)
  } else {
    plot(dates,
         SWC_vals_complex,
         ylab = 'SWC', xlab = 'Date', type = 'n',
         main = 'Complex vs. Measured',
         ylim = y_limits)
    lines(dates,
          SWC_vals_complex,
          col = 'green')
    lines(dates,
          SWC_vals_meas, col = 'red', lwd = 2)

    plot(SWC_vals_complex,
         SWC_vals_meas, pch = 20, xlab = 'Complex', ylab = 'Measured',
         ylim = y_limits, xlim = y_limits)
  }
}

#' Internal function to plot results for SWC comparing both models
#'
#' To use inside \code{\link{plot_res}} function
#'
#' @param models List with the data frames with the results of the models as
#'   generated by \code{\link{saveRes}}
#'
#' @param soil Soil object used in the model
#'
#' @param measured_data Measured data dataframe
#'
#' @export

plot_swc_both <- function(models, soil, measured_data) {
  dates <- models[['simple']][['Dates']]
  SWC_vals_complex <- models[['complex']][['W.1']] * soil$Theta_FC[[1]]
  SWC_vals_simple <- models[['simple']][['W.1']] * soil$Theta_FC[[1]]
  y_limits <- c(
    min(min(SWC_vals_complex, na.rm = TRUE), min(SWC_vals_simple, na.rm = TRUE)) - (min(min(SWC_vals_complex, na.rm = TRUE), min(SWC_vals_simple, na.rm = TRUE)))*0.05,
    max(max(SWC_vals_complex, na.rm = TRUE), max(SWC_vals_simple, na.rm = TRUE)) + (max(max(SWC_vals_complex, na.rm = TRUE), max(SWC_vals_simple, na.rm = TRUE)))*0.05
  )

  plot(dates,
       SWC_vals_simple,
       ylab = 'SWC', xlab = 'Date', type = 'n',
       main = 'Simple vs. Complex',
       ylim = y_limits)
  lines(dates,
        SWC_vals_simple,
        col = 'blue')
  lines(dates,
        SWC_vals_complex, col = 'green')

  plot(SWC_vals_simple,
       SWC_vals_complex, pch = 20, xlab = 'Simple', ylab = 'Complex',
       ylim = y_limits, xlim = y_limits)
}

#' Internal function to plot results for Eplanttot in the simple model
#'
#' To use inside \code{\link{plot_res}} function
#'
#' @param models List with the data frames with the results of the models as
#'   generated by \code{\link{saveRes}}
#'
#' @param measured_data Measured data dataframe
#'
#' @export

plot_eplanttot_simple <- function(models, measured_data) {
  dates <- models[['simple']][['Dates']]
  Eplanttot_vals_simple <- models[['simple']][['Eplanttot']]
  Eplanttot_vals_meas <- measured_data[['Eplanttot']]
  y_limits <- c(
    min(min(Eplanttot_vals_simple, na.rm = TRUE), min(Eplanttot_vals_meas, na.rm = TRUE)) - (min(min(Eplanttot_vals_simple, na.rm = TRUE), min(Eplanttot_vals_meas, na.rm = TRUE)))*0.05,
    max(max(Eplanttot_vals_simple, na.rm = TRUE), max(Eplanttot_vals_meas, na.rm = TRUE)) + (max(max(Eplanttot_vals_simple, na.rm = TRUE), max(Eplanttot_vals_meas, na.rm = TRUE)))*0.05
  )

  # check if there are measured values, if not, return an empty plot
  if (all(is.na(Eplanttot_vals_meas))) {
    plot(dates,
         Eplanttot_vals_simple,
         ylab = 'Eplanttot', xlab = 'Date', type = 'n',
         main = 'Simple vs. Measured',
         ylim = y_limits)

    plot(Eplanttot_vals_simple,
         Eplanttot_vals_meas, pch = 20, xlab = 'Simple', ylab = 'Measured',
         type = 'n',
         ylim = y_limits, xlim = y_limits)
  } else {
    plot(dates,
         Eplanttot_vals_simple,
         ylab = 'Eplanttot', xlab = 'Date', type = 'n',
         main = 'Simple vs. Measured',
         ylim = y_limits)
    lines(dates,
          Eplanttot_vals_simple,
          col = 'blue')
    lines(dates,
          Eplanttot_vals_meas, col = 'red', lwd = 2)

    plot(Eplanttot_vals_simple,
         Eplanttot_vals_meas, pch = 20, xlab = 'Simple', ylab = 'Measured',
         ylim = y_limits, xlim = y_limits)
  }
}

#' Internal function to plot results for Eplanttot in the complex model
#'
#' To use inside \code{\link{plot_res}} function
#'
#' @param models List with the data frames with the results of the models as
#'   generated by \code{\link{saveRes}}
#'
#' @param measured_data Measured data dataframe
#'
#' @export

plot_eplanttot_complex <- function(models, measured_data) {
  dates <- models[['complex']][['Dates']]
  Eplanttot_vals_complex <- models[['complex']][['Eplanttot']]
  Eplanttot_vals_meas <- measured_data[['Eplanttot']]
  y_limits <- c(
    min(min(Eplanttot_vals_complex, na.rm = TRUE), min(Eplanttot_vals_meas, na.rm = TRUE)) - (min(min(Eplanttot_vals_complex, na.rm = TRUE), min(Eplanttot_vals_meas, na.rm = TRUE)))*0.05,
    max(max(Eplanttot_vals_complex, na.rm = TRUE), max(Eplanttot_vals_meas, na.rm = TRUE)) + (max(max(Eplanttot_vals_complex, na.rm = TRUE), max(Eplanttot_vals_meas, na.rm = TRUE)))*0.05
  )

  # check if there are measured values, if not, return an empty plot
  if (all(is.na(Eplanttot_vals_meas))) {
    plot(dates,
         Eplanttot_vals_complex,
         ylab = 'Eplanttot', xlab = 'Date', type = 'n',
         main = 'Complex vs. Measured',
         ylim = y_limits)

    plot(Eplanttot_vals_complex,
         Eplanttot_vals_meas, pch = 20, xlab = 'Complex', ylab = 'Measured',
         type = 'n',
         ylim = y_limits, xlim = y_limits)
  } else {
    plot(dates,
         Eplanttot_vals_complex,
         ylab = 'Eplanttot', xlab = 'Date', type = 'n',
         main = 'Complex vs. Measured',
         ylim = y_limits)
    lines(dates,
          Eplanttot_vals_complex,
          col = 'green')
    lines(dates,
          Eplanttot_vals_meas, col = 'red', lwd = 2)

    plot(Eplanttot_vals_complex,
         Eplanttot_vals_meas, pch = 20, xlab = 'Complex', ylab = 'Measured',
         ylim = y_limits, xlim = y_limits)
  }
}

#' Internal function to plot results for Eplanttot comparing both models
#'
#' To use inside \code{\link{plot_res}} function
#'
#' @param models List with the data frames with the results of the models as
#'   generated by \code{\link{saveRes}}
#'
#' @param soil Soil object used in the model
#'
#' @param measured_data Measured data dataframe
#'
#' @export

plot_eplanttot_both <- function(models, soil, measured_data) {
  dates <- models[['simple']][['Dates']]
  Eplanttot_vals_complex <- models[['complex']][['Eplanttot']]
  Eplanttot_vals_simple <- models[['simple']][['Eplanttot']]
  y_limits <- c(
    min(min(Eplanttot_vals_complex, na.rm = TRUE), min(Eplanttot_vals_simple, na.rm = TRUE)) - (min(min(Eplanttot_vals_complex, na.rm = TRUE), min(Eplanttot_vals_simple, na.rm = TRUE)))*0.05,
    max(max(Eplanttot_vals_complex, na.rm = TRUE), max(Eplanttot_vals_simple, na.rm = TRUE)) + (max(max(Eplanttot_vals_complex, na.rm = TRUE), max(Eplanttot_vals_simple, na.rm = TRUE)))*0.05
  )

  plot(dates,
       Eplanttot_vals_simple,
       ylab = 'Eplanttot', xlab = 'Date', type = 'n',
       main = 'Simple vs. Complex',
       ylim = y_limits)
  lines(dates,
        Eplanttot_vals_simple,
        col = 'blue')
  lines(dates,
        Eplanttot_vals_complex, col = 'green')

  plot(Eplanttot_vals_simple,
       Eplanttot_vals_complex, pch = 20, xlab = 'Simple', ylab = 'Complex',
       ylim = y_limits, xlim = y_limits)
}

#' Internal function to plot results for Ecohort in the simple model
#'
#' To use inside \code{\link{plot_res}} function
#'
#' @param models List with the data frames with the results of the models as
#'   generated by \code{\link{saveRes}}
#'
#' @param measured_data Measured data dataframe
#'
#' @export

plot_cohorts_simple <- function(models, measured_data) {
  coh_names <- as.character(
    na.omit(stringr::str_extract(names(models[['simple']]), '^E_.+'))
  )
  dates <- models[['simple']][['Dates']]

  for (cohort in coh_names) {
    E_vals_simple <- models[['simple']][[cohort]]
    E_vals_meas <- measured_data[[cohort]]

    y_limits <- c(
      min(min(E_vals_simple, na.rm = TRUE), min(E_vals_meas, na.rm = TRUE)) - (min(min(E_vals_simple, na.rm = TRUE), min(E_vals_meas, na.rm = TRUE)))*0.05,
      max(max(E_vals_simple, na.rm = TRUE), max(E_vals_meas, na.rm = TRUE)) + (max(max(E_vals_simple, na.rm = TRUE), max(E_vals_meas, na.rm = TRUE)))*0.05
    )

    if (all(is.na(E_vals_meas))) {
      plot(dates,
           E_vals_simple,
           ylab = cohort, xlab = 'Date', type = 'n',
           main = 'Simple vs. Measured',
           ylim = y_limits)
      plot(E_vals_simple,
           E_vals_meas, pch = 20, xlab = 'Simple', ylab = 'Measured', type = 'n',
           ylim = y_limits, xlim = y_limits)
    } else {
      plot(dates,
           E_vals_simple,
           ylab = cohort, xlab = 'Date', type = 'n',
           main = 'Simple vs. Measured',
           ylim = y_limits)
      lines(dates,
            E_vals_simple,
            col = 'blue')
      lines(dates,
            E_vals_meas, col = 'red', lwd = 2)

      plot(E_vals_simple,
           E_vals_meas, pch = 20, xlab = 'Simple', ylab = 'Measured',
           ylim = y_limits, xlim = y_limits)
    }
  }
}

#' Internal function to plot results for Ecohort in the complex model
#'
#' To use inside \code{\link{plot_res}} function
#'
#' @param models List with the data frames with the results of the models as
#'   generated by \code{\link{saveRes}}
#'
#' @param measured_data Measured data dataframe
#'
#' @export

plot_cohorts_complex <- function(models, measured_data) {
  coh_names <- as.character(
    na.omit(stringr::str_extract(names(models[['complex']]), '^E_.+'))
  )
  dates <- models[['complex']][['Dates']]

  for (cohort in coh_names) {
    E_vals_complex <- models[['complex']][[cohort]]
    E_vals_meas <- measured_data[[cohort]]

    y_limits <- c(
      min(min(E_vals_complex, na.rm = TRUE), min(E_vals_meas, na.rm = TRUE)) - (min(min(E_vals_complex, na.rm = TRUE), min(E_vals_meas, na.rm = TRUE)))*0.05,
      max(max(E_vals_complex, na.rm = TRUE), max(E_vals_meas, na.rm = TRUE)) + (max(max(E_vals_complex, na.rm = TRUE), max(E_vals_meas, na.rm = TRUE)))*0.05
    )

    if (all(is.na(E_vals_meas))) {
      plot(dates,
           E_vals_complex,
           ylab = cohort, xlab = 'Date', type = 'n',
           main = 'Complex vs. Measured',
           ylim = y_limits)
      plot(E_vals_complex,
           E_vals_meas, pch = 20, xlab = 'Complex', ylab = 'Measured', type = 'n',
           ylim = y_limits, xlim = y_limits)
    } else {
      plot(dates,
           E_vals_complex,
           ylab = cohort, xlab = 'Date', type = 'n',
           main = 'Complex vs. Measured',
           ylim = y_limits)
      lines(dates,
            E_vals_complex,
            col = 'green')
      lines(dates,
            E_vals_meas, col = 'red', lwd = 2)

      plot(E_vals_complex,
           E_vals_meas, pch = 20, xlab = 'Complex', ylab = 'Measured',
           ylim = y_limits, xlim = y_limits)
    }
  }
}

#' Internal function to plot results for Ecohort comparing both models
#'
#' To use inside \code{\link{plot_res}} function
#'
#' @param models List with the data frames with the results of the models as
#'   generated by \code{\link{saveRes}}
#'
#' @param soil Soil object used in the model
#'
#' @param measured_data Measured data dataframe
#'
#' @export

plot_cohorts_both <- function(models, measured_data) {
  coh_names <- as.character(
    na.omit(stringr::str_extract(names(models[['simple']]), '^E_.+'))
  )
  dates <- models[['simple']][['Dates']]

  for (cohort in coh_names) {
    E_vals_complex <- models[['complex']][[cohort]]
    E_vals_simple <- models[['simple']][[cohort]]

    y_limits <- c(
      min(min(E_vals_complex, na.rm = TRUE), min(E_vals_simple, na.rm = TRUE)) - (min(min(E_vals_complex, na.rm = TRUE), min(E_vals_simple, na.rm = TRUE)))*0.05,
      max(max(E_vals_complex, na.rm = TRUE), max(E_vals_simple, na.rm = TRUE)) + (max(max(E_vals_complex, na.rm = TRUE), max(E_vals_simple, na.rm = TRUE)))*0.05
    )

    plot(dates,
         E_vals_complex,
         ylab = cohort, xlab = 'Date', type = 'n',
         main = 'Simple vs. Complex',
         ylim = y_limits)
    lines(dates,
          E_vals_simple,
          col = 'blue')
    lines(dates,
          E_vals_complex, col = 'green')

    plot(E_vals_simple,
         E_vals_complex, pch = 20, xlab = 'Simple', ylab = 'Complex',
         ylim = y_limits, xlim = y_limits)
  }
}

#' Plotting the models results
#'
#' Function to generate models plots in order to inspect the results in a
#' visual way
#'
#' @param variable Character indicating the variable to visualize
#'
#' @param models List with the data frames with the results of the models as
#'   generated by \code{\link{saveRes}}
#'
#' @param soil Soil object used in the model
#'
#' @param measured_data Measured data dataframe
#'
#' @param mode Transpiration mode indicated in the models ('simple', 'complex' or 'both')
#'
#' @export

plot_res <- function(variable, models, soil, measured_data, mode) {

  # par
  par(mfcol = c(2,3))

  # SWC
  if (variable == 'SWC') {

    if (mode == 'both') {
      plot_swc_simple(models, soil, measured_data)
      plot_swc_complex(models, soil, measured_data)
      plot_swc_both(models, soil, measured_data)
    }

    if (mode == 'simple') {
      plot_swc_simple(models, soil, measured_data)
    }

    if (mode == 'complex') {
      plot_swc_complex(models, soil, measured_data)
    }
  }

  # E plant total
  if (variable == 'Eplanttot') {

    if (mode == 'both') {
      plot_eplanttot_simple(models, measured_data)
      plot_eplanttot_complex(models, measured_data)
      plot_eplanttot_both(models, measured_data)
    }

    if (mode == 'simple') {
      plot_eplanttot_simple(models, measured_data)
    }

    if (mode == 'complex') {
      plot_eplanttot_complex(models, measured_data)
    }
  }

  # E by cohort
  if (variable == 'E_by_Cohort') {

    if (mode == 'both') {
      plot_cohorts_simple(models, measured_data)
      plot_cohorts_complex(models, measured_data)
      plot_cohorts_both(models, measured_data)
    }

    if (mode == 'simple') {
      plot_cohorts_simple(models, measured_data)
    }

    if (mode == 'complex') {
      plot_cohorts_complex(models, measured_data)
    }
  }

  # reset layout
  par(mfrow = c(1,1))
}