5 New mizer Plotting functions

R/comparedietmatrix.R

@@ -0,0 +1,74 @@
+#'I am not sure if this code works as intended. It appears so, as I have checked it,
+#'but the results are not what it is expected.
+#'
+#'I will also need to see if edge cases, where the timerange is 1:1, work.
+
+
+
+#This function plots the diet matrix from the mizersim objects.
+#' Plot Relative Diet Proportion of each Prey/Predator  
+#'
+#' This function takes two mizerSim objects and calculates the relative 
+#' change in the proportion of a given prey species in a predators diet. This 
+#' is done for every prey/predator in the model. 
+#'
+#' @param harvested A mizerSim object
+#' @param unharvested A mizerSim object - to compare to.
+#' @param chosentime The year range to plot (example 1:2).
+#'
+#' @return A ggplot object of a matrix of predator species on the X axis, 
+#' prey species on the Y axis. The colour of the box indicates the change 
+#' of the proportion in the predator's diet of the given prey species.
+#' 
+#'
+#' @examples
+#' harvested <- getBiomass(NS_sim)
+#' unharvested <- getBiomass(NS_sim)
+#' comparedietmatrix(harvested, unharvested, 5)
+#'
+#' @export
+comparedietmatrix <- function(unharvestedprojection, harvestedprojection, timerange){
+
+  #THE TIMERANGE SHOULD BE X:Y
+  dietunharv <- getDiet(unharvestedprojection@params, 
+                        n = apply(unharvestedprojection@n[timerange,,], c(2, 3), mean),
+                        n_pp = apply(unharvestedprojection@n_pp[timerange,], 2, mean),
+                        n_other = apply(unharvestedprojection@n_other[timerange,], 2, mean),
+                        proportion = TRUE) %>%
+    as.table()%>%
+    as.data.frame()%>%
+    group_by(predator, prey)%>%
+    summarise(Proportion=mean(Freq))
+
+  dietharv <- getDiet(unharvestedprojection@params, 
+                      n = apply(unharvestedprojection@n[timerange,,], c(2, 3), mean),
+                      n_pp = apply(unharvestedprojection@n_pp[timerange,], 2, mean),
+                      n_other = apply(unharvestedprojection@n_other[timerange,], 2, mean),
+                      proportion = TRUE) %>%
+    as.table()%>%
+    as.data.frame()%>%
+    group_by(predator, prey)%>%
+    summarise(Proportion=mean(Freq))
+
+  joindiet <- left_join(dietharv, dietunharv, by = c("prey", "predator"))%>%
+    mutate(Difference = ((Proportion.x - Proportion.y) / Proportion.y) * 100) %>%  # Calculate percentage change
+    select(predator, prey, Difference)%>%
+    filter(!predator %in% c("2", "4", "6", "8", "16", "17", "18", "19", "20", "Resource"), 
+           !prey %in% c("2", "4", "6", "8", "16", "17", "18", "19", "20", "Resource"))
+
+  dietplot <- ggplot(joindiet, aes(x = predator, y = prey, fill = Difference)) +
+    geom_tile() +  
+    scale_fill_gradient2() +  
+    labs(x = "Predator",
+         y = "Prey",
+         fill = "Difference") +  
+    theme_minimal()+
+    theme(axis.text.x = element_text(angle = 45, hjust = 1,size = 14),
+          axis.text.y = element_text(size = 14),
+          axis.title.x = element_text(size = 16),
+          axis.title.y = element_text(size = 16))
+
+
+  return(dietplot)
+
+}

R/guildplot.R

