figures.R

# Read, Q. D., et al. 2019. Assessing the environmental impacts of halving food loss and waste along the food supply chain. Science of the Total Environment.   #
#                                                                                                                                                               #
# Code to reproduce all data figures in manuscript (Figures 2, 3, and 4), beginning with output generated by analysis.R                                         #
# Last modified by QDR on 07 January 2020                                                                                                                       #
# Tested with R version 3.6.2 on Windows and Linux                                                                                                              #                                                                                                   
# Contact: qread@sesync.org                                                                                                                                     #
#################################################################################################################################################################

# ===========================================================
# Load all packaages and data to create Figures 2, 3, and 4 =
# ===========================================================

library(tidyverse)
library(units)
library(gridExtra)
library(ggrepel)
library(cowplot)

grid_result <- read.csv('output/sixstage_scenario_grid_lcia_results.csv', stringsAsFactors = FALSE)
base_and_nowaste <- grid_result %>%
  filter(rowSums(.[,1:6]) %in% c(0, 6))
grid_uncertainty_CIs <- read.csv('output/uncertainty_grid_CIs.csv', stringsAsFactors = FALSE)
results_by_commodity_df <- read.csv('output/bestpathway_bycommodity.csv', stringsAsFactors = FALSE)

############
# Figure 2 #
############

category_labels_expression <- list("Eutrophication potential (kg N eq.)",
                                   expression(paste("GHG emissions (kg ", CO[2], " eq.)")),
                                   "Energy use (GJ)",
                                   expression(paste("Land used (", m^2, ")")),
                                   expression(paste("Water used ", (m^3))))

waste_impacts <- base_and_nowaste %>%
  mutate(scenario = if_else(rowSums(.[,1:6]) == 0, 'baseline', 'zero_waste')) %>%
  select(scenario, impact_category, value) %>%
  spread(scenario, value) %>%
  mutate(waste_impact = baseline - zero_waste,
         proportion_impact = waste_impact / baseline)

# US population in 2012 was 314 million
pop <- 314e6

waste_impacts <- mutate(waste_impacts, percapita_waste_impact = waste_impact / pop)

# Hardcoded literature values for comparison.
energy_litvalues <- c(1740, 2030, 2559)
energy_mj <- set_units(set_units(energy_litvalues*1e12, 'BTU'), 'MJ') / 314e6
          
values <- c(405,498,1117,2.7,42,53,54,352,energy_mj,329)

lit_df <- data.frame(impact_category = c('land','land','land','eutrophication','water','water','water','water','energy','energy','energy','GHG'),
                     source = 'literature',
                     percapita_waste_impact = values)

waste_impacts <- waste_impacts %>% 
  filter(grepl('land|watr|enrg|eutr|co2', impact_category)) %>%
  mutate(impact_category = c('eutrophication','GHG','energy','land','water'),
         source = 'this study') %>%
  select(impact_category, source, percapita_waste_impact) %>%
  rbind(lit_df) %>%
  mutate(percapita_waste_impact = if_else(impact_category=='energy',percapita_waste_impact/1000, percapita_waste_impact))
  
plotargs <- tibble(impact_category = unique(waste_impacts$impact_category),
                       maxy = c(3, 600, 10, 2000, 80),
                       labelexpr = category_labels_expression)

waste_plots <- waste_impacts %>%
  left_join(plotargs) %>%
  group_by(impact_category) %>%
  do(p = ggplot(., aes(x = 1, y = percapita_waste_impact, fill = source, alpha = source, size = source)) +
       geom_point(pch = 21) +
       geom_hline(yintercept = 0, color = 'white') +
       scale_x_continuous(limits=c(0.95,1.05), expand=c(0,0)) +
       scale_y_continuous(limits=c(0,.$maxy[1]), expand=c(0,0), name = .$labelexpr[[1]]) +
       scale_alpha_manual(values = c(0.3, 1)) +
       scale_size_manual(values = c(3, 5)) +
       scale_fill_manual(values = c('black', 'red')) +
       theme_classic() +
       theme(axis.ticks.x = element_blank(), axis.text.x = element_blank(), axis.title.x = element_blank(), axis.line.x = element_blank(), legend.position = 'none'))

