---
title: "Sustainable commodity sourcing requires measuring and governing land use change at multiple scales"
author: "Erasmus zu Ermgassen"
date: '2023-08-28'
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, message = F, warning = F, fig.align = 'center', fig.height = 8)

# Workstation set-up ------------------------------------------------------

library(tidyverse)
library(scales)
library(sf)
library(extrafont)
library(patchwork)


# Load the data file
# >>> Update the directory to your local directory
load(here::here("SCALE.rdata"))


# The above loads the following objects:
# cov_muni - Mapbiomas v7.1 data on the land cover per municipality
# trans_muni - Mapbiomas v7.1 data on land cover transitions per municipality
# property - the land cover changes in properties in Brazil from 1985-2021, calculated by intersecting Mapbiomas v7.1 with data from Sparovek et al. 2019, Land Use Policy
# property_2000 - the land cover in properties in Brazil from 2000-2021, calculated by intersecting Mapbiomas v7.1 with data from Sparovek et al. 2019, Land Use Policy
# geocode_dic - dictionary of IBGE municipal codes and names
# state_dic - dictionary of IBGE state codes and names
# non_veg_lu - vector of non-vegetative land uses
# natural_lu - vector of natural vegetation types
# ag_lu - vector of agricultural land use types
# biome - municipailities classified according to their dominant biome
# mt_sf - the IBGE municipal code of each property polygon in Mato Grosso
# str_trans - function to remove special characters from names
# state_abbr - list of state abbreviations
# cbPalette - the colour scale
```

``` {r data-processing}
# Make a dictionary of class_id -> land use class name
class_dic <- cov_muni %>% distinct(class_id, level_4)
```

## Figure 1a 

Plot changes in major land uses over time.
``` {r}
p_line <- cov_muni %>%
  mutate(region = str_split(city, pattern = " - "),
         city = map_chr(region,1),
         state = map_chr(region,2), 
         city = str_trans(city)) %>%
  select(-class_id,-city, -region) %>%
  gather(year, area, `1985`:`2021`) %>%
  filter(level_4 %in% ag_lu) %>%
  group_by(level_4, year) %>%
  summarise(sum_area = sum(area, na.rm = T)) %>%
  ungroup() %>%
  mutate(level_4 = as_factor(level_4), 
         level_4 = fct_relevel(level_4, 
                               levels = ag_lu)) %>%
  ggplot(aes(as.numeric(year), sum_area, col = level_4)) + 
  geom_line(size = 2) +
  theme_classic() +
  theme(legend.title = element_blank(), 
        axis.title.x = element_blank(),
        text=element_text(size=16,  family= "DM Sans"), 
        legend.position = "bottom") +
  scale_color_manual(values = cbPalette) +
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6, sep = "")) +
  labs(y = "Area (hectares)") +
  guides(color = guide_legend(nrow = 3))
```

## Figure 1B

#### 1985-2021 PAA CALCULATION

Here we calculate PAA at five different scales: pixel, property, municipality, state, and country.

``` {r}
# PIXEL
paa_pixel_1985 <- 
  trans_muni %>%
  mutate(region = str_split(city, pattern = " - "),
         state = map_chr(region,2)) %>%
  filter(from_level_4 %in% natural_lu, 
         !to_level_4 %in% natural_lu, 
         from_level_4 != to_level_4) %>%
  filter(!to_level_4 %in% non_veg_lu) %>%
  group_by(to_level_4) %>%
  summarise(PAA = sum(`1985-2021`, na.rm = T)) %>%
  ungroup() %>%
  rename(level_4 = to_level_4)


