R/SimulationPlots.R

Defines functions SimulationPlots

Documented in SimulationPlots

#' Plot figures associated with simulation of landcover spread
#'
#' This plots the results of the landcover simulation data created by `LandCoverSpread`.
#'
#' @import sp tidyr rgdal raster ggplot2 viridis ggpubr
#'
#' @param landcover_sim_results list. The results of the `LandCoverSpread` function.
#' @param depvar_sim_results list. The results of the `ChangeLandcover_to_ChangeDepVar` function.
#' @param infest_val numerical. The value of the invasive landcover.
#' @param suscep_val numerical. Vector of landcover values that are susceptible to the spread.
#' @param dep_var_modified logical. If `TRUE`, creates line graphs with added modified dependent variable from the simulation. Requires `dep_var_modifier` to be specified in `LandCoverSpread`.
#' @param line_color_labels character vector. If `dep_var_modified == TRUE`, specifies the labels for the dependent variable and the modified dependent variable.
#' @param font_size value. Font size of text in the figure.
#' @param line_thickness value. For the line graph, the thickness of the lines.
#' @param flip_colors logical. Flips the colors of the raster plots
#' @param decimal_places numerical. Number of decimal places to report on the legend of continuous plots.
#' @param positive_vals_only logical. When plotting the line graph, should the negative values be removed when aggregating by year?
#' @param n_grid numerical. If you want a grid of timelapse figures, specify the number of plots.
#'
#' @return A list of landcover plots and line plots associated with the landcover spread simulation.
#'
#' @details Included in the resulting list is a list of yearly landcover plots, a list of yearly dependent variable plots, and a line graph of yearly changes both the dependent variable and, if specified, the modified dependent variable.
#'
#' @examples
#' # load packages
#' library(LandCover); library(foreach)
#'
#' # initialize data.frame with coordinates
#' dat <- expand.grid(x = 1:20, y = 1:20, KEEP.OUT.ATTRS = FALSE)
#'
#' # create some data: elevation, landcover, and temp/ET dependent on elevation and landcover
#' dat$elevation <- with(dat, 50 + 2*x + 5*y + rnorm(nrow(dat), sd = 7))
#' dat$landcover <- ifelse(dat$elevation < median(dat$elevation), 1, 2)
#' dat$temp      <- with(dat, (120-0.7*(0.5*elevation + 0.3*y - 0.5*x + ifelse(landcover == 'lc1', -30, 0) + rnorm(nrow(dat)))))
#' dat$ET        <- with(dat, (   -0.4*(-2*temp       + 0.5*y - 1.0*x + ifelse(landcover == 'lc1', +20, 0) + rnorm(nrow(dat)))))
#'
#' # run the gls model
#' regression_results <- gls_spatial(data = dat, landcover_varname = 'landcover', landcover_vec = c(1,2),
#'                                   reg_formula = ET ~ elevation + temp, error_formula = ~ x + y)
#'
#' # predict values of ET before and after invasion
#' pred_values <- gls_spatial_predict(data = dat, regression_results = regression_results, landcover_varname = 'landcover', landcover_invasive = 1, landcover_susceptible = 2,
#'                                    dep_varname = 'ET', x_coords_varname = 'x', y_coords_varname = 'y')
#'
#' # get landcover raster
#' lc_raster <- rasterFromXYZ(dat[c('x', 'y', 'landcover')])
#'
#' # run landcover simulation
#' landcover_sim <- LandCoverSpread(infest_val = 1, suscep_val = 2, spread_rate = 0.05, birdcell = 0, simlength = 15, simulation_count = 100,
#'                                  lc_raster = lc_raster, dep_var_raster_initial = pred_values$`Predicted values raster, current landcover`,
#'                                  dep_var_raster_pred = pred_values$`Predicted values raster, post-invasion`,
#'                                  dep_var_modifier = 0.80, silent = TRUE)
#'
#' @export


### FUNCTION:

SimulationPlots <- function(landcover_sim_results, depvar_sim_results, infest_val, suscep_val, dep_var_modified = TRUE,
                            dep_var_label = 'dep_var',
                            dep_var_modified_label = 'dep_var_modified',
                            font_size = 15, line_thickness = 0.8,
                            line_colors = viridis(3)[1:2], line_color_labels = c('Dep. var.', 'Modified dep. var.'),
                            line_plot_axis_label = 'Dep. var.',
                            flip_colors = FALSE,
                            decimal_places = 2,
                            infest_label = 'Invasive',
                            suscep_label = 'Susceptible',
                            positive_vals_only = TRUE,
                            n_grid = NA){




  ##### create list of landcover plots over time #####

  # initiate list
  lc_plot_timelapse <- list()

  # fill list
  for(i in 1:length(landcover_sim_results$list_of_landcover_rasters)){

    # create data.frame from raster
    df <- as.data.frame(landcover_sim_results$list_of_landcover_rasters[[i]], xy = TRUE)


    # convert to invasive or non-invasive, then remove cells that are NA
    df$landcover_category <- ifelse(df[,3] %in% infest_val, infest_label,
                             ifelse(df[,3] %in% suscep_val, suscep_label, 'Other'))

    # get colors for landcovers
    lc_colors <- viridis(length(unique(df$landcover_category)))
    if(isTRUE(flip_colors)){ lc_colors = rev(lc_colors) }

    # create landcover color scheme
    lc_groups <- factor(c(infest_label, suscep_label, 'Other'), levels = c(infest_label, suscep_label, 'Other'))
    lc_groups_colors <- c('#FDE725FF', '#440154FF', 'lightgray')

    # plot
    lc_plot_timelapse[[i]] <- ggplot(data = df[!is.na(df[,3]),]) +

      geom_raster(aes(x, y, fill = landcover_category)) +

      scale_fill_manual(limits = lc_groups, values = lc_groups_colors, name = 'Landcover') +

      coord_equal() +

      labs(title = paste0('Year ', i-1)) +

      theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.background = element_blank(), text = element_text(size = font_size))

  }




  ##### create list of dep_var plots over time #####

  # get range of values for limits of color scale
  color_limits <- range(c(sapply(depvar_sim_results$depvar_list_change_from_year_0, function(ras){range(values(ras), na.rm = TRUE)})), na.rm = TRUE)

  # set color scale according to whether the colors are flipped
  if(isFALSE(flip_colors)){

    fill_low  = '#440154FF'
    fill_mid  = 'lightgray'
    fill_high = '#FDE725FF'

  } else {

    fill_low  = '#FDE725FF'
    fill_mid  = 'lightgray'
    fill_high = '#440154FF'

  }

  # initiate list
  dep_var_plot_timelapse <- list()
  df_list                <- list()

  # fill list
  for(i in 1:length(depvar_sim_results$depvar_list_change_from_year_0)){

    # create data.frame from raster
    df_list[[i]] <- as.data.frame(depvar_sim_results$depvar_list_change_from_year_0[[i]], xy = TRUE)

    # create breaks and labels for the plots, while checking to see if any change in dependent variable
    legend_breaks = c(min(color_limits), mean(color_limits), max(color_limits))
    if(length(unique(legend_breaks))==1){ legend_breaks[[1]] = legend_breaks[[1]]-1*10^{-1*decimal_places}
                                          legend_breaks[[3]] = legend_breaks[[3]]+1*10^{-1*decimal_places} }
    legend_labels = format(round(legend_breaks, digits = decimal_places), nsmall = decimal_places)


    # plot
    df <- df_list[[i]]
    colnames(df)[[3]] <- 'layer'

    dep_var_plot_timelapse[[i]] <- ggplot(data = df[!is.na(df$layer),]) +

      geom_raster(aes(x, y, fill = layer)) +

      scale_fill_gradient2(name = dep_var_label, low = fill_low, mid = fill_mid, high = fill_high,
                           limits = range(legend_breaks, na.rm = TRUE),
                           breaks = legend_breaks,
                           labels = ) +

      coord_equal() +

      labs(title = paste0('Year ', i-1), fill = paste0('Change in ', dep_var_label, 'from year 0')) +

      theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.background = element_blank(), text = element_text(size = font_size))

  }




  ##### create list of dep_var_modified plots over time, if specified #####

  if(isTRUE(dep_var_modified)){

    # get range of values for limits of color scale
    color_limits <- range(c(sapply(depvar_sim_results$depvar_list_change_from_year_0_modified, function(ras){range(values(ras), na.rm = TRUE)})), na.rm = TRUE)

    # set color scale according to whether the colors are flipped
    if(isFALSE(flip_colors)){

      fill_low  = '#440154FF'
      fill_mid  = 'lightgray'
      fill_high = '#FDE725FF'

    } else {

      fill_low  = '#FDE725FF'
      fill_mid  = 'lightgray'
      fill_high = '#440154FF'

    }

    # initiate list
    dep_var_modified_plot_timelapse <- list()
    df_list                         <- list()

    # fill list
    for(i in 1:length(depvar_sim_results$depvar_list_change_from_year_0_modified)){

      # create data.frame from raster
      df_list[[i]] <- as.data.frame(depvar_sim_results$depvar_list_change_from_year_0_modified[[i]], xy = TRUE)

      # create breaks and labels for the plots, while checking to see if any change in dependent variable
      legend_breaks = c(min(color_limits), mean(color_limits), max(color_limits))
      if(length(unique(legend_breaks))==1){ legend_breaks[[1]] = legend_breaks[[1]]-1*10^{-1*decimal_places}
      legend_breaks[[3]] = legend_breaks[[3]]+1*10^{-1*decimal_places} }
      legend_labels = format(round(legend_breaks, digits = decimal_places), nsmall = decimal_places)

      # plot
      df <- df_list[[i]]
      colnames(df)[[3]] <- 'layer'

      dep_var_modified_plot_timelapse[[i]] <- ggplot(data = df[!is.na(df$layer),]) +

        geom_raster(aes(x, y, fill = layer)) +

        scale_fill_gradient2(name = dep_var_modified_label, low = fill_low, mid = fill_mid, high = fill_high,
                             limits = range(legend_breaks, na.rm = TRUE),
                             breaks = legend_breaks,
                             labels = legend_labels) +

        coord_equal() +

        labs(title = paste0('Year ', i-1), fill = paste0('Change in ', dep_var_label, 'from year 0')) +

        theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.background = element_blank(), text = element_text(size = font_size))

    }

  }





  ##### create line graphs #####

  # if depvar is modified, create line graph with depvar and modified depvar
  if(isTRUE(dep_var_modified)){

    # get data - sum of raster values depends on whether only positive values are wanted
    if(positive_vals_only){
      LineGraphData <- data.frame(Year             = 0:(length(landcover_sim_results$list_of_landcover_rasters)-1),
                                  dep_var          = sapply(depvar_sim_results$depvar_list_change_from_year_0,          function(r){ sum(values(r)[values(r)>0], na.rm = TRUE) }),
                                  dep_var_modified = sapply(depvar_sim_results$depvar_list_change_from_year_0_modified, function(r){ sum(values(r)[values(r)>0], na.rm = TRUE) }))
    } else {
      LineGraphData <- data.frame(Year             = 0:(length(landcover_sim_results$list_of_landcover_rasters)-1),
                                  dep_var          = sapply(depvar_sim_results$depvar_list_change_from_year_0,          function(r){ sum(values(r), na.rm = TRUE) }),
                                  dep_var_modified = sapply(depvar_sim_results$depvar_list_change_from_year_0_modified, function(r){ sum(values(r), na.rm = TRUE) }))
    }


    # melt data
    LineGraphData_melt <- with(LineGraphData, data.frame(Year     = c(Year, Year),
                                                         variable = rep(line_color_labels, each = nrow(LineGraphData)),
                                                         value    = c(dep_var, dep_var_modified)))


    # plot
    LineGraph <- ggplot(data = LineGraphData_melt) + geom_line(aes(x = Year, y = value, color = variable), size = line_thickness) +

      scale_color_manual(values = line_colors, name = NULL) +

      theme(text = element_text(size = font_size)) +

      labs(y = line_plot_axis_label)


  # if depvar isn't modified, just create line graph of depvar
  } else {

    # adjust colors if they're specified
    if(is.na(line_colors[[1]])){ line_colors = viridis(2)} else { line_colors = line_colors}

    # get data
    LineGraphData <- data.frame(Year             = 0:(length(landcover_sim_results$list_of_landcover_rasters)-1                                                   ),
                                dep_var          = sapply(depvar_sim_results$depvar_list_change_from_year_0,          function(r){ sum(values(r), na.rm = TRUE) }))


    # plot
    LineGraph <- ggplot(data = LineGraphData) + geom_line(aes(x = Year, y = dep_var), size = line_thickness) +

      scale_color_manual(values = line_colors) +

      theme(legend.position = NULL, text = element_text(size = font_size)) +

      labs(y = line_plot_axis_label)

  }




  ##### grid plots #####

  # if n_grid is not NA, plot the lc grid
  if(!is.na(n_grid)){

    # get sequence from 1 to n_grid for plot indicies
    plot_indices <- round(seq(1, length(lc_plot_timelapse), length.out = n_grid), digits = 0)

    # create list of plots
    lc_plot_list <- list()
    for(i in 1:length(plot_indices)){

      lc_plot_list[[i]] <- lc_plot_timelapse[[ plot_indices[[i]] ]]

    }

    # plot grid
    lc_timelapse_grid <- ggpubr::ggarrange(plotlist = lc_plot_list, common.legend = TRUE, legend = 'right')

  }

  # if n_grid is not NA, plot the dep var grid
  if(!is.na(n_grid)){

    # get sequence from 1 to n_grid for plot indicies
    plot_indices <- round(seq(1, length(dep_var_plot_timelapse), length.out = n_grid), digits = 0)

    # create list of plots
    dep_var_plot_list <- list()
    for(i in 1:length(plot_indices)){

      dep_var_plot_list[[i]] <- dep_var_plot_timelapse[[ plot_indices[[i]] ]]

    }

    # plot grid
    dep_var_timelapse_grid <- ggpubr::ggarrange(plotlist = dep_var_plot_list, common.legend = TRUE, legend = 'right')


  }

  # if n_grid is not NA and dep_var_modified is not NA, plot the dep var modified grid
  if(!is.na(n_grid) & isTRUE(dep_var_modified)){

    # get sequence from 1 to n_grid for plot indicies
    plot_indices <- round(seq(1, length(dep_var_modified_plot_timelapse), length.out = n_grid), digits = 0)

    # create list of plots
    dep_var_modified_plot_list <- list()
    for(i in 1:length(plot_indices)){

      dep_var_modified_plot_list[[i]] <- dep_var_modified_plot_timelapse[[ plot_indices[[i]] ]]

    }

    # plot grid
    dep_var_modified_timelapse_grid <- ggpubr::ggarrange(plotlist = dep_var_modified_plot_list, common.legend = TRUE, legend = 'right')


  }





  ##### return results #####

  # return results no modified dep var and no grid
  if(isFALSE(dep_var_modified) & is.na(n_grid)){

    sim_plots_results <- list(lc_plot_timelapse               = lc_plot_timelapse,
                              dep_var_plot_timelapse          = dep_var_plot_timelapse,
                              line_graph                      = LineGraph,
                              line_graph_data                 = LineGraphData)

    return(sim_plots_results)

  }

  # return results no modified dep var but with grid
  if(isFALSE(dep_var_modified) & !is.na(n_grid)){

    sim_plots_results <- list(lc_plot_timelapse               = lc_plot_timelapse,
                              lc_timelapse_grid               = lc_timelapse_grid,
                              dep_var_plot_timelapse          = dep_var_plot_timelapse,
                              dep_var_timelapse_grid          = dep_var_timelapse_grid,
                              line_graph                      = LineGraph,
                              line_graph_data                 = LineGraphData)

    return(sim_plots_results)

  }

  # return results with modified dep_var but no grid
  if(isTRUE(dep_var_modified) & is.na(n_grid)){

    sim_plots_results <- list(lc_plot_timelapse               = lc_plot_timelapse,
                              dep_var_plot_timelapse          = dep_var_plot_timelapse,
                              dep_var_modified_plot_timelapse = dep_var_modified_plot_timelapse,
                              line_graph                      = LineGraph,
                              line_graph_data                 = LineGraphData_melt)

    return(sim_plots_results)

  }

  # return results with modified dep_var and with grid
  if(isTRUE(dep_var_modified) & !is.na(n_grid)){

    sim_plots_results <- list(lc_plot_timelapse               = lc_plot_timelapse,
                              lc_timelapse_grid               = lc_timelapse_grid,
                              dep_var_plot_timelapse          = dep_var_plot_timelapse,
                              dep_var_timelapse_grid          = dep_var_timelapse_grid,
                              dep_var_modified_plot_timelapse = dep_var_modified_plot_timelapse,
                              dep_var_modified_timelapse_grid = dep_var_modified_timelapse_grid,
                              line_graph                      = LineGraph,
                              line_graph_data                 = LineGraphData_melt)

    return(sim_plots_results)

  }

}
natedemaagd/LandCover documentation built on April 1, 2021, 4:14 p.m.