potential

This page synthesizes results on the agronomic potential of N in human excreta in Europe. We compare the total amount of N in human excreta to typical N inputs in agroecosystems: synthetic fertilizers, manure, total N inputs…

knitr::opts_chunk$set(message=F, warning=F, fig.align = "center",  dev='svg')

# Load the function file
source("functions.R")
#included:
#-function to save csv f_save_csv_files
#-set the default scale_color and scale_fill to viridis theme
#-loads the core tidyverse package

#set theme for graphs
theme_set(
  theme_classic() +
  theme(
    panel.grid.major.y = element_line(), #no vertical lines by default
    #text = element_text(family = "Times New Roman"), #default font
    plot.title = element_text(face="bold"), #graphs titles in bolds
    )
  )


#additional packages loaded
library(ggrepel) #for geom_text_repel
library(patchwork) #for multiple plots using +

Source <- "based on FAOSTAT data for 2018-2022"


percent_food_waste_by_region <- tibble(
  world_region = c(
    "North America and Oceania", 
    "Industrialized Asia",
    "Europe", 
    "North Africa, West and Central Asia", 
    "Latin America", 
    "South and Southeat Asia", 
    "Sub-Saharan Africa"
    ),
  #percent_food_waste = c(42, 25, 22, 19, 15, 17, 23)
  percent_food_waste = c(
    round(42*(61+7)/10^4, 2),
    round(25*(46+11)/10^4, 2),
    round(22*(52+9)/10^4, 2),
    round(19*(34+18)/10^4, 2),
    round(15*(28+17)/10^4, 2),
    round(17*(13+15)/10^4, 2),
    round(23*(5+13)/10^4, 2)
    )
)

percent_food_waste_by_region

# A tibble: 7 × 2
  world_region                        percent_food_waste
  <chr>                                            <dbl>
1 North America and Oceania                         0.29
2 Industrialized Asia                               0.14
3 Europe                                            0.13
4 North Africa, West and Central Asia               0.1 
5 Latin America                                     0.07
6 South and Southeat Asia                           0.05
7 Sub-Saharan Africa                                0.04

#percent_food_waste <- 0.25
#percent_food_waste_max <- 0.3
#percent_food_waste_min <- 0.04
percent_N_urine <- 0.85

resolution_dpi <- 350
grahps_height <- 8
grahps_width <- 8

#function to select some countries in our dataset
f_select_countries <- function(dataset, countries_to_select){
  dataset <- dataset %>%
    filter(country %in% countries_to_select)
  return(dataset)
}

Load and prepare data

Load data

Load the cleaned FAO data (see in Data source presentation -> FAO data).

#load data (national for OneEarth, national and regional for STOTEN)
file_FAO_countries <- read_csv("output/fao_manure_fertilizers_balance_food_supply_countries.csv")
file_FAO_world <- read_csv("output/fao_manure_fertilizers_balance_food_supply_world.csv")
file_countries_world_regions <- read_csv("source/countries_by_world_regions.csv")


#simplify name and remove some entities (mainly for our large countries plots)
FAO_countries <- file_FAO_countries %>%
  # Merge world region information into FAO_countries
  left_join(file_countries_world_regions, by = "country") %>%
  # Merge the percent food waste by region
  left_join(percent_food_waste_by_region, by = "world_region") %>%
  #remove particular entities
  filter(
    !country %in% c("China") #this entity gathers mainland China (then simply renamed "China") and Taiwan
    ) %>%
  #simplify some countries names
  mutate(
    country = case_when(
      country == "United Kingdom of Great Britain and Northern Ireland" ~ "UK",
      country == "United States of America" ~ "USA",
      country == "Republic of Korea" ~ "Korea",
      country == "Iran (Islamic Republic of)" ~ "Iran",
      country == "United Republic of Tanzania" ~ "Tanzania",
      country == "Democratic Republic of the Congo" ~ "Dem. Rep. Congo",
      country == "Russian Federation" ~ "Russia",
      country == "China, Taiwan Province of" ~ "Taiwan",
      country == "China, mainland" ~ "China",
      country == "Netherlands (Kingdom of the)" ~ "Netherlands",
      T~country
      )
    )

# assign a global food waste and losses, based on year 2018
temp <- FAO_countries %>% 
  filter(
    Year == 2020, 
    #remove empty values of population and food waste
    is.na(percent_food_waste)==F,
    is.na(`Population (Mhab)`)==F
    )
# check that population summed over all countire is coherent (7.6-8 billion globally on the 2018-2022 period)
sum(temp$`Population (Mhab)`)/1000 # 7.7

[1] 7.75505

# copute the population weighted food waste and losses
percent_food_waste_global <-  weighted.mean(temp$percent_food_waste, temp$`Population (Mhab)`, na.rm = T)
percent_food_waste_global # about 10%

[1] 0.09415807

# add global food waste for the world
file_FAO_world <- file_FAO_world %>%
  mutate(percent_food_waste = percent_food_waste_global)

New variables analyzed

Compute additional variables for analysis.

#function to compute some new variables (excretions or urine as share of N inputs, fertilizers or manue)
f_compute_new_variables <- function(dataset_FAO){
  dataset_FAO <- dataset_FAO %>%
  # compute new variables of interest
    mutate(
  
      # food supply and excretion after food waste
      `human excreta (ktN)` = `N supply (kt/y)`*(1-percent_food_waste),
      `human urine (ktN)` = `human excreta (ktN)`*percent_N_urine,
      
      #   % fertilization (synthetic fertilizers, manure, total)
      #fertilizers
      `human excreta % of synthetic fertilization` = `human excreta (ktN)`/`N mineral fertilization (ktN)`,
      `human urine % of synthetic fertilization` = `human urine (ktN)`/`N mineral fertilization (ktN)`,
      #manure
      `human excreta % of manure applied to soils` = `human excreta (ktN)`/`All Animals Manure applied to soils (ktN)`,
      `human excreta % of ruminants applied to soils` = `human excreta (ktN)`/`ruminants Manure applied to soils (ktN)`,
      `human excreta % of non ruminants applied to soils` = `human excreta (ktN)`/`non ruminants Manure applied to soils (ktN)`,
      #total inputs
      `human excreta % of croplands inputs` = `human excreta (ktN)`/`N croplands Input (ktN)`,
      
      # % livestock excretion
      `human excreta % of livestock excretion` = `human excreta (ktN)`/`All Animals Amount excreted in manure (ktN)`,
      `human excreta % of ruminant excretion` = `human excreta (ktN)`/`ruminants Amount excreted in manure (ktN)`,
      `human excreta % of monogastric excretion` = `human excreta (ktN)`/`non ruminants Amount excreted in manure (ktN)`
      )
  
  return(dataset_FAO)
}