# Property level
# 'gid' is the polygon ID
# NB we run the analysis per state to speed up the operation
# ... takes ca 5 mins on a laptop
property_analysis <- function(state_abbr, df){
  df %>%
    filter(state == state_abbr) %>%
    rename(class_id = class) %>%
    left_join(class_dic, by = "class_id") %>%
    select(-`.geo`,-`system:index`) %>%
    filter(!level_4 %in% non_veg_lu) %>%
    mutate(level_4 = ifelse(level_4 %in% natural_lu, "Native Vegetation", level_4)) %>%
    spread(year, area, fill = 0) %>%
    group_by(gid, level_4) %>%
    summarise(`2021` = sum(`2021`, na.rm = T), 
              `1985` = sum(`1985`, na.rm = T)) %>%
    ungroup() %>%
    mutate(
      change = `2021` - `1985`,
      direction = ifelse(change < 0, "contraction", "expansion")
    )  %>%
      filter(change != 0)
}
change_property_1985 <- map_df(
  .x = state_abbr,
  .f = property_analysis,
  df = property
  ) 
paf_property_1985 <- change_property_1985 %>%
  filter(direction == "expansion") %>%
  filter(level_4 %in% ag_lu) %>%
  group_by(gid) %>%
  mutate(total_change = sum(change, na.rm = T)) %>%
  group_by(gid, level_4) %>%
  summarise(PAF = sum(change, na.rm = T) / unique(total_change)) %>%
  ungroup()
paa_property_1985 <- change_property_1985 %>%
  filter(direction == "contraction") %>%
  filter(level_4  == "Native Vegetation") %>%
  group_by(gid) %>%
  summarise(native_veg_change = sum(change, na.rm = T) * -1) %>%
  ungroup() %>%
  left_join(paf_property_1985, by = "gid") %>%
  mutate(PAA = native_veg_change * PAF) %>%
  drop_na(PAA)


# MUNI
paf_muni_1985 <- cov_muni %>%
  filter(level_4 %in% ag_lu) %>%
  mutate(change = `2021`- `1985`,
         direction = ifelse(change <= 0, "contraction", "expansion")) %>%
  filter(direction == "expansion") %>%
  group_by(city) %>%
  mutate(total_change = sum(change, na.rm = T)) %>%
  group_by(city, level_4) %>%
  summarise(PAF = sum(change, na.rm = T) / unique(total_change)) %>%
  ungroup()
paa_muni_1985 <- 
  cov_muni %>%
  filter(level_4 %in% natural_lu) %>%
  mutate(change = `1985` - `2021`,
         direction = ifelse(change > 0, "contraction", "expansion")) %>%
  filter(direction == "contraction") %>%
  group_by(city) %>%
  summarise(native_veg_change = sum(change, na.rm = T)) %>%
  ungroup() %>%
  left_join(paf_muni_1985, by = "city") %>%
  mutate(PAA = native_veg_change * PAF) 


# STATE
paf_state_1985 <- cov_muni %>%
  mutate(region = str_split(city, pattern = " - "),
         state = map_chr(region,2)) %>%
  filter(level_4 %in% ag_lu) %>%
  group_by(state, level_4) %>%
  summarise(`2021` = sum(`2021`, na.rm = T),
            `1985` = sum(`1985`, na.rm = T)) %>%
  ungroup() %>%
  mutate(change = `2021`- `1985`,
         direction = ifelse(change <= 0, "contraction", "expansion")) %>%
  filter(direction == "expansion") %>%
  group_by(state) %>%
  mutate(total_change = sum(change, na.rm = T)) %>%
  group_by(state, level_4) %>%
  summarise(PAF = sum(change, na.rm = T) / unique(total_change)) %>%
  ungroup()
paa_state_1985 <- 
  cov_muni %>%
  mutate(region = str_split(city, pattern = " - "),
         state = map_chr(region,2)) %>%
  filter(level_4 %in% natural_lu) %>%
  mutate(change = `1985` - `2021`,
         direction = ifelse(change > 0, "contraction", "expansion")) %>%
  filter(direction == "contraction") %>%
  group_by(state) %>%
  summarise(native_veg_change = sum(change, na.rm = T)) %>%
  ungroup() %>%
  left_join(paf_state_1985, by = "state") %>%
  mutate(PAA = native_veg_change * PAF) 