@@ -0,0 +1,151 @@
+
+#' Plot Guild Relative Change Across Timescales  
+#'
+#' This function takes two mizerSim objects and calculates the relative % 
+#' change in each given feeding guilde in the chosen year, short term (1/2 of the 
+#' chosen year) and the long term (2x the chosen year) 
+#' This function requires a dataframe in the environment titled guildparams - this dataframe should have
+#' a column for minw (minimum weight of guild), maxw (maxmimum weight), guild (
+#' the guild for the given weight), and a column for the species (which matches the mizersim species).
+#' The mizerSim objects must also have a tmax of 2 * year2. 
+#' 
+#'
+#' @param harvested A mizerSim object
+#' @param unharvested A mizerSim object - to compare to.
+#' @param year1 The lower year to plot in the range 
+#' @param year2 The higher year to plot in the range
+#'
+#' @return A ggplot object that plots 3 bars per species - in the short, 
+#' chosen and long time - it plots the relative biomass of each feeding guild 
+#' in comparison to the unharvested.
+#' 
+#'
+#' @examples
+#' harvested <- getBiomass(NS_sim)
+#' unharvested <- getBiomass(NS_sim)
+#' guildplot(harvested, unharvested, 1, 2)
+#'
+#' @export
+guildplot <- function(harvestedprojection, unharvestedprojection, year1, year2) {
+
+  harvestedshort <- plotSpectra(harvestedprojection, time_range = max(1, round(year1 * (1/2))):max(1, round(year2 * (1/2))), return_data = TRUE)
+  harvested <- plotSpectra(harvestedprojection, time_range = year1:year2, return_data = TRUE)
+  harvestedlong <- plotSpectra(harvestedprojection, time_range = (year1 * 2):(year2 * 2), return_data = TRUE)
+
+  unharvestedshort <- plotSpectra(unharvestedprojection, time_range = max(1, round(year1 * (1/2))):max(1, round(year2 * (1/2))), return_data = TRUE)
+  unharvested <- plotSpectra(unharvestedprojection, time_range = year1:year2, return_data = TRUE)
+  unharvestedlong <- plotSpectra(unharvestedprojection, time_range = (year1 * 2):(year2 * 2), return_data = TRUE)
+
+  process_guilds <- function(mizerprojection) {
+
+
+    assign_guild <- function(data, rules) {
+      data <- data %>%
+        mutate(Guild = NA_character_)  # Initialize Guild column with NA
+
+      # Loop through each rule in the rules dataframe
+      for (i in 1:nrow(rules)) {
+        data <- data %>%
+          mutate(
+            #THIS CODE ASSUMES THAT ANYTHING UNDER W 0.05 IS PLANKTIVOROUS, AS IT IS VERY SMALL 
+            Guild = ifelse(w < 0.05, "Plank",
+                           ifelse(
+                             is.na(Guild) & w >= rules$minw[i] & w < rules$maxw[i], 
+                             rules$Feeding.guild[i], Guild)
+            )
+          )
+      }
+
+      return(data)
+    }
+
+
+    mizerprojection <- mizerprojection %>%
+      group_by(Species) %>%
+      group_modify(~ {
+        species_data <- .x  
+
+        species_name <- unique(species_data$Legend)
+
+        species_rules <- guildparams %>%
+          filter(Species == species_name)
+
+        if (nrow(species_rules) == 0) {
+          return(species_data)
+        }
+
+        assign_guild(species_data, species_rules)
+
+      }) %>%
+      ungroup() %>%
+      #this next step takes out anything without an assigned guild, but you might not choose to do this
+      #and then you can have a column of the change in biomass of species/sizes that we do not have guild rules for
+      #this would be useful to observe where the biomass change is going, but would be confusing to interpret and explain.
+      #(as its a possibility that all 3 guilds show a negative decrease, which looks like a decrease in biomass,
+      #but may just be due to other sizes/species taking this biomass)
+      drop_na(Guild)%>%
+      group_by(Guild) %>%
+      summarise(value = mean(value))
+
+    return(mizerprojection)
+
+  }
+
+  #for the harvested - 
+  guildsshort <- process_guilds(harvestedshort)
+  guilds <- process_guilds(harvested)
+  guildslong <- process_guilds(harvestedlong)
+  #for the unharvested - 
+  unguildsshort <- process_guilds(unharvestedshort)
+  unguilds <- process_guilds(unharvested)
+  unguildslong <- process_guilds(unharvestedlong)
+
+  #now joining them together
+  guildsshort$time <- "short"
+  guilds$time <- "chosen"
+  guildslong$time <- "long"
+  unguildsshort$time <- "short"
+  unguilds$time <- "chosen"
+  unguildslong$time <- "long"
+
+  joinedguilds <- bind_rows(guildsshort, guilds, guildslong) %>%
+    group_by(Guild, time) %>%
+    summarise(value = sum(value, na.rm = TRUE), .groups = "drop")
+
+  unjoinedguilds <- bind_rows(unguildsshort, unguilds, unguildslong) %>%
+    group_by(Guild, time) %>%
+    summarise(value = sum(value, na.rm = TRUE), .groups = "drop")
+
+  joinedguilds <- joinedguilds%>%
+    full_join(unjoinedguilds, by = c("Guild", "time"),relationship = "many-to-many") %>%
+    mutate(percentage_diff = ((value.x-value.y)/value.y))%>%
+    select(Guild, time, percentage_diff)
+
+  #this sets the colours correctly
+  joinedguilds$time <- factor(joinedguilds$time, levels = c("short", "chosen", "long"))
+  joinedguilds$fill_group <- interaction(joinedguilds$percentage_diff >= 0, joinedguilds$time)
+
+  #plotting
+  ggplot(joinedguilds, aes(x = Guild, y = percentage_diff, fill = fill_group)) +
+    geom_bar(stat = "identity", position = position_dodge(width = 0.9)) +
+    scale_fill_manual(values = c(
+      "FALSE.short" = "#E76F51",  
+      "FALSE.chosen" = "#E98C6B",  
+      "FALSE.long" = "#F2A488",   
+      "TRUE.short" = "#2FA4E7", 
+      "TRUE.chosen" = "#2FA4E7cc",
+      "TRUE.long" = "#2FA4E799" 
+    )) +
+    labs(title = "Percentage Change by Guild", 
+         x = "Guild", 
+         y = "Percentage Change") +
+    theme_minimal() +
+    theme(
+      axis.text.x = element_text(size = 14, angle = 90, hjust = 1, vjust = 0.5),
+      axis.text.y = element_text(size = 14),
+      legend.position = "none",
+      axis.title.x = element_text(size = 16),
+      axis.title.y = element_text(size = 16)
+    )
+
+}