#function to synthesize 2018-2022 mean
f_synthesize_2018_2022 <- function(dataset_FAO){
  
  dataset_FAO <- dataset_FAO %>%
    filter(Year>=2018 & Year<=2022) %>%
    select(-any_of(c("Year", "world_region"))) %>%
    group_by(country) %>%
    summarise_all(~mean(., na.rm=T))
  
  return(dataset_FAO)
}

#compute new variable for countries and world files
FAO_countries <- f_compute_new_variables(FAO_countries)
FAO_world <- f_compute_new_variables(file_FAO_world)

#synthesize for countries and world files
FAO_countries_2018_2022 <- f_synthesize_2018_2022(FAO_countries)
FAO_world_2018_2022 <- f_synthesize_2018_2022(FAO_world)

#get large countrie
countries_to_select_large <- FAO_countries_2018_2022 %>%
  select(country, `Population (Mhab)`) %>%
  distinct() %>%
  filter(`Population (Mhab)`>40) %>%
  pull(country)

#function to plot national bar graphs 
f_graph_countries_bar <- function(dataset, variable){
  ggplot(dataset) +
    
    #countries bars and print % 
    geom_col(
      aes(
        reorder(country, {{ variable }}), 
        {{ variable }}
        )
      ) +
    geom_text(
      data = dataset,
      aes(
        x=country, y={{ variable }}, 
        label=paste0(round({{ variable }}*100), "%")
        ),
      vjust=-0.5, fontface="bold.italic"
      ) +
    
    #world ba and print % 
    annotate(
      geom="text", x="World", y= FAO_world_2018_2022 %>% pull({{ variable }}), 
      label = paste0(round(FAO_world_2018_2022 %>% pull({{ variable }})*100), "%"),
      vjust=-0.5, fontface="bold.italic"
      ) +
    geom_col(
      data = FAO_world_2018_2022,
      aes(country, {{ variable }}),
      fill="red"
      ) +
    
    # rest of graph options
    scale_y_continuous(
      limits = c(0, NA),
      labels = scales::percent,
      breaks = seq(0, 1.5, by=0.1),
      minor_breaks = seq(.1, 1.5, by=0.05)
      ) +
    theme(
      legend.position = "none",
      panel.grid.minor.y= element_line(linewidth=.3),
      axis.text.x = element_text(angle = 45, hjust=1)
      ) +
    labs(
      fill="large country >35M habitants",
      x="",
      y="human excreta as % of\nmineral N fertilizers consumption"
    ) +
    guides(size="none")
}

#function to add extended countried to bar graph
f_graph_countries_bar_extended_countries <- function(dataset, variable){
    #countries bars and print % 
    gg <- geom_col(
      data = dataset,
      aes(
        reorder(country, {{ variable }}), 
        {{ variable }}
        ),
      fill="#440154"
      ) 
    
    return(gg)
}

#prepare inputs file
f_prepare_inputs <- function(dataset){
  dataset <- dataset %>%select(
    #ountry and population
    country, `Population (Mhab)`, 
    
    #total excretions and urine,
    `human excreta (ktN)`, `human urine (ktN)`, percent_food_waste,
    `human excreta % of croplands inputs`, `human urine % of synthetic fertilization`, 
    
    #total fertilizers used
    `N mineral fertilization (ktN)`,
    
    #inputs to croplands
    `N croplands Input (ktN)`,
    `N croplands Atmospheric Deposition (ktN)`, `N croplands Biological Fixation (ktN)`,
    `N croplands Manure applied to Soils (ktN)`, `N croplands Mineral fertilizers (ktN)`,
    
    ) %>%
  
  gather(
    input_flow, value, 
    c(`N croplands Atmospheric Deposition (ktN)`, 
      `N croplands Biological Fixation (ktN)`, 
      `N croplands Manure applied to Soils (ktN)`, 
      `N croplands Mineral fertilizers (ktN)`)
    ) %>%
  
  mutate(
    input_flow = gsub("N croplands ", "", input_flow),
    input_flow = gsub(" \\(ktN\\)", "", input_flow)
    ) 
  
  return(dataset)
}
temp_croplands_inputs <- f_prepare_inputs(FAO_countries_2018_2022)
temp_croplands_inputs_world <- f_prepare_inputs(FAO_world_2018_2022)