# NATIONAL
paf_national_1985 <- cov_muni %>%
  filter(level_4 %in% ag_lu) %>%
  mutate(level_4 = ifelse(level_4 %in% natural_lu, "Native Vegetation", level_4)) %>%
  group_by(level_4) %>%
  summarise(`2021` = sum(`2021`, na.rm = T),
            `1985` = sum(`1985`, na.rm = T)) %>%
  ungroup() %>%
  mutate(change = `2021`- `1985`,
         direction = ifelse(change <= 0, "contraction", "expansion")) %>%
  filter(direction == "expansion") %>%
  mutate(total_change = sum(change, na.rm = T)) %>%
  group_by(level_4) %>%
  summarise(PAF = sum(change, na.rm = T) / unique(total_change)) %>%
  ungroup() %>%
  mutate(resolution = "national")
paa_national_1985 <- 
  cov_muni %>%
  filter(level_4 %in% natural_lu) %>%
  mutate(change = `1985` - `2021`,
         direction = ifelse(change > 0, "contraction", "expansion")) %>%
  filter(direction == "contraction") %>%
  summarise(native_veg_change = sum(change, na.rm = T)) %>%
  mutate(resolution = "national") %>%
  left_join(paf_national_1985, by = "resolution") %>%
  mutate(PAA = native_veg_change * PAF)  %>%
  select(-resolution)


# Here simplifying the original land use 
pan_1985 <- 
  bind_rows(
    "National" = paa_national_1985,
    "State" = paa_state_1985,
    "Municipality" = paa_muni_1985,
    "Property" = paa_property_1985,
    "Pixel" = paa_pixel_1985,
    .id = "resolution"
  ) %>%
    mutate(resolution = as_factor(resolution),
           resolution = fct_relevel(resolution, levels = c("National", "State", "Municipality", "Property","Pixel"))) %>%
    group_by(resolution) %>%
    mutate(total_paa = sum(PAA),
           level_4 = ifelse(level_4 %in% c("Aquaculture", "Citrus", "Mining"), "Other", level_4)) %>%
    group_by(resolution) %>%
    mutate(total_paa = sum(PAA, na.rm = T)) %>%
    group_by(level_4, resolution) %>%
    summarise(sum_paa = sum(PAA),
              PAN = sum_paa / unique(total_paa) * 100) %>%
    ungroup()


# Visualise %
p_pan_1985 <- pan_1985 %>%
  ggplot(aes(resolution, PAN, fill = level_4)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = cbPalette) +
  theme_minimal() +
  theme(legend.title = element_blank(),
  text=element_text(size=16,  family= "DM Sans"), 
  legend.position = "none") +
  # theme_iluc() +
  labs(x = "Scale of analysis", 
       y = "% deforestation attributed to each driver\nin Brazil (1985-2021)") +
  coord_flip()
p_pan_1985


# tidy up 
rm(property)
```

## Figure 1C

#### 2000-2021 PAA CALCULATION

``` {r}
# PIXEL
paa_pixel_2000 <- 
  trans_muni %>%
  mutate(region = str_split(city, pattern = " - "),
         state = map_chr(region,2)) %>%
  filter(from_level_4 %in% natural_lu, 
         !to_level_4 %in% natural_lu, 
         from_level_4 != to_level_4) %>%
  filter(!to_level_4 %in% non_veg_lu) %>%
  group_by(to_level_4) %>%
  summarise(PAA = sum(`2000-2021`, na.rm = T)) %>%
  ungroup() %>%
  rename(level_4 = to_level_4)