# Change ticks on plot 1
waste_plots$p[[1]] <- waste_plots$p[[1]] + scale_y_continuous(limits=c(0, 10), expand=c(0,0), name = 'Energy use (GJ)', breaks = c(0,5,10))

# Draw figure 2
grid.arrange(grobs = waste_plots$p, nrow = 1)

############
# Figure 3 #
############

stage_full_names <- c('production', 'processing', 'retail', 'consumption: foodservice', 'consumption: institutional', 'consumption: household')
stage_full_names_lookup <- c(none = '', L1 = 'production', L2 = 'processing', L3 = 'retail', L4a = 'consumption:\nfoodservice', L4b = 'consumption:\ninstitutional', L5 = 'consumption:\nhousehold')

# Calculate sequences for each iteration of uncertainty analysis.
create_sequence <- function(dat, value_id, proportion) {
  value_id <- enquo(value_id)
  dat %>%
    ungroup %>%
    mutate(nbypct = rowSums(.[,1:6] == proportion)) %>%
    filter(rowSums(.[,1:6] > 0) == nbypct) %>%
    group_by(nbypct) %>%
    filter(!!value_id == min(!!value_id)) %>%
    arrange(nbypct) %>%
    ungroup %>%
    mutate(stage_reduced = c('none', names(sort(apply(.[,1:6],2,function(x) which(diff(x) > 0))))))
}

# Find the sequence for the actual parameter values.
trueseq <- grid_result %>%
  filter(grepl('co2|watr|land|enrg|eutr', impact_category)) %>%
  group_by(impact_category) %>%
  do(create_sequence(., value, 0.5))

# Match the true sequence values with the error bar values.
trueseq_withcis <- trueseq %>% left_join(grid_uncertainty_CIs)

# Function to create reduction plot with error bars
reduction_plot_with_errorbars <- function(sequence, yval, yminval, ymaxval, plot_title, plot_subtitle) {
  yval <- enquo(yval)
  yminval <- enquo(yminval)
  ymaxval <- enquo(ymaxval)
  sequence %>%
    mutate(norm = max(!!yval)) %>%
    mutate(!!yval := !!yval/norm, 
           !!yminval := !!yminval/norm,
           !!ymaxval := !!ymaxval/norm) %>%
    mutate(stage_reduced = stage_full_names_lookup[stage_reduced]) %>%
    ggplot(aes(x = nbypct, y = !!yval, ymin = !!yminval, ymax = !!ymaxval)) +
    geom_line(size = 1) +
    geom_errorbar(width = 0.2, color = 'gray60') +
    geom_point(size = 3, color = 'white', fill = 'black', shape = 21, stroke = 2) +
    geom_text(aes(label = stage_reduced), angle = 45, hjust = 1.25) +
    scale_x_continuous(name = 'Number of stages where waste is reduced', breaks = 0:6) +
    scale_y_continuous(name = 'Impact relative to baseline', labels = scales::percent_format(accuracy=1)) +
    coord_cartesian(ylim = c(0.77, 1.02), expand = TRUE) +
    ggtitle(plot_title, plot_subtitle) +
    theme_bw() +
    theme(panel.grid.minor.x = element_blank(),
          panel.grid.minor.y = element_blank(),
          panel.grid.major.x = element_blank(),
          plot.subtitle = element_text(size = 10))
  
}

# Create data frame with data to draw faceted plot.
facetplotdat <- trueseq_withcis %>%
  filter(grepl('co2|eutr|land|watr|enrg', impact_category)) %>%
  group_by(impact_category) %>%
  mutate(value_norm = value/max(value),
         q025_norm = q025/max(value),
         q975_norm = q975/max(value)) %>%
  ungroup %>%
  mutate(stage_reduced = stage_full_names_lookup[stage_reduced],
         impact_category = factor(impact_category, labels = c('energy use','eutrophication potential', 'greenhouse warming potential', 'land use', 'water use'))) %>%
  mutate(impact_category = factor(impact_category, levels = levels(impact_category)[c(3,1,2,4,5)]))
  