f_graph_comparison_inputs_per_cap <- function(dataset, variable_inputs, variable_order, variable_excretions){
  data_excretions <- dataset %>%
    select(country, `Population (Mhab)`, {{ variable_excretions }}, percent_food_waste) %>%
    distinct()
  
  #graph comparison absolute per cap : excretions vs inputs
  gg <- ggplot(dataset) +
    
    #bars columns inputs N in different colors
    geom_col(
      aes(
        {{ variable_inputs }}/`Population (Mhab)`, reorder(country, {{ variable_order }}/`Population (Mhab)`),
        fill=input_flow
        ),
      alpha=.85
      ) +
    
    #error bars of excretions, based on min and max values of excretions
    # geom_linerange(
    #   data = data_excretions,
    #   aes(
    #     xmin={{ variable_excretions }}/`Population (Mhab)`*(1-percent_food_waste_min)/(1-percent_food_waste),
    #     xmax={{ variable_excretions }}/`Population (Mhab)`*(1-percent_food_waste_max)/(1-percent_food_waste),
    #     y=country,
    #     color="human excreta"
    #     )
    #   ) +
    
    #points excretions N in red
    geom_point(
      data = data_excretions,
      aes({{ variable_excretions }}/`Population (Mhab)`, country, color="human excreta")
      ) +
    scale_color_manual(values="red") +
    
    #other options
    scale_x_reverse(
      breaks = seq(0, 1000, 10),
      minor_breaks = seq(0, 1000, 5)
      ) +
    theme(
      axis.text.y = element_text(hjust=0.5),
      legend.position = c(0.35, 0.35),
      legend.background = element_rect(
        fill   = "transparent",  
        colour = NA),
      axis.ticks.y = element_blank(),
      axis.line.y = element_blank(),
      panel.grid.major.y= element_blank(),
      panel.grid.major.x= element_line(),
      panel.grid.minor.x= element_line(),
      ) +
    labs(
      y="", x=expression("kgN" %.% capita^-1 %.% year^-1), color="", fill=""
      ) 
  
  return(gg)
  
}

f_graph_comparison_inputs_per_cap_add_world <- function(graph_countries, dataset_world, variable_inputs, variable_excretions){
    
    #columns
    gg <- graph_countries +
    
    #bars columns inputs N in different colors
    geom_col(
      data = dataset_world,
      aes(
        {{ variable_inputs }}/`Population (Mhab)`, country,
        fill=input_flow
        ),
      alpha=.85
      ) +
    
    #error bars of excretions, based on min and max values of excretions
    # geom_linerange(
    #   data = dataset_world,
    #   aes(
    #     xmin={{ variable_excretions }}/`Population (Mhab)`*(1-percent_food_waste_min)/(1-percent_food_waste),
    #     xmax={{ variable_excretions }}/`Population (Mhab)`*(1-percent_food_waste_max)/(1-percent_food_waste),
    #     y=country,
    #     color="human excreta"
    #     )
    #   ) +

    #points excretions N in red
    geom_point(
      data = dataset_world,
      aes({{ variable_excretions }}/`Population (Mhab)`, country, color="human excreta")
      ) 
  return(gg)
}


f_graph_potential_percent <- function(dataset, variable, variable_order){
  
  gg <- ggplot(dataset) +
    
    #columns
    geom_col(
      aes(
        {{ variable }}, reorder(country, {{ variable_order }}/`Population (Mhab)`)
        ),
      alpha=.85
      ) +
    
    #add % value atht the end of the column
    geom_text(
      aes(
        {{ variable }}, country, 
        label = paste0(round({{ variable }}*100), "%")
        ),
      hjust=0
      ) +
    # geom_text(
    #   aes(
    #     {{ variable }}, country, 
    #     label = paste0(
    #       "(",
    #       round({{ variable }}*(1-percent_food_waste_max)/(1-percent_food_waste)*100), 
    #       "-",
    #       round({{ variable }}*(1-percent_food_waste_min)/(1-percent_food_waste)*100),
    #       ")"
    #       )
    #     ),
    #   hjust=-1.3, vjust=0.5, size=2, fontface="italic"
    #   ) +
    
    
    #other options
    theme(
      axis.text.y = element_text(hjust=0.5),
      axis.ticks.y = element_blank(),
      axis.line.y = element_blank(),
      panel.grid.major.y= element_blank(),
      panel.grid.major.x= element_line(),
      panel.grid.minor.x= element_line()
      ) +
    labs(
      x="", y=""
      ) +
    coord_cartesian(
      xlim = c(0, 1)
      )  +
    scale_x_continuous(
      labels=scales::percent,
      breaks=seq(0, 10, .2),
      minor_breaks = seq(0, 10, .05)
      ) 
    
  return(gg)
  
}

f_graph_potential_percent_add_world <- function(graph_countries, dataset_world, variable){
    
    #columns
    gg <- graph_countries +
      geom_col(
      data = dataset_world,
      aes(
        {{ variable }}, country
        ),
      alpha=.85
      ) +

    #add % value atht the end of the column
    geom_text(
      data = dataset_world,
      aes(
        {{ variable }}, country,
        label = paste0(round({{ variable }}*100), "%")
        ),
      hjust=0
      ) 
    # geom_text(
    #   data = dataset_world,
    #   aes(
    #     {{ variable }}, country,
    #     label = paste0(
    #       "(",
    #       round({{ variable }}*(1-percent_food_waste_max)/(1-percent_food_waste)*100),
    #       "-",
    #       round({{ variable }}*(1-percent_food_waste_min)/(1-percent_food_waste)*100),
    #       ")"
    #       )
    #     ),
    #   hjust=-1.3, vjust=0.5, size=2, fontface="italic"
    #   )
  return(gg)
}

% of excreted manure applied to soils