# Property level
# NB we run the analysis per state to speed up the operation
# ... takes ca 5 mins on a laptop
property_analysis_2000 <- function(state_abbr, df){
  df %>%
    filter(state == state_abbr) %>%
    rename(class_id = class) %>%
    left_join(class_dic, by = "class_id") %>%
    select(-`.geo`,-`system:index`) %>%
    filter(!level_4 %in% non_veg_lu) %>%
    mutate(level_4 = ifelse(level_4 %in% natural_lu, "Native Vegetation", level_4)) %>%
    spread(year, area, fill = 0) %>%
    group_by(gid, level_4) %>%
    summarise(`2021` = sum(`2021`, na.rm = T), 
              `2000` = sum(`2000`, na.rm = T)) %>%
    ungroup() %>%
    mutate(
      change = `2021` - `2000`,
      direction = ifelse(change < 0, "contraction", "expansion")
    )  %>%
      filter(change != 0)
}
change_property_2000 <- map_df(
  .x = state_abbr,
  .f = property_analysis_2000,
  df = property_2000
  ) 
paf_property_2000 <- change_property_2000 %>%
  filter(direction == "expansion") %>%
  filter(level_4 %in% ag_lu) %>%
  group_by(gid) %>%
  mutate(total_change = sum(change, na.rm = T)) %>%
  group_by(gid, level_4) %>%
  summarise(PAF = sum(change, na.rm = T) / unique(total_change)) %>%
  ungroup()
paa_property_2000 <- change_property_2000 %>%
  filter(direction == "contraction") %>%
  filter(level_4  == "Native Vegetation") %>%
  group_by(gid) %>%
  summarise(native_veg_change = sum(change, na.rm = T) * -1) %>%
  ungroup() %>%
  left_join(paf_property_2000, by = "gid") %>%
  mutate(PAA = native_veg_change * PAF) %>%
  drop_na(PAA)


# MUNI
paf_muni_2000 <- cov_muni %>%
  filter(level_4 %in% ag_lu) %>%
  mutate(level_4 = ifelse(level_4 %in% natural_lu, "Native Vegetation", level_4)) %>%
  mutate(change = `2021`- `2000`,
         direction = ifelse(change <= 0, "contraction", "expansion")) %>%
  filter(direction == "expansion") %>%
  group_by(city) %>%
  mutate(total_change = sum(change, na.rm = T)) %>%
  group_by(city, level_4) %>%
  summarise(PAF = sum(change, na.rm = T) / unique(total_change)) %>%
  ungroup()
paa_muni_2000 <- 
  cov_muni %>%
  filter(level_4 %in% natural_lu) %>%
  mutate(change = `2000` - `2021`,
         direction = ifelse(change > 0, "contraction", "expansion")) %>%
  filter(direction == "contraction") %>%
  group_by(city) %>%
  summarise(native_veg_change = sum(change, na.rm = T)) %>%
  ungroup() %>%
  left_join(paf_muni_2000, by = "city") %>%
  mutate(PAA = native_veg_change * PAF) 


# STATE
paf_state_2000 <- cov_muni %>%
  mutate(region = str_split(city, pattern = " - "),
         state = map_chr(region,2)) %>%
  filter(level_4 %in% ag_lu) %>%
  mutate(level_4 = ifelse(level_4 %in% natural_lu, "Native Vegetation", level_4)) %>%
  group_by(state, level_4) %>%
  summarise(`2021` = sum(`2021`, na.rm = T),
            `2000` = sum(`2000`, na.rm = T)) %>%
  ungroup() %>%
  mutate(change = `2021`- `2000`,
         direction = ifelse(change <= 0, "contraction", "expansion")) %>%
  filter(direction == "expansion") %>%
  group_by(state) %>%
  mutate(total_change = sum(change, na.rm = T)) %>%
  group_by(state, level_4) %>%
  summarise(PAF = sum(change, na.rm = T) / unique(total_change)) %>%
  ungroup()
paa_state_2000 <- 
  cov_muni %>%
  mutate(region = str_split(city, pattern = " - "),
         state = map_chr(region,2)) %>%
  filter(level_4 %in% natural_lu) %>%
  mutate(change = `2000` - `2021`,
         direction = ifelse(change > 0, "contraction", "expansion")) %>%
  filter(direction == "contraction") %>%
  group_by(state) %>%
  summarise(native_veg_change = sum(change, na.rm = T)) %>%
  ungroup() %>%
  left_join(paf_state_2000, by = "state") %>%
  mutate(PAA = native_veg_change * PAF) 