R/plotSpeciesWithTimeRange.R

@@ -0,0 +1,142 @@
+#This function plots the species plot - which the change in species for a given 
+#year, and also for 2x in future and 1/2 year in the past.
+
+
+#' Plot MizerSim Relative Biomass per Species Across Varying Timescales 
+#'
+#' This function takes two mizerSim objects and calculates the relative % 
+#' change in each given species in the chosen year, short term (1/2 of the 
+#' chosen year) and the long term (2x the chosen year). The mizerSim
+#' objects must have tmax = 2 * chosentime2. 
+#'
+#' @param harvested A mizerSim object
+#' @param unharvested A mizerSim object - to compare to.
+#' @param chosentime The year to plot 
+#'
+#' @return A ggplot object that plots 3 bars per species - in the short, 
+#' chosen and long time - it plots the relative biomass of each species in
+#' comparison to the unharvested.
+#' 
+#'
+#' @examples
+#' harvested <- getBiomass(NS_sim)
+#' unharvested <- getBiomass(NS_sim)
+#' plotSpeciesWithTimeRange(harvested, unharvested, 1, 2)
+#'
+#' @export
+plotSpeciesWithTimeRange <- function(harvestedprojection, unharvestedprojection, chosentime1, chosentime2) {
+
+  #get the biomass of the species
+  unharvestedbio <- getBiomass(unharvestedprojection) %>%
+    .[chosentime1:chosentime2, ] %>% 
+    melt() %>%
+    group_by(sp) %>%
+    summarize(value = mean(value, na.rm = TRUE))
+
+  harvestedbio <- getBiomass(harvestedprojection) %>%
+    .[chosentime1:chosentime2, ] %>% 
+    melt() %>%
+    group_by(sp) %>%
+    summarize(value = mean(value, na.rm = TRUE))
+
+  #calculate percentage change in species in the chosen year
+  percentage_diff <-  harvestedbio %>%
+    left_join(unharvestedbio, by = "sp") %>%
+    mutate(percentage_diff = ((value.x - value.y) / value.y) * 100,
+           Species = sp) %>%
+    select(Species, percentage_diff) %>%
+    filter(!Species %in% c("2", "4", "6", "8", "16", "17", "18", "19", "20", "Resource"))%>%
+    mutate(class = "chosen")
+
+  calculate_biomass_triples <- function(unharvestedprojection, harvestedprojection, year1, year2) {
+
+    # Calculate unharvested biomass at different time points
+    unharvestedbiotriple <- getBiomass(unharvestedprojection)
+
+    #the range has to be 1-2, becuase if it is 1-1 it messes with the way the data is formatted.
+    #i have also used ceiling here, because using round means they round to nearest even number, so some cases (11:13)
+    #end up as 6:6 for the lowbiotrip - this does not work for the code format.
+    lowunbiotrip <- unharvestedbiotriple[max(1, ceiling(year1 * (1/2))):max(2, ceiling(year2 * (1/2))), ] %>%
+      melt() %>%
+      group_by(sp) %>%
+      summarize(value = mean(value, na.rm = TRUE))
+
+    highunbiotrip <- unharvestedbiotriple[(year1 * 2):(year2 * 2), ] %>%
+      melt() %>%
+      group_by(sp) %>%
+      summarize(value = mean(value, na.rm = TRUE))
+
+    # Calculate harvested biomass at different time points
+    harvestedbiotriple <- getBiomass(harvestedprojection)
+
+    lowbiotrip <- harvestedbiotriple[max(1, ceiling(year1 * (1/2))):max(2, ceiling(year2 * (1/2))),] %>%
+      melt() %>%
+      group_by(sp) %>%
+      summarize(value = mean(value, na.rm = TRUE))
+
+    highbiotrip <- harvestedbiotriple[(year1 * 2):(year2 * 2), ] %>%
+      melt() %>%
+      group_by(sp) %>%
+      summarize(value = mean(value, na.rm = TRUE))
+
+    # Return the results as a list
+    list(
+      lowunbiotrip,
+      highunbiotrip,
+      lowbiotrip,
+      highbiotrip
+    )
+  }
+
+  #calculate percentage change in other years
+  biorange <- calculate_biomass_triples(unharvestedprojection, harvestedprojection, chosentime1, chosentime2)
+
+  #percentage_difflow <- percentdiff(biorange[[3]], biorange[[1]])
+  #percentage_difflow$class <- "short"
+
+  percentage_difflow <-  biorange[[3]] %>%
+    left_join(biorange[[1]], by = "sp") %>%
+    mutate(percentage_diff = ((value.x - value.y) / value.y) * 100,
+           Species = sp) %>%
+    select(Species, percentage_diff) %>%
+    filter(!Species %in% c("2", "4", "6", "8", "16", "17", "18", "19", "20", "Resource"))%>%
+    mutate(class = "short")
+
+  #percentage_diffhigh <- percentdiff(biorange[[4]], biorange[[2]])
+  #percentage_diffhigh$class <- "long"
+
+  percentage_diffhigh <-  biorange[[4]] %>%
+    left_join(biorange[[2]], by = "sp") %>%
+    mutate(percentage_diff = ((value.x - value.y) / value.y) * 100,
+           Species = sp) %>%
+    select(Species, percentage_diff) %>%
+    filter(!Species %in% c("2", "4", "6", "8", "16", "17", "18", "19", "20", "Resource"))%>%
+    mutate(class = "long")
+
+  percentage_diff <- rbind(percentage_difflow, percentage_diff, percentage_diffhigh)
+
+  #now plot them together - the first lines sort out the colors of the bars
+  percentage_diff$class <- factor(percentage_diff$class, levels = c("short", "chosen", "long"))
+  percentage_diff$fill_group <- interaction(percentage_diff$percentage_diff >= 0, percentage_diff$class)
+
+  ggplot(percentage_diff, aes(x = Species, y = percentage_diff, fill = fill_group)) +
+    geom_bar(stat = "identity", position = position_dodge(width = 0.9)) +
+    geom_hline(yintercept = 0, color = "grey", linetype = "dashed", size = 0.5)+
+    labs(x = "Species", y = "Percentage Change") +
+    scale_fill_manual(values = c(
+      "FALSE.short" = "#E76F51",  
+      "FALSE.chosen" = "#E98C6B",  
+      "FALSE.long" = "#F2A488",   
+      "TRUE.short" = "#2FA4E7", 
+      "TRUE.chosen" = "#2FA4E7cc",
+      "TRUE.long" = "#2FA4E799" 
+    )) +
+    theme_minimal() +
+    theme(
+      axis.text.x = element_text(size = 16, angle = 90, hjust = 1, vjust = 0.5),
+      axis.text.y = element_text(size = 14),
+      legend.position = "none",
+      axis.title.x = element_text(size = 16),
+      axis.title.y = element_text(size = 16)
+    )
+}