f_prepare_share_manure_to_soil <- function(dataset){
  variable_subcategories_ruminants <- c(
  "ruminants Amount excreted in manure (ktN)", 
  "ruminants Losses from manure treated (ktN)", 
  "ruminants Manure applied to soils (ktN)", 
  "ruminants Manure left on pasture (ktN)"
  )
  variable_subcategories_non_ruminants <- c(
    "non ruminants Amount excreted in manure (ktN)", 
    "non ruminants Losses from manure treated (ktN)", 
    "non ruminants Manure applied to soils (ktN)", 
    "non ruminants Manure left on pasture (ktN)"
  )
  variable_subcategories_all <- c(
    "All Animals Amount excreted in manure (ktN)", 
    "All Animals Losses from manure treated (ktN)", 
    "All Animals Manure applied to soils (ktN)", 
    "All Animals Manure left on pasture (ktN)"
  )
  
  #prepare ruminants
  temp1 <- dataset %>%
    select(country, all_of(variable_subcategories_ruminants)) %>%
    mutate(
      percent_applied = `ruminants Manure applied to soils (ktN)`/`ruminants Amount excreted in manure (ktN)`,
      species = "ruminants") %>%
    gather(subgroup, subgroup_value, -c(country, percent_applied, species)) 
  temp1$subgroup <- gsub(" \\(ktN\\)", "", temp1$subgroup)
  temp1$subgroup <- gsub("ruminants ", "", temp1$subgroup)
  
  #prepare non ruminants
  temp2 <- dataset %>%
    select(country, all_of(variable_subcategories_non_ruminants)) %>%
    mutate(
      percent_applied = `non ruminants Manure applied to soils (ktN)`/`non ruminants Amount excreted in manure (ktN)`,
      species = "non ruminants") %>%
    gather(subgroup, subgroup_value, -c(country, percent_applied, species)) 
  temp2$subgroup <- gsub(" \\(ktN\\)", "", temp2$subgroup)
  temp2$subgroup <- gsub("non ruminants ", "", temp2$subgroup)
  
  #prepare all
  temp3 <- dataset %>%
    select(country, all_of(variable_subcategories_all)) %>%
    mutate(
      percent_applied = `All Animals Manure applied to soils (ktN)`/`All Animals Amount excreted in manure (ktN)`,
      species = "all livestock") %>%
    gather(subgroup, subgroup_value, -c(country, percent_applied, species)) 
  temp3$subgroup <- gsub(" \\(ktN\\)", "", temp3$subgroup)
  temp3$subgroup <- gsub("All Animals ", "", temp3$subgroup)
  
  #binds them together
  temp <- bind_rows(temp1, temp2, temp3)
  
  
  temp <- temp %>%
  filter(subgroup!="Amount excreted in manure") %>%
  mutate(
    subgroup = case_when(
      subgroup == "Manure left on pasture" ~ "Left on pasture",
      subgroup == "Losses from manure treated" ~ "Collected but volatilized\nduring storage",
      subgroup == "Manure applied to soils" ~ "Collected and spread\non agricultural land"
    )
  )
  
  
  #order categories for legend
  variable_subcategories <- c("Left on pasture", "Collected but volatilized\nduring storage", "Collected and spread\non agricultural land")

  temp$subgroup <- factor(
    temp$subgroup,
    levels = c(variable_subcategories)
  )
  
  return(temp)
  
}
temp_share_manure_soils <- f_prepare_share_manure_to_soil(FAO_countries_2018_2022)
temp_share_manure_soils_world <- f_prepare_share_manure_to_soil(FAO_world_2018_2022)

f_graph_manure_applied_losses <- function(dataset, dataset_world, countries_to_select){
  
  temp <- dataset %>% filter(country %in% countries_to_select)
  
  gg <- ggplot() +
    #countries
    geom_col(
      data = temp %>% filter(subgroup != "Amount excreted in manure"),
      aes(subgroup_value, reorder(country, percent_applied), fill=subgroup),
      position = "fill", alpha=.85
      ) +
    geom_text(
      data = temp %>% select(country, percent_applied, species) %>% distinct(),
      aes(0, country, label=paste0(round(percent_applied*100), "%")),
      hjust=0, fontface="italic", size=3
      ) +
    
    #world
    geom_col(
      data = dataset_world %>% filter(subgroup != "Amount excreted in manure"),
      aes(subgroup_value, country, fill=subgroup),
      position = "fill"
      ) +
    geom_text(
      data = dataset_world %>% select(country, percent_applied, species) %>% distinct(),
      aes(0, country, label=paste0(round(percent_applied*100), "%")),
      hjust=0, fontface="italic", size=3
      ) +
    
    
    
    scale_x_continuous(
      labels = scales::percent,
      breaks = seq(0, 1, 0.2)
      ) +
    theme(legend.position = "top") +
    
    #labels
    labs(
        color="", x="", y="",
        fill="repartition of N in\nexcreted manure",
        title = "Share of excreted manure collected and applied to agricultural land",
        caption = Source
      ) +
      facet_wrap(vars(species))
    
  return(gg)
}

gg <- f_graph_manure_applied_losses(
  temp_share_manure_soils, temp_share_manure_soils_world,
  countries_to_select_large
  )
gg

f_save_graph_pdf_png(
  gg, "graph/potential/", "fate_manure",
  resolution_dpi, grahps_height, grahps_width
)

% of synthetic fertilizers covered by urine

2018-2022

#select countries > 50 Mhab
temp <- temp_croplands_inputs %>% filter(country %in% countries_to_select_large)

gg1 <- f_graph_comparison_inputs_per_cap(
  temp %>% filter(input_flow=="Mineral fertilizers"),
  `N mineral fertilization (ktN)`, 
  `N mineral fertilization (ktN)`, 
  `human urine (ktN)`
  ) +
  scale_fill_manual(values="#fde725") +
  labs(subtitle = "N in urine and fertilizers consumption")
gg1 <- f_graph_comparison_inputs_per_cap_add_world(
  gg1,
  temp_croplands_inputs_world %>% filter(input_flow=="Mineral fertilizers"),
  `N mineral fertilization (ktN)`, 
  `human urine (ktN)`
)


gg2 <- f_graph_potential_percent(
  temp %>% select(-c(value, input_flow)) %>% distinct(), 
  `human urine % of synthetic fertilization`, 
  `N mineral fertilization (ktN)` 
  )  + 
  labs(subtitle ="% of fertilizers replaced by urine")
gg2 <- f_graph_potential_percent_add_world(
  gg2,
  temp_croplands_inputs_world %>% select(-c(value, input_flow)) %>% distinct(), 
  `human urine % of synthetic fertilization`
  ) +
  coord_cartesian(
    xlim = c(0, 2)
    ) +
  scale_x_continuous(
    labels=scales::percent,
    breaks=seq(0, 10, .2),
    minor_breaks = seq(0, 10, .05)
    )   