# NATIONAL
paf_national_2000 <- cov_muni %>%
  filter(level_4 %in% ag_lu) %>%
  mutate(level_4 = ifelse(level_4 %in% natural_lu, "Native Vegetation", level_4)) %>%
  group_by(level_4) %>%
  summarise(`2021` = sum(`2021`, na.rm = T),
            `2000` = sum(`2000`, na.rm = T)) %>%
  ungroup() %>%
  mutate(change = `2021`- `2000`,
         direction = ifelse(change <= 0, "contraction", "expansion")) %>%
  filter(direction == "expansion") %>%
  mutate(total_change = sum(change, na.rm = T)) %>%
  group_by(level_4) %>%
  summarise(PAF = sum(change, na.rm = T) / unique(total_change)) %>%
  ungroup() %>%
  mutate(resolution = "national")
paa_national_2000 <- 
  cov_muni %>%
  filter(level_4 %in% natural_lu) %>%
  mutate(change = `2000` - `2021`,
         direction = ifelse(change > 0, "contraction", "expansion")) %>%
  filter(direction == "contraction") %>%
  summarise(native_veg_change = sum(change, na.rm = T)) %>%
  mutate(resolution = "national") %>%
  left_join(paf_national_2000, by = "resolution") %>%
  mutate(PAA = native_veg_change * PAF)  %>%
  select(-resolution)


# Here simplifying the original land use 
pan_2000 <- 
  bind_rows(
    "National" = paa_national_2000,
    "State" = paa_state_2000, 
    "Municipality" = paa_muni_2000,
    "Property" = paa_property_2000,
    "Pixel" = paa_pixel_2000,
    .id = "resolution"
  ) %>%
  mutate(resolution = as_factor(resolution), 
         resolution = fct_relevel(resolution, levels = c("National", "State", "Municipality", "Property","Pixel"))) %>%
  group_by(resolution) %>%
  mutate(total_paa = sum(PAA),
         level_4 = ifelse(level_4 %in% c("Aquaculture", "Citrus", "Mining"), "Other", level_4)) %>%
  group_by(resolution) %>%
  mutate(total_paa = sum(PAA, na.rm = T)) %>%
  group_by(level_4, resolution) %>%
  summarise(sum_paa = sum(PAA), 
            PAN = sum_paa / unique(total_paa) * 100) %>%
  ungroup()
p_pan_2000 <- pan_2000 %>%
  ggplot(aes(resolution, PAN, fill = level_4)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = cbPalette) +
  theme_minimal() +
  theme(legend.title = element_blank(), 
        text=element_text(size=16,  family= "DM Sans"), 
        legend.position = "none") +
  labs(x = "Scale of analysis", 
       y = "% deforestation attributed to each driver\nin Brazil (2000-2021)") +
  coord_flip()
p_pan_2000
```

## Figure 1

``` {r, fig.width = 6, fig.height = 14}
# Group legends
# ... simplify axis titles
# ... add A, B, C
# ... export/save a long, thin version
fig1 <- p_line / p_pan_1985 / p_pan_2000 + plot_annotation(tag_levels = 'A')
fig1
ggsave(filename = "fig1.tiff",
       plot = fig1,
       width = 18,
       height = 22,
       units = "cm")
```

## Stats quoted

"the proportion of deforestation attributed to soy is 3-4 times greater when calculated at national- rather than pixel-level"
``` {r}
pan_1985 %>%
  filter(level_4 == "Soybeans")
pan_2000 %>%
  filter(level_4 == "Soybeans")