# Draw Figure 3.
ggplot(facetplotdat %>% mutate(letter = letters[1:5][impact_category]), aes(x = nbypct, y = value_norm, ymin = q025_norm, ymax = q975_norm)) +
  geom_line(size = 1) +
  geom_errorbar(width = 0.2, color = 'gray60') +
  geom_point(size = 3, color = 'white', fill = 'black', shape = 21, stroke = 2) +
  geom_text(aes(label = stage_reduced), angle = 45, hjust = 1.25, size = 2.4) +
  geom_text(aes(label = letter), x = 0, y = 0.79, size = 10) +
  scale_x_continuous(name = 'Number of stages where waste is reduced', breaks = 0:6) +
  scale_y_continuous(name = 'Impact relative to baseline', labels = scales::percent_format(accuracy = 1L)) +
  facet_wrap(~ impact_category, nrow = 3) +
  coord_cartesian(ylim = c(0.77, 1.02), expand = TRUE) +
  theme_bw() +
  theme(panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        plot.subtitle = element_text(size = 10),
        strip.background = element_blank())

############
# Figure 4 #
############

# Map stage, category, and food commodity group labels to more descriptive labels
stage_full_names_lookup <- c(none = '', L1 = 'production', L2 = 'processing', L3 = 'retail', L4a = 'foodservice', L4b = 'institutions', L5 = 'households')
categories <- c("impact potential/gcc/kg co2 eq", "resource use/land/m2*yr", "resource use/watr/m3", "resource use/enrg/mj", "impact potential/eutr/kg n eq")
category_short <- c('GHG', 'land', 'water', 'energy', 'eutrophication')
food_categories <- c("cereals", "roots_tubers_fresh", "roots_tubers_processed", 
                     "oilseeds_pulses", "fruit_veg_fresh", "fruit_veg_processed", 
                     "meat", "fish_fresh", "fish_processed", "milk", "eggs", "sugar", 
                     "beverages")
food_categories_labels <- c('cereals', 'roots and tubers (fresh)', 'roots and tubers (processed)', 'oilseeds and pulses', 'fruits and vegetables (fresh)', 'fruits and vegetables (processed)', 'meat', 'fish (fresh)', 'fish (processed)', 'milk', 'eggs', 'sugar', 'beverages')

results_by_commodity_df <- results_by_commodity_df %>%
  mutate(category = category_short[match(category, categories)],
         stage_reduced = stage_full_names_lookup[stage_reduced],
         commodity = food_categories_labels[match(commodity, food_categories)]) 

# Add baseline case for plotting
baseline_df <- data.frame(category = 'all', commodity = unique(results_by_commodity_df$commodity), n_stages_reduced = 0, stage_reduced = 'none', impact = 1, baseline_impact = 1, impact_norm = 1, stringsAsFactors = FALSE)

results_by_commodity_df <- rbind(results_by_commodity_df, baseline_df)

# Function to pull the top 5 food commodity groups
get_top_n <- function(data, the_category, n_comm = 5) {
  data %>% 
    filter(category %in% the_category) %>%
    group_by(commodity) %>%
    filter(n_stages_reduced == max(n_stages_reduced)) %>%
    ungroup %>%
    arrange(impact_norm) %>%
    slice(1:n_comm) %>%
    pull(commodity)
}

# Figures showing top 5 as impact averted ---------------------------------

results_by_commodity_df <- results_by_commodity_df %>% 
  mutate(impact_averted = baseline_impact - impact)

# Convert axis values to a more legible unit by dividing by the 2012 USA population to get per capita impact
pop2012 <- 314e6
percapita_results <- results_by_commodity_df %>%
  mutate_at(vars(impact, baseline_impact, impact_averted), ~ if_else(category%in% 'all', ., .x/pop2012))