gg <- gg1 + theme(axis.text.y = element_blank()) + gg2 + 
  plot_annotation(
    title = "Comparison of N in urine to synthetic fertilizers consumption",
    caption = Source
    )
gg

f_save_graph_pdf_png(
  gg, "graph/potential/", "share_urine_vs_synth_ferti",
  resolution_dpi, grahps_height, grahps_width
)

% of cropland inputs

2018-2022

#select countries > 50 Mhab 
temp <- temp_croplands_inputs %>% filter(country %in% countries_to_select_large)

#Graphs
gg1 <- f_graph_comparison_inputs_per_cap(
  temp, value, `N croplands Input (ktN)`, `human excreta (ktN)` 
  ) + labs(fill="inputs to croplands", subtitle = "N in excretions vs cropland inputs")
gg1 <- f_graph_comparison_inputs_per_cap_add_world(
  gg1,
  temp_croplands_inputs_world, value, `human excreta (ktN)` 
)

gg2 <- f_graph_potential_percent(
  temp %>% select(-c(value, input_flow)) %>% distinct(), 
  `human excreta % of croplands inputs`, `N croplands Input (ktN)`
  ) + 
  labs(subtitle = "excretions as % of cropland inputs") +
  coord_cartesian(
      xlim = c(0, 1.35)
      )  
gg2 <- f_graph_potential_percent_add_world(
  gg2,
   temp_croplands_inputs_world %>% select(-c(value, input_flow)) %>% distinct(), 
  `human excreta % of croplands inputs`
)


gg <- gg1 + theme(axis.text.y = element_blank()) + gg2 + 
  plot_annotation(
    title = "human excreta compared to croplands N inputs",
    caption = Source
    )
gg

f_save_graph_pdf_png(
  gg, "graph/potential/", "share_excretions_vs_cropland_ferti",
  resolution_dpi, grahps_height, grahps_width
)

Comparison cropland inputs / total N

f_graph_ratio <- function(dataset, dataset_world, var_numerator, var_denominator){
  
  dataset <- dataset %>%
    mutate(
      ratio = {{ var_numerator }}/{{ var_denominator }}
    )
  
  dataset_world <- dataset_world %>%
    mutate(
      ratio = {{ var_numerator }}/{{ var_denominator }}
    )
  
  ggplot(dataset) +
    
    #countries
    geom_col(
      aes(ratio, reorder(country, ratio)),
      fill="#440154"
      ) +
    
    #world
    geom_col(
      data = dataset_world,
      aes(ratio, country),
      fill="#440154"
      ) +
    
    #other options
    scale_x_continuous(
      labels=scales::percent,
      breaks = seq(0, 2, .2),
      minor_breaks = seq(0, 2, .1)
    ) +
    geom_vline(xintercept = 1, linetype="dashed") +
    theme(
      panel.grid.major.y = element_blank(),
      panel.grid.major.x = element_line(),
      panel.grid.minor.x = element_line(),
    ) +
    labs(
      y="", x=""
    )
  
}

#select countries > 50 Mhab 
temp <- f_select_countries(FAO_countries_2018_2022, countries_to_select_large)

gg1 <- f_graph_ratio(temp, FAO_world_2018_2022, `N croplands Manure applied to Soils (ktN)`, `All Animals Manure applied to soils (ktN)`) +
  labs(
    title = "...manure"
    )
gg2 <- f_graph_ratio(temp, FAO_world_2018_2022, `N croplands Mineral fertilizers (ktN)`, `N mineral fertilization (ktN)`) +
  labs(
    title  = "...synthetic fertilizers",
    )
gg <- gg1 + gg2 +
  plot_annotation(
      title = "Cropland N fertilization as % of Agricultural Land N fertilization, for...",
      caption = Source
      )
gg

f_save_graph_pdf_png(
  gg, "graph/potential/", "total_vs_cropland_ferti",
  resolution_dpi, grahps_height, grahps_width
)

% manure handled and applied to soils (croplands and grasslands)

recovery_human_excretions <- 0.8

f_prepare_manure_to_soils <- function(dataset){
  #prepare inputs file
  dataset <- dataset %>%
    select(
      #country and population
      country, `Population (Mhab)`, 
      
      #total excretions,
      `human excreta (ktN)`, `human excreta % of manure applied to soils`, percent_food_waste,
      
      #inputs to croplands
      `All Animals Manure applied to soils (ktN)`,
      `non ruminants Manure applied to soils (ktN)`,
      `ruminants Manure applied to soils (ktN)` 
  
      
      ) %>%
    
    gather(
      input_flow, value, 
      c(`non ruminants Manure applied to soils (ktN)`, 
        `ruminants Manure applied to soils (ktN)`)
      ) %>%
    
    mutate(
      input_flow = gsub(" Manure applied to soils \\(ktN\\)", "", input_flow)
      ) 
  
  #order ruminants / non ruminants legend
  dataset$input_flow <- factor(
    dataset$input_flow,
    levels = c("ruminants", "non ruminants")
  )
  
  #apply recovery rate to human excreta
  dataset <- dataset %>%
    mutate(
      `human excreta (ktN)` = `human excreta (ktN)` * recovery_human_excretions,
      `human excreta % of manure applied to soils` = `human excreta % of manure applied to soils` * recovery_human_excretions
      )
  
  return(dataset)
  
}

temp_manure_to_soils <- f_prepare_manure_to_soils(FAO_countries_2018_2022)
temp_manure_to_soils_world <- f_prepare_manure_to_soils(FAO_world_2018_2022)

2018-2022

#select countries > 50 Mhab
temp <- f_select_countries(temp_manure_to_soils, countries_to_select_large)

#Graphs
gg1 <- f_graph_comparison_inputs_per_cap(
  temp, value, `All Animals Manure applied to soils (ktN)`, `human excreta (ktN)` 
  ) +
  #coord_cartesian(xlim=c(15, 0)) +
  labs(
    fill="manure applied to soils", subtitle = "human excreta compared to manure"
  )