```

"the proportion of deforestation attributed to soy is almost double when analyzed at the property- rather than the pixel-level (20.3% vs 11.4%)." 
``` {r}
# PROPERTY
change_property_mt_2008 <- property_mt %>%
  rename(class_id = class) %>%
  left_join(class_dic, by = "class_id") %>%
  left_join(mt_sf, by = "gid") %>%
  left_join(biome, by = c("cd_mun" = "GEOCODMU")) %>%
  select(-`.geo`,-`system:index`) %>%
  filter(!level_4 %in% non_veg_lu) %>%
  mutate(level_4 = ifelse(level_4 %in% natural_lu, "Native Vegetation", level_4)) %>%
  spread(year, area, fill = 0) %>%
  group_by(gid, BIOME, level_4) %>%
  summarise(`2021` = sum(`2021`, na.rm = T), 
            `2008` = sum(`2008`, na.rm = T)) %>%
  ungroup() %>%
  mutate(
    change = `2021` - `2008`,
    direction = ifelse(change < 0, "contraction", "expansion")
  )  %>%
    filter(change != 0)
paf_property_mt_2008 <- change_property_mt_2008 %>%
  filter(direction == "expansion") %>%
  filter(level_4 %in% ag_lu) %>%
  group_by(gid) %>%
  mutate(total_change = sum(change, na.rm = T)) %>%
  group_by(gid, level_4) %>%
  summarise(PAF = sum(change, na.rm = T) / unique(total_change)) %>%
  ungroup()
paa_property_mt_2008 <- change_property_mt_2008 %>%
  filter(direction == "contraction") %>%
  filter(level_4  == "Native Vegetation") %>%
  group_by(BIOME, gid) %>%
  summarise(native_veg_change = sum(change, na.rm = T) * -1) %>%
  ungroup() %>%
  left_join(paf_property_mt_2008, by = "gid") %>%
  mutate(PAA = native_veg_change * PAF) %>%
  drop_na(PAA)


# PIXEL
change_pixel_mt_biome <- 
  trans_muni %>%
  mutate(region = str_split(city, pattern = " - "),
         city = map_chr(region,1),
         state = map_chr(region,2), 
         city = str_trans(city)) %>%
  filter(state == "MT") %>%
  left_join(geocode_dic, by = c("city" = "MUNICIPALITY", state = "STATE")) %>%
  left_join(biome, by = c("GEOCODE" = "GEOCODMU")) %>%
  filter(from_level_4 %in% natural_lu, 
         !to_level_4 %in% natural_lu, 
         from_level_4 != to_level_4) %>%
  filter(!to_level_4 %in% non_veg_lu) %>%
  group_by(BIOME, to_level_4) %>%
  summarise(PAA = sum(`2008-2021`, na.rm = T)) %>%
  ungroup() %>%
  rename(level_4 = to_level_4)


# Merge property/pixel-level data
pan_mt_2008 <- 
  bind_rows(
    "Property" = paa_property_mt_2008,
    "Pixel" = change_pixel_mt_biome,
    .id = "resolution"
  ) %>%
  group_by(resolution) %>%
  filter(BIOME == "AMAZONIA") %>%
  group_by(resolution) %>%
  mutate(total_paa = sum(PAA, na.rm = T)) %>%
  group_by(level_4, resolution) %>%
  summarise(sum_paa = sum(PAA), 
            PAN = sum_paa / unique(total_paa) * 100) %>%
  ungroup()


# Compare the area of soy expansion at different scales
pan_mt_2008 %>%
  filter(level_4 == "Soybeans")
```

## Supplementary material

``` {r}
# There are three approaches for sLUC calculation
# [A] Shared responsibility allocation
# [B] product expansion allocation
# [C] Spatially explicit sLUC
# Note: these are calculated based on consecutive years of data
# ... i.e. year y, and y-1.