# Get all top 5 names to assign colors to them
alltop5 <- map(category_short, ~ get_top_n(results_by_commodity_df, .x, 5)) %>% reduce(c) %>% unique
alltop5_colormap <- setNames(RColorBrewer::brewer.pal(length(alltop5), 'Dark2'), alltop5)

plot_impactaverted <- function(dat, the_category, unit_name, n = 5) {
  dat <- dat %>% 
    filter(commodity %in% get_top_n(., the_category, n), category %in% c('all', the_category))
  
  ggplot(dat, aes(x = n_stages_reduced, y = impact_averted, group = commodity)) + 
    geom_line(size = 1, alpha = 0.33) +
    geom_point(aes(color = commodity), size = 3) +
    geom_text_repel(data = dat %>% filter(n_stages_reduced > 0), aes(label = stage_reduced), size = 2) + # Label to see which stages are being averted
    scale_x_continuous(name = 'Number of stages where waste is reduced', breaks = 0:6) +
    scale_y_continuous(name = parse(text = paste('Per~capita~reduction', unit_name, sep = '~')), 
                       limits = c(0, max(dat$impact_averted) * 1.05), expand = c(0, 0),
                       sec.axis = sec_axis(trans = ~ ./dat$baseline_impact[1], name = 'relative to baseline', labels = scales::percent)) +
    scale_color_manual(name = 'Food group', values = alltop5_colormap, guide = guide_legend(nrow = 2)) +
    theme_bw() +
    theme(panel.grid.minor.x = element_blank(),
          panel.grid.minor.y = element_blank(),
          panel.grid.major.x = element_blank(),
          legend.position = 'bottom')
}

lsize = 10

p_ghg <- plot_impactaverted(percapita_results, 'GHG', '(kg~CO^2~eq.~y^-1)') + annotate('text',x=-Inf,y=Inf,hjust=-0.2,vjust=1.2,label='c', size = lsize) + ggtitle('greenhouse warming potential')
p_land <- plot_impactaverted(percapita_results, 'land', '(m^2~y^-1)') + annotate('text',x=-Inf,y=Inf,hjust=-0.2,vjust=1.2,label='d', size = lsize) + ggtitle('land use')
p_water <- plot_impactaverted(percapita_results, 'water', '(m^3~y^-1)') + annotate('text',x=-Inf,y=Inf,hjust=-0.2,vjust=1.2,label='e', size = lsize) + ggtitle('water use')
p_energy <- plot_impactaverted(percapita_results, 'energy', '(MJ~y^-1)') + annotate('text',x=-Inf,y=Inf,hjust=-0.2,vjust=1.2,label='a', size = lsize) + ggtitle('energy use')
p_eutr <- plot_impactaverted(percapita_results, 'eutrophication', '(kg~N~eq.~y^-1)') + annotate('text',x=-Inf,y=Inf,hjust=-0.2,vjust=1.2,label='b', size = lsize) + ggtitle('eutrophication potential')

# Make a legend for all the plots.
leg_plot <- ggplot(percapita_results %>% filter(commodity %in% alltop5), aes(x=1, y=1, color=commodity)) +
  geom_point() +
  theme_bw() +
  scale_color_manual(name = 'Food group', values = alltop5_colormap, guide = guide_legend(nrow = 4))

leg <- get_legend(leg_plot)

# Lay out the 5 plots plus legend as a 3 x 2 grid.

th_top <- theme(legend.position = 'none',
                axis.text.x = element_blank(),
                axis.title.x = element_blank(),
                axis.ticks.x = element_blank(), 
                plot.title = element_text(hjust = 0.5))
th_bottom <- theme(legend.position = 'none',
                   plot.title = element_text(hjust = 0.5))

top_row <- plot_grid(p_energy + th_top, p_eutr + th_top, nrow = 1)
middle_row <- plot_grid(p_ghg + th_top, p_land + th_bottom, nrow = 1, align = 'h')
bottom_row <- plot_grid(p_water + th_bottom, leg, nrow = 1)

# Draw Figure 4
plot_grid(top_row, middle_row, bottom_row, nrow = 3)