gg1 <- f_graph_comparison_inputs_per_cap_add_world(
  gg1,
  temp_manure_to_soils_world, value, `human excreta (ktN)` 
)


gg2 <- f_graph_potential_percent(
  temp %>% select(-c(value, input_flow)) %>% distinct(), 
  `human excreta % of manure applied to soils`, `All Animals Manure applied to soils (ktN)`
  ) +
  coord_cartesian(xlim = c(0, 2))  +
  labs(
    subtitle = "human excreta as % of manure"
  )
gg2 <- f_graph_potential_percent_add_world(
  gg2,
  temp_manure_to_soils_world %>% select(-c(value, input_flow)) %>% distinct(), 
  `human excreta % of manure applied to soils`
)



gg <- gg1 + theme(axis.text.y = element_blank()) + gg2 + 
  plot_annotation(
    title = "human excreta (80% recovery) compared to manure applied to soils",
    caption = Source
    )
gg

f_save_graph_pdf_png(
  gg, "graph/potential/", "share_excretions_vs_manure_applied",
  resolution_dpi, grahps_height, grahps_width
)

% manure applied to croplands

f_prepare_manure_to_croplands <- function(dataset){
  #prepare inputs file
  dataset <- dataset %>%
    select(
      #country and population
      country, `Population (Mhab)`, 
      
      #total excretions,
      `human excreta (ktN)`, percent_food_waste,
      
      #manure to croplands
      `applied to croplands`  = `N croplands Manure applied to Soils (ktN)` ,
      `Manure applied to croplands`  = `N croplands Manure applied to Soils (ktN)` ,
      
      # all excretions
      `All Animals Amount excreted in manure (ktN)` 
  
      ) %>%
    
    mutate(
      `human excreta % of manure applied to croplands` = `human excreta (ktN)`/`Manure applied to croplands` ,
      `not applied to croplands` = `All Animals Amount excreted in manure (ktN)` -`Manure applied to croplands`
    ) %>%
    
    gather(
      input_flow, value, `applied to croplands` , `not applied to croplands`
      ) 
  
  #apply recovery rate to human excreta
  dataset <- dataset %>%
    mutate(
      `human excreta (ktN)` = `human excreta (ktN)` * recovery_human_excretions,
      `human excreta % of manure applied to croplands` = `human excreta % of manure applied to croplands` * recovery_human_excretions
      )
  
  return(dataset)
  
}

2018-2022

80% recovery

recovery_human_excretions <- 0.8
temp_manure_to_croplands <- f_prepare_manure_to_croplands(FAO_countries_2018_2022)
temp_manure_to_croplands_world <- f_prepare_manure_to_croplands(FAO_world_2018_2022)

temp_manure_to_croplands$input_flow <- factor(
  temp_manure_to_croplands$input_flow,
  levels = c(
    "not applied to croplands",
    "applied to croplands"
  )
)
temp_manure_to_croplands_world$input_flow <- factor(
  temp_manure_to_croplands_world$input_flow,
  levels = c(
    "not applied to croplands",
    "applied to croplands"
  )
)

temp <- temp_manure_to_croplands %>% filter(country %in% countries_to_select_large)

#Graphs
gg1 <- f_graph_comparison_inputs_per_cap(
  temp, value, `Manure applied to croplands`, `human excreta (ktN)` 
  ) + labs(fill="manure excreted and...", subtitle = "human excreta compared to manure")
gg1 <- f_graph_comparison_inputs_per_cap_add_world(
  gg1,
  temp_manure_to_croplands_world, value, `human excreta (ktN)` 
) +
  scale_fill_manual(
    values = c("applied to croplands" = "#35b779",   
               "not applied to croplands" = "#c4eed9")    
    ) +
  theme(legend.position = c(.3, .51))
gg1

gg2 <- f_graph_potential_percent(
  temp %>% select(-c(value, input_flow)) %>% distinct(), 
  `human excreta % of manure applied to croplands`, `Manure applied to croplands`
  ) + 
  labs(subtitle = "excretions as % of manure applied to cropland") +
  coord_cartesian(
      xlim = c(0, 2.5)
      ) 
gg2 <- f_graph_potential_percent_add_world(
  gg2,
   temp_manure_to_croplands_world %>% select(-c(value, input_flow)) %>% distinct(), 
  `human excreta % of manure applied to croplands`
) +
  scale_x_continuous(
      labels=scales::percent,
      breaks=seq(0, 10, .5),
      minor_breaks = seq(0, 10, .1)
      ) 

gg <- gg1 + theme(axis.text.y = element_blank()) + gg2 + 
  plot_annotation(
    title = "human excreta (80% recovery) compared to manure applied to croplands",
    caption = Source
    )
gg

f_save_graph_pdf_png(
  gg, "graph/potential/", "share_excretions_vs_manure_on_croplands_80perc",
  resolution_dpi, grahps_height, grahps_width
)

100% recovery

recovery_human_excretions <- 1
temp_manure_to_croplands <- f_prepare_manure_to_croplands(FAO_countries_2018_2022)
temp_manure_to_croplands_world <- f_prepare_manure_to_croplands(FAO_world_2018_2022)

temp_manure_to_croplands$input_flow <- factor(
  temp_manure_to_croplands$input_flow,
  levels = c(
    "not applied to croplands",
    "applied to croplands"
  )
)
temp_manure_to_croplands_world$input_flow <- factor(
  temp_manure_to_croplands_world$input_flow,
  levels = c(
    "not applied to croplands",
    "applied to croplands"
  )
)

temp <- temp_manure_to_croplands %>% filter(country %in% countries_to_select_large)

#Graphs
gg1 <- f_graph_comparison_inputs_per_cap(
  temp, value, `Manure applied to croplands`, `human excreta (ktN)` 
  ) + labs(fill="manure excreted and...", subtitle = "human excreta compared to manure")