# [A] Shared responsibility allocation
p_sr <- 
  cov_muni %>%
  filter(!level_4 %in% c(non_veg_lu, natural_lu)) %>%
  select(-feature_id, -city, -class_id, -level_0, -level_1, -level_2, -level_3, -color) %>%
  group_by(level_4) %>%
  summarise_all(sum, na.rm = T) %>%
  ungroup() %>%
  gather(year, value, `1985`:`2021`) %>%
  group_by(year) %>%
  mutate(total_value = sum(value, na.rm = T)) %>%
  group_by(year, level_4) %>%
  summarise(sum_value = sum(value, na.rm = T), 
            perc_value = sum_value / unique(total_value) * 100) %>%
  # group small flows
  # mutate(small_flows = ifelse(perc_value < 1.5, TRUE, FALSE), 
  #        level_4 = ifelse(small_flows == T, "Other", level_4)) %>%
  group_by(year) %>%
  mutate(total_value = sum(sum_value, na.rm = T)) %>%
  group_by(year, level_4, total_value) %>%
  summarise(value = sum(sum_value, na.rm = T), 
            perc_value = value/ unique(total_value) * 100) %>%
  ungroup() %>%
  mutate(year = as.numeric(year)) %>%
  ggplot(aes(year, perc_value, fill = level_4)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = cbPalette) +
  labs(title = "Shared responsibility") +
  theme_classic() +
  theme(legend.position = "none",
        text=element_text(size=16,  family= "DM Sans"),
        axis.title.x = element_blank())


# [B] Product expansion allocation
# ... first I disambiguate that the column names are characters, not numbers
# ... then create a function to calculate the expansion in any two given years
# ... then I loop across sequential pairs of years to calculate the result
cov_muni_simple <- cov_muni %>%
  select(-feature_id, class_id, -`level_0`, -`level_1`, -`level_2`, -level_3, -color)
pe <- function(y, y1){
  col_name1 = y
  col_name2 = y1
  col_name3 = paste0(y ,"-", y1)
  change_muni <- cov_muni_simple %>%
    mutate(change = get(!!col_name2) - get(!!col_name1),
           direction = ifelse(change < 0, "contraction", "expansion")) %>%
    group_by(city, direction) %>%
    mutate(total_change = sum(change, na.rm = T)) %>%
    group_by(city, direction, level_4) %>%
    summarise(sum_change = sum(change, na.rm = T),
              prop_change = sum_change / unique(total_change)) %>%
    ungroup()
  expansion_muni <- change_muni %>%
    drop_na(prop_change) %>%
    filter(direction == "expansion") %>%
    select(-prop_change, -direction)
  change_muni <-
    change_muni %>%
    drop_na(prop_change) %>%
    filter(direction == "contraction") %>%
    select(-sum_change, -direction) %>%
    left_join(expansion_muni, by = "city") %>%
    rename(from_level_4 = level_4.x, 
           to_level_4 = level_4.y) %>%
    mutate(transition = sum_change * prop_change) %>%
    group_by(from_level_4, to_level_4) %>%
    summarise(value = sum(transition)) %>%
    ungroup() %>%
    mutate(years = col_name3)
  return(change_muni)
}
timeframe <- c(1985:2021)
concat_pairs <- function(x) {
  # Concatenate paris of items in a list
  # from https://stackoverflow.com/questions/12866418/how-to-loop-through-an-array-in-pairs
  s <- seq(2, length(x), 2)
  paste(x[s-1], x[s], sep="-")
}
year_pairs_odd <- concat_pairs(timeframe)
year_pairs_even <- concat_pairs(timeframe[-1])
year_pairs <- c(year_pairs_odd, year_pairs_even)
year_pairs_df <-
  tibble("YEARS" = year_pairs) %>%
  mutate(YEARS = str_split(year_pairs, "-")) %>%
  mutate(y = map_chr(YEARS, 1),
         y1 = map_chr(YEARS, 2)) %>%
  select(-YEARS) %>%
  arrange(y)
pe_res <- pmap_df(
  .l = year_pairs_df, 
  .f = pe
)
rm(year_pairs_odd, year_pairs_even, year_pairs, year_pairs_df, pe)