gg1 <- f_graph_comparison_inputs_per_cap_add_world(
  gg1,
  temp_manure_to_croplands_world, value, `human excreta (ktN)` 
) +
  scale_fill_manual(
    values = c("applied to croplands" = "#35b779",   
               "not applied to croplands" = "#c4eed9")    
    ) +
  theme(legend.position = c(.3, .51))
gg1

gg2 <- f_graph_potential_percent(
  temp %>% select(-c(value, input_flow)) %>% distinct(), 
  `human excreta % of manure applied to croplands`, `Manure applied to croplands`
  ) + 
  labs(subtitle = "excretions as % of manure applied to cropland") +
  coord_cartesian(
      xlim = c(0, 3.2)
      ) 
gg2 <- f_graph_potential_percent_add_world(
  gg2,
   temp_manure_to_croplands_world %>% select(-c(value, input_flow)) %>% distinct(), 
  `human excreta % of manure applied to croplands`
) +
  scale_x_continuous(
      labels=scales::percent,
      breaks=seq(0, 10, .5),
      minor_breaks = seq(0, 10, .1)
      ) 

gg <- gg1 + theme(axis.text.y = element_blank()) + gg2 + 
  plot_annotation(
    title = "human excreta compared to manure excreted and applied to croplands",
    caption = Source
    )
gg

f_save_graph_pdf_png(
  gg, "graph/potential/", "share_excretions_vs_manure_on_croplands_100perc",
  resolution_dpi, grahps_height, grahps_width
)

% manure excreted

f_prepare_manue_excreted <- function(dataset){
  #prepare inputs file
  dataset <- dataset %>%
    select(
      #country and population
      country, `Population (Mhab)`, 
      
      #total excretions,
      `human excreta (ktN)`, `human excreta % of livestock excretion`, percent_food_waste,
      
      #inputs to croplands
      `All Animals Amount excreted in manure (ktN)`,
      `non ruminants Amount excreted in manure (ktN)`,
      `ruminants Amount excreted in manure (ktN)`
  
      
      ) %>%
    
    gather(
      input_flow, value, 
      c(`non ruminants Amount excreted in manure (ktN)`,
        `ruminants Amount excreted in manure (ktN)`)
      ) %>%
    
    mutate(
      input_flow = gsub(" Amount excreted in manure \\(ktN\\)", "", input_flow)
      ) 
  
  #order ruminants / non ruminants legend
  dataset$input_flow <- factor(
    dataset$input_flow,
    levels = c("ruminants", "non ruminants")
  )
  
  return(dataset)
  
}
temp_manure_excreted <- f_prepare_manue_excreted(FAO_countries_2018_2022)
temp_manure_excreted_world <- f_prepare_manue_excreted(FAO_world_2018_2022)

2018-2022

#select countries > 50 Mhab
temp <- f_select_countries(temp_manure_excreted, countries_to_select_large)

#Graphs
gg1 <- 
  f_graph_comparison_inputs_per_cap(
    temp, value, 
    `All Animals Amount excreted in manure (ktN)`, 
    `human excreta (ktN)` 
    ) +
  coord_cartesian(xlim = c(55, 0)) +
  labs(
    fill="manure excreted", subtitle = "human excreta compared to manure"
    )
gg1 <- f_graph_comparison_inputs_per_cap_add_world(
  gg1,
  temp_manure_excreted_world, value,`human excreta (ktN)`
  )

gg2 <- 
  f_graph_potential_percent(
    temp %>% select(-c(value, input_flow)) %>% distinct(), 
    `human excreta % of livestock excretion`, 
    `All Animals Amount excreted in manure (ktN)`
    ) +
  coord_cartesian(xlim = c(0, 1.1))  +
  labs(
    subtitle = "human excreta as % of manure"
    )
gg2 <- f_graph_potential_percent_add_world(
  gg2, temp_manure_excreted_world %>% select(-c(value, input_flow)) %>% distinct(), 
  `human excreta % of livestock excretion`
)



gg <- gg1 + theme(axis.text.y = element_blank()) + gg2 + 
  plot_annotation(
    title = "human excreta compared to manure excretions",
    caption = Source
    )
gg

f_save_graph_pdf_png(
  gg, "graph/potential/", "share_excretions_vs_manure",
  resolution_dpi, grahps_height, grahps_width
)

Share animal / vegetal food supply

https://ourworldindata.org/grapher/daily-protein-supply-from-animal-and-plant-based-foods

supply_veg_anim <- read_csv("source/daily-protein-supply-from-animal-and-plant-based-foods.csv")
supply_veg_anim_world <- supply_veg_anim %>% filter(Entity=="World")

f_prepare_food_supply <- function(dataset){
  
  dataset <- dataset %>%
    select(
      country = Entity, Year, 
      animal = `Animal Products | 00002941 || Food available for consumption | 0674pc || grams of protein per day per capita`,
      vegetal = `Vegetal Products | 00002903 || Food available for consumption | 0674pc || grams of protein per day per capita`
      ) %>%
    mutate(
      animal = animal*365/1000/6.25,
      vegetal = vegetal*365/1000/6.25,
      total_kgN_per_cap = animal + vegetal,
      share_vegetal = vegetal/total_kgN_per_cap
      )
  
  dataset <- dataset %>%
    filter(Year>=2018 & Year <=2022) %>%
    select(-Year) %>%
    group_by(country) %>%
    summarise_all(mean)
  
  dataset <- dataset %>%
    gather(
      animal_or_vegetal, animal_or_vegetal_value, vegetal, animal
    )
  
  return(dataset)
  
}
supply_veg_anim <- f_prepare_food_supply(supply_veg_anim)
supply_veg_anim_world <- f_prepare_food_supply(supply_veg_anim_world)



supply_veg_anim <- supply_veg_anim %>%
  mutate(
    country = case_when(
      country=="United Kingdom" ~ "UK",
      country=="United States" ~ "USA",
      country=="South Korea" ~ "Korea",
      country=="Vietnam" ~ "Viet Nam",
      country=="Democratic Republic of Congo" ~ "Dem. Rep. Congo",
      T ~ country
    )
  )

f_graph_N_supply <- function(dataset, dataset_world, var_animal_or_vegetal_value, var_share_vegetal){
  
  gg1 <- ggplot(dataset) +
    #countries
    geom_col(
      aes(animal_or_vegetal_value, reorder(country, {{ var_share_vegetal }}), fill={{ var_animal_or_vegetal_value }}),
      position = "fill", color="black"
      ) +
    geom_text(
      data = dataset %>% select(country, {{ var_share_vegetal }}) %>% distinct(),
      aes(0, country, label = paste0(round({{ var_share_vegetal }}*100), "%")),
      hjust=1, fontface="italic"
    ) +
    
    #world
    geom_col(
      data = dataset_world,
      aes(animal_or_vegetal_value, country, fill={{ var_animal_or_vegetal_value }}),
      position = "fill", color="black"
      ) +
    geom_text(
      data = dataset_world %>% select(country, {{ var_share_vegetal }}) %>% distinct(),
      aes(0, country, label = paste0(round({{ var_share_vegetal }}*100), "%")),
      hjust=1, fontface="italic"
    ) +
    
    
    scale_x_reverse(
      labels = scales::percent
    ) +
    theme(
      axis.ticks.y = element_blank(),
      axis.text.y = element_blank(),
      axis.line.y = element_blank(),
      legend.position = "top"
    ) +
    labs(
      y="", x="", fill="", subtitle = "relative shares"
      )
  
  gg2 <- ggplot(dataset) +
    #countries
    geom_col(
      aes(animal_or_vegetal_value, reorder(country, {{ var_share_vegetal }}), fill={{ var_animal_or_vegetal_value }})
      )+
    geom_col(
      data = dataset %>% select(country, total_kgN_per_cap, {{ var_share_vegetal }}) %>% distinct(),
      aes(total_kgN_per_cap, reorder(country, {{ var_share_vegetal }})),
      fill="transparent", color="black"
      ) +
    geom_text(
      data = dataset %>% select(country, total_kgN_per_cap) %>% distinct(),
      aes(total_kgN_per_cap, country, label = round(total_kgN_per_cap, 1)),
      hjust=-.3, fontface="italic"
    ) +
    
    #world
    geom_col(
      data = dataset_world,
      aes(animal_or_vegetal_value, country, fill={{ var_animal_or_vegetal_value }})
      )+
    geom_col(
      data = dataset_world %>% select(country, total_kgN_per_cap, {{ var_share_vegetal }}) %>% distinct(),
      aes(total_kgN_per_cap, country),
      fill="transparent", color="black"
      ) +
    geom_text(
      data = dataset_world %>% select(country, total_kgN_per_cap) %>% distinct(),
      aes(total_kgN_per_cap, country, label = round(total_kgN_per_cap, 1)),
      hjust=-.3, fontface="italic"
    ) +
    
    
    
    theme(
      legend.position = "none",
      axis.text.y = element_text(hjust=0.5)
      ) +
    labs(y="", x=expression("kgN" %.% capita^-1 %.% year^-1), subtitle = "absolute quantities") +
    xlim(0, 8.5)
  
  gg <- gg1 + gg2 + 
    plot_annotation(
      title = "Repartition of vegetal and animal products in N food supply",
      caption = Source
      )
  
  return(gg)
}

temp <- supply_veg_anim %>% filter(country %in% countries_to_select_large)

gg <- f_graph_N_supply(temp, supply_veg_anim_world, animal_or_vegetal, share_vegetal)
gg

f_save_graph_pdf_png(
  gg, "graph/potential/", "food_supply_animal_vegetal",
  resolution_dpi, grahps_height, grahps_width
)

Share ruminants vs non ruminants

temp <- f_select_countries(FAO_countries_2018_2022, countries_to_select_large)


f_prepare_share_ruminants <- function(dataset) {
  temp <- dataset %>%
    select(
      country,
      `All Animals Amount excreted in manure (ktN)`, 
      `non ruminants Amount excreted in manure (ktN)`, 
      `ruminants Amount excreted in manure (ktN)`
      ) %>%
    mutate(
      `non ruminants` = `non ruminants Amount excreted in manure (ktN)`/`All Animals Amount excreted in manure (ktN)`,
      ruminants = `ruminants Amount excreted in manure (ktN)`/`All Animals Amount excreted in manure (ktN)`
    ) %>%
    select(country, `non ruminants`, ruminants)
  
  temp$country <- factor(
    temp$country,
    levels = arrange(temp, `non ruminants`) %>% pull(country)
  )
  
  temp <- temp %>%
    gather(
      livestock_type, value, `non ruminants`, ruminants
    )
  
  temp$livestock_type <- factor(
    temp$livestock_type,
    levels = c("ruminants", "non ruminants")
  )
  
  return(temp)
  
}
temp <- f_prepare_share_ruminants(temp)
temp_world <- f_prepare_share_ruminants(FAO_world_2018_2022)




gg <- ggplot(temp) +
  
  #countries bars and labels
  geom_col(aes(value, country, fill=livestock_type)) +
  geom_text(
    data = temp %>% filter(livestock_type=="non ruminants"),
    aes(0, country, label = paste(round(value*100), "%")),
    hjust=0
  ) +
  
  #world
  geom_col(
    data = temp_world,
    aes(value, country, fill=livestock_type)
    ) +
  geom_text(
    data = temp_world %>% filter(livestock_type=="non ruminants"),
    aes(0, country, label = paste(round(value*100), "%")),
    hjust=0
  ) +
  
  #other options
  scale_x_continuous(
    labels = scales::percent
  ) +
  labs(
    x="", y="", fill="",
    title = "% of non ruminants excretions in total livestock excretions",
    caption = Source
    ) 
gg

f_save_graph_pdf_png(
  gg, "graph/potential/", "share_non_ruminants_excretions",
  resolution_dpi, grahps_height, 7
)

rm(list = ls())