# visualise 
# ... NB we limit analyses to (expanding) agricultural land uses
# ... i.e. land uses which produce products
p_pe <- pe_res %>%
  filter(!to_level_4 %in% c(non_veg_lu, natural_lu)) %>%
  group_by(years) %>%
  mutate(total_value = sum(value, na.rm = T)) %>%
  group_by(years, to_level_4) %>%
  summarise(sum_value = sum(value, na.rm = T), 
            perc_value = sum_value/ unique(total_value) * 100) %>%
  ungroup() %>%
  # mutate(transition = paste(from_level_4, "-", to_level_4)) %>%
  # mutate(small_flows = ifelse(perc_value < 1.5, TRUE, FALSE), 
  #        to_level_4 = ifelse(small_flows == T, "Other", to_level_4)) %>%
  group_by(years) %>%
  mutate(total_value = sum(sum_value, na.rm = T)) %>%
  group_by(years, to_level_4, total_value) %>%
  summarise(value = sum(sum_value, na.rm = T), 
            perc_value = value/ unique(total_value) * 100) %>%
  ungroup() %>%
  mutate(years = str_sub(years, start = 6, end = 9),
         years = as.numeric(years)) %>%
  ggplot(aes(years, perc_value, fill = to_level_4)) +
  geom_bar(stat = "identity") +
  theme(axis.text = element_text(angle = 45, hjust = 1), 
        legend.title = element_blank()) +
  scale_fill_manual(values = cbPalette) +
  labs(title = "Product expansion", 
       y = "%") +
  theme_classic() +
  theme(legend.position = "bottom",
        text=element_text(size=16,  family= "DM Sans"),
        axis.title.x = element_blank(), 
        legend.title = element_blank()) +
  guides(fill = guide_legend(nrow = 3))


# [C] Spatially explicit sLUC
p_sluc_annual <- 
  trans_muni %>% 
  filter(from_level_4 %in% natural_lu, 
         !to_level_4 %in% c(non_veg_lu, natural_lu), 
         from_level_4 != to_level_4) %>%
  select(to_level_4,
         `1985-1986`, `1986-1987`, `1987-1988`, `1988-1989`, 
         `1989-1990`, `1990-1991`, `1991-1992`, `1992-1993`,
         `1993-1994`, `1994-1995`, `1995-1996`, `1996-1997`,
         `1997-1998`, `1998-1999`, `1999-2000`, `2000-2001`, 
         `2001-2002`, `2002-2003`, `2003-2004`, `2004-2005`, 
         `2005-2006`, `2006-2007`, `2007-2008`, `2008-2009`,
         `2009-2010`, `2010-2011`, `2011-2012`, `2012-2013`,
         `2013-2014`, `2014-2015`, `2015-2016`, `2016-2017`,
         `2017-2018`, `2018-2019`, `2019-2020`, `2020-2021`) %>%
  group_by(to_level_4) %>%
  summarise_all(sum, na.rm = T) %>%
  ungroup() %>%
  gather(years, value, `1985-1986`:`2020-2021`) %>%
  group_by(years) %>%
  mutate(total_value = sum(value)) %>%
  group_by(years, to_level_4) %>%
  summarise(sum_value = sum(value), 
            perc_value = sum_value / unique(total_value) * 100) %>%
  # group small flows
  mutate(small_flows = ifelse(perc_value < 1.5, TRUE, FALSE), 
         to_level_4 = ifelse(small_flows == T, "Other", to_level_4)) %>%
  group_by(years) %>%
  mutate(total_value = sum(sum_value, na.rm = T)) %>%
  group_by(years, to_level_4, total_value) %>%
  summarise(value = sum(sum_value, na.rm = T), 
            perc_value = value/ unique(total_value) * 100) %>%
  ungroup() %>%
   mutate(years = str_sub(years, start = 6, end = 9),
         years = as.numeric(years)) %>%
  ggplot(aes(years, perc_value, fill = to_level_4)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = cbPalette) +
  labs(title = "Spatially explicit sLUC", 
       y = "%") +
  theme_classic() +
  theme(legend.position = "none",
        text=element_text(size=16,  family= "DM Sans"), 
        axis.title.x = element_blank())


# Composite plot
p_si <- p_sr / p_sluc_annual / p_pe
p_si
ggsave(filename = "si.tiff",
       plot = p_si,
       width = 18,
       height = 22,
       units = "cm")
```