---
title: "Ghana - Direct and indirect deforestation"
author: "Cecile Renier"
date: "09/02/2025"
output:
  html_document:
    keep_md: yes
  word_document: default
editor_options:
  chunk_output_type: inline
---

Script creating the output of the manuscript: Renier C, Addoah T, Guye V, Garrett R, Van den Broeck G, zu Ermgassen E K H J, Meyfroidt P, 2025. Direct and indirect deforestation for cocoa in the tropical moist forests of Ghana, Envir. Res. Food Systems.

Input data created with the Google Earth Engine script.

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = F, warning = F,  fig.align = 'center',fig.height=5, fig.width=7, dpi = 200) 
 

options(encoding = "utf8")
library(raster)
library(terra)
library(tidyverse)
library(knitr)
library(sf)
library(scales)
library(kableExtra)
library(ggplot2)
library(dplyr)

```


# Summary of the Methods

First, we quantified the direct deforestation and degradation footprint of cocoa, oil palm, rubber, coconut, mining, and urban sprawl using the most accurate and up-to-date satellite-based datasets available, thus covering the period 2000-2019 (section 2.2 in the manuscript). 

Second, with the same data we calculated their expansion outside what was forest in 2000 (section 2.3.1). 

Third, absent a comprehensive food crop satellite-based dataset, we used statistical records and a land balance approach to estimate how much of that expansion occurred on food crops, thus displacing food crops elsewhere, based on the relative importance of the food crop area in the region compared to other pools of land ‘available’ for expansion (section 2.3.2.a). 

Fourth, following the same approach, we estimated the area of food crops which expanded or were displaced onto forests (section 2.3.2.b-c). 

Fifth, we attributed the total deforestation and degradation of land used for food crops in 2019 to the other land uses based on their relative expansion, to calculate their indirect footprint on forests through food crop displacement (section 2.3.3) (see SI figure S1 for the detailed steps of the land balance). 

We performed the calculations at the region level, to capture the specific dynamics in each region, and then repeated the land balance and indirect deforestation attribution with data aggregated directly at the level of the whole study area. 



# LOAD DATA

```{r load}

# Load the .rdata file (produced from the script "R_script_load_data.R")
# >>> Specify here the directory where the data are located

load(here::here(".../input_data/GHA_LOAD_DATA.rdata"))

```


## Color scales

```{r color palettes}

def_scale <- c(
  "#EB5E00", "#cc79a7","#9F4500", 
  "#9BDBFF",
  "#338ABC",
  "#FFD700","#DAA520", 
  "#646464","#474747",
  "#00FFFF","#00B9B9"
)


def_scale_direct <- c(
  "#EB5E00",
  "#cc79a7",
  "#00b3ff",
  "#FFD700",
  "#646464",
  "#00FFFF"
)

def_scale_direct_ccmerged <- c(
  "#EB5E00",
  "#EB5E00",
  "#00b3ff",
  "#FFD700",
  "#646464",
  "#00FFFF"
)

fc_def_scale_levels <- c(
  "#00B9B9",
  "#474747",
  "#DAA520", 
  "#9F4500",
  "#9BDBFF"
)

def_scale_levels <- c(
  "#00FFFF",
  "#646464",
  "#FFD700",
  "#cc79a7",
  "#EB5E00",
  "#338ABC",
  "#00B9B9",
  "#474747",
  "#DAA520", 
  "#9F4500",
  "#9BDBFF"
)


def_scale_levels_ccmerged <- c(
  "#00FFFF",
  "#646464",
  "#FFD700",
  "#EB5E00",
  "#EB5E00",
  "#338ABC",
  "#00B9B9",
  "#474747",
  "#DAA520", 
  "#9F4500",
  "#9BDBFF"
)

exp_scale <- c(
  "#8da0cb",
  "#66c2a5",
  "#fc8d62"
)

lu2019_scale <- c(
  "#DDDDDD", "#88CCEE","#999933","#8FBC8F","#44AA99","#117733","#332288", "#CC6677","#882255", "#DDCC77","#A0522D" )

```


## 2019 land use areas

```{r}
# In regions with a negative "OTHERLAND" (i.e. the sum of land use areas > region area), we need to adjust the area. We do this by subtracting the area surplus to the "OTHER TMF" area (we assume that the food crop is produced on land with a cover detected as forest (which is <8m and/or <1ha)). For 2019 it is only the case for the Western Region (-107Kha).

overlap_fc_forest_2019 <- merged_landuses_rg_2019_notadjusted %>%
  mutate(NOT_LUi_HA = REGION_HA - LANDUSE_HA) %>% 
  left_join(crops_rg_2019, by = "REGION") %>% 
  mutate(OTHERLAND_HA = NOT_LUi_HA - CROP_HA) %>% 
  mutate(OVERLAP_FC_F_HA = case_when(
    OTHERLAND_HA < 0 ~ -OTHERLAND_HA, 
    TRUE ~ 0)) 

landuses_rg_2019 <- landuses_rg_2019_notadjusted %>% 
  left_join(overlap_fc_forest_2019 %>% select(REGION, OVERLAP_FC_F_HA), by = "REGION") %>% 
  mutate(LANDUSE_HA = case_when(
    LANDUSE == "OTHERTMF" ~ LANDUSE_HA - OVERLAP_FC_F_HA, 
    TRUE ~ LANDUSE_HA
  )) %>% select(-OVERLAP_FC_F_HA)
  
merged_landuses_rg_2019 <- landuses_rg_2019 %>% 
  group_by(REGION) %>% 
  reframe(LANDUSE_HA = sum(LANDUSE_HA)) %>% 
  left_join(gha_regions_2007, by = "REGION")

landuses_hfz_2019 <- landuses_rg_2019 %>% 
  group_by(LANDUSE) %>% 
  reframe(LANDUSE_HA = sum(LANDUSE_HA))
```

## 2000 land use areas

### Estimation of 2000 Cocoa area

```{r}
cocoa_tot_50y <- cocoa_rg_50y %>% 
  group_by(YEAR) %>% 
  reframe(TOTAL = sum(PRODUCTION_MT))
```

```{r}
# Calculate 2000 cocoa area and expansion using Ruf data and CIAT cocoa per district + Adding unproductive land (estimate Teal 2002)

# Percentage of variation 
perc_var <- 0.1

# Ruf data on yields in 2000/01
# Does not include yield for Volta region so we use the average of the other regions
cocoa_yield_ruf_2000 <- data.frame(
  REGION  = c("WESTERN", "BRONG AHAFO", "ASHANTI", "CENTRAL", "EASTERN", "VOLTA"),
  YIELD_MID = c(356, 484, 181, 259, 158, 313)) %>% 
  mutate(YIELD_LO = YIELD_MID - YIELD_MID*perc_var) %>% 
  mutate(YIELD_HI = YIELD_MID + YIELD_MID*perc_var) 

# Unproductive land 2001/02 as per Teal 2006
perc_unprod_cocoa_area_teal <- 0.236

# cocoa_ciat %>% filter(YEAR == 1999) %>% {sum(.$VALUE)}
# cocoa_ciat %>% filter(YEAR == 2000) %>% {sum(.$VALUE)}
# cocoa_ciat %>% filter(YEAR == 2001) %>% {sum(.$VALUE)}

```

```{r}
# Per district
# Average cocoa production per district in 1999-2001
cocoa_dst_2000 <- cocoa_ciat %>% 
  filter(YEAR %in% c(1999:2001)) %>%
  group_by(REGION, DISTRICT) %>% 
  reframe(PRODUCTION_MT = mean(VALUE), YEAR = "1999-2001") %>% 
  ungroup()

# Estimate of 2000 cocoa area (total)
cocoa_area_dst_2000 <- cocoa_dst_2000  %>% 
  left_join(cocoa_yield_ruf_2000, by = "REGION") %>% 
  mutate(AREA_HA_MID = ((PRODUCTION_MT*1000)/YIELD_MID)* (1+ perc_unprod_cocoa_area_teal)) %>% 
  mutate(AREA_HA_MIN = ((PRODUCTION_MT*1000)/YIELD_HI)* (1+ perc_unprod_cocoa_area_teal)) %>% 
  mutate(AREA_HA_MAX = ((PRODUCTION_MT*1000)/YIELD_LO)* (1+ perc_unprod_cocoa_area_teal))

# Total cocoa area in 2000
# sum(cocoa_area_dst_2000$AREA_HA_MID)
# sum(cocoa_area_dst_2000$AREA_HA_MIN)
# sum(cocoa_area_dst_2000$AREA_HA_MAX)
```

```{r}
# Per region
# Estimate of 2000 cocoa area (total)
cocoa_area_rg_2000 <- cocoa_dst_2000 %>% 
  group_by(REGION) %>% 
  reframe(PRODUCTION_MT = sum(PRODUCTION_MT), YEAR = "1999-2001") %>% 
  left_join(cocoa_yield_ruf_2000, by = "REGION") %>% 
  mutate(MID_EXP = (PRODUCTION_MT*1000)/YIELD_MID* (1+ perc_unprod_cocoa_area_teal)) %>% 
  # A high yield gives a smaller area in 2000 => thus a higher expansion area between 2000-2019
  mutate(MAX_EXP = (PRODUCTION_MT*1000)/YIELD_HI* (1+ perc_unprod_cocoa_area_teal)) %>% 
  # A low yield gives a larger area in 2000 => thus a smaller expansion area between 2000-2019
  mutate(MIN_EXP = (PRODUCTION_MT*1000)/YIELD_LO* (1+ perc_unprod_cocoa_area_teal)) %>% 
  pivot_longer(cols = c("MID_EXP", "MIN_EXP", "MAX_EXP"),
               names_to = "SCENARIO",
               values_to = "AREA_HA")

# Total cocoa area in 2000
# cocoa_area_rg_2000 %>% group_by(SCENARIO) %>% reframe(AREA_HA = sum(AREA_HA))
# sum(cocoa_area_rg_2019$LANDUSE_HA)

cocoa_area_rg_2000 <- cocoa_area_rg_2000 %>% st_drop_geometry() %>% reframe(REGION, AREA_HA, SCENARIO, YEAR)
temp2 <- cocoa_area_rg_2019 %>% 
  reframe(REGION, AREA_HA = LANDUSE_HA, SCENARIO = "MID_EXP", YEAR = "2019") 
temp3 <- cocoa_area_rg_2000 %>% rbind(temp2)

ggplot(data = temp3 %>% 
         pivot_wider(names_from = SCENARIO, values_from = AREA_HA), 
       aes(x = YEAR, y = MID_EXP, fill = YEAR)) + 
    geom_bar(stat="identity", position = position_dodge(width=0.9)) +
    geom_errorbar(aes(ymin = MAX_EXP, ymax = MIN_EXP), width = 0.2) +
  facet_wrap(~REGION) +
    theme(axis.title.y = element_text(angle = 90, size=9),
        axis.text.y = element_text(size = 9),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(), 
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 9),
        strip.background = element_rect(colour = "grey50"),
        legend.position = "bottom", 
        legend.text = element_text(size = 9), 
        legend.title = element_text(size = 9),
        panel.spacing = unit(2, "lines"), 
        panel.background = element_rect(fill = "white", colour = "grey50"), 
        panel.grid.major.y = element_line(colour = "grey50"), 
        panel.grid.major.x = element_blank()) +
  labs(y = "Area (ha)", x = "Cocoa campaigns") +
  #scale_y_continuous(label = unit_format(unit = "K", scale = 1e-3, sep = "")) +
  scale_color_brewer(palette = "Dark2", "Cocoa area", labels=c('2000 (estimated)', '2019')) 
```

### Merge land uses

```{r}
# Negative "OTHERLAND" adjustment.
# In 2000, no region has an area surplus (though the code here runs the adjustment anyway in case input data change)

overlap_fc_forest <- merged_landuses_rg_2000_notadjusted %>%
  left_join(cocoa_area_rg_2000 %>% select(REGION, SCENARIO, COCOA_HA = AREA_HA), by = c("REGION", "SCENARIO")) %>% 
  group_by(REGION, SCENARIO) %>% 
  mutate(NOT_LUi_HA = REGION_HA - LANDUSE_HA_2000 - COCOA_HA) %>% 
  left_join(crops_rg_2000, by = "REGION") %>% 
  group_by(REGION, SCENARIO) %>% 
  mutate(OTHERLAND_HA = NOT_LUi_HA - CROP_HA) %>% 
  mutate(OVERLAP_FC_F_HA = case_when(
    OTHERLAND_HA < 0 ~ -OTHERLAND_HA, 
    TRUE ~ 0)) 

landuses_rg_2000 <- landuses_rg_2000_notadjusted %>% 
  left_join(overlap_fc_forest %>% 
              select(REGION, SCENARIO, OVERLAP_FC_F_HA), by = c("REGION", "SCENARIO")) %>% 
  mutate(LANDUSE_HA_2000 = case_when(
    LANDUSE == "UTMF" ~ LANDUSE_HA_2000 - OVERLAP_FC_F_HA, 
    TRUE ~ LANDUSE_HA_2000
  )) %>% select(-OVERLAP_FC_F_HA)
  
merged_landuses_rg_2000 <- landuses_rg_2000 %>% 
  group_by(REGION, SCENARIO) %>% 
  reframe(LANDUSE_HA_2000 = sum(LANDUSE_HA_2000)) %>% 
  left_join(gha_regions_2007, by = "REGION")

landuses_hfz_2000 <- landuses_rg_2000 %>% 
  group_by(LANDUSE, SCENARIO) %>% 
  reframe(LANDUSE_HA_2000 = sum(LANDUSE_HA_2000))

# Merge at regional level
merged_landuses_rg_2000_step2 <- merged_landuses_rg_2000 %>%
  left_join(cocoa_area_rg_2000 %>% select(REGION, COCOA_HA = AREA_HA, SCENARIO), by = c("REGION", "SCENARIO")) %>% 
  group_by(REGION, SCENARIO) %>% 
  mutate(NOT_LUi_HA = REGION_HA - LANDUSE_HA_2000 - COCOA_HA)

merged_landuses_rg_2000_step3 <- merged_landuses_rg_2000_step2 %>% 
  left_join(crops_rg_2000, by = "REGION") %>% 
  group_by(REGION, SCENARIO) %>% 
  mutate(OTHERLAND_HA = NOT_LUi_HA - CROP_HA) %>% 
  left_join(gha_regions_2007basedon2019 %>% select(-REGION_HA), by = "REGION")

# Merge at HFZ level
merged_landuses_hfz_2000_step2 <- merged_landuses_rg_2000_step2 %>%
  mutate(HFZ = "HFZ") %>% 
  group_by(HFZ, SCENARIO) %>% 
  summarize(LANDUSE_HA_2000 = sum(LANDUSE_HA_2000), REGION_HA = sum(REGION_HA), COCOA_HA = sum(COCOA_HA), NOT_LUi_HA = REGION_HA - LANDUSE_HA_2000 - COCOA_HA)

merged_landuses_hfz_2000_step3 <- merged_landuses_hfz_2000_step2 %>% 
  left_join(crops_rg_2000 %>% mutate(HFZ = "HFZ") %>% group_by(HFZ) %>% summarize(CROP_HA = sum(CROP_HA)), by = "HFZ") %>% 
  group_by(HFZ, SCENARIO) %>% 
  mutate(OTHERLAND_HA = NOT_LUi_HA - CROP_HA) 
```

# SECTION 2.2. DIRECT DEFORESTATION AND DEGRADATION

Here we calculate the expansion area (EA) between 2000 and 2019 on the 2000 forests of each land use for which we have Earth Observations.

```{r}
# Undisturbed Forest cover 2000: 
utmf2000 <- landuses_hfz_2000 %>% 
  filter(LANDUSE == "UTMF") %>% 
  group_by(SCENARIO) %>% #, LANDUSE
  reframe(LANDUSE_HA_2000 = sum(LANDUSE_HA_2000)) %>% 
  filter(SCENARIO == "MID_EXP") %>% pull(LANDUSE_HA_2000)

# Undisturbed Forest cover 2000 in protected areas:
utmf2000_cropped_pa <- tmf_cropped_pa_2000 %>%
  mutate(AREA_HA = UTMF) %>% 
  reframe(AREA_HA = sum(AREA_HA)) %>% pull(AREA_HA)

UTMF_2000_cropped_pa_cat <- tmf_cropped_pa_2000 %>% 
  group_by(CATEGORY) %>% 
  summarize(UTMF = sum(UTMF), DTMF = sum(DTMF))

# Remaining UTMF in 2019:
leftUTMF <- landuses_hfz_2019 %>%  
  filter(LANDUSE == "UTMF") %>% 
  reframe(AREA_HA = sum(LANDUSE_HA)) %>%  
  pull(AREA_HA)

leftDTMF <- landuses_hfz_2019 %>%  
  filter(LANDUSE == "DTMF") %>% 
  reframe(AREA_HA = sum(LANDUSE_HA)) %>%  
  pull(AREA_HA)

# Remaining UTMF in 2019 in protected areas (= 887Kha, similar to uncropped PAs (889Kha)):
leftUTMF_cropped_pa_2019_notadjusted <- landuses_cropped_pa_2019_notadjusted %>% 
  filter(LANDUSE == "UTMF") %>% 
  reframe(AREA_HA = sum(AREA_HA)) %>% 
  pull(AREA_HA)

# remaining DEGRADED forest (restricted to >8m and >1ha) in cropped Protected areas = 143Kha
leftDTMF_cropped_pa_2019_notadjusted <- landuses_cropped_pa_2019_notadjusted %>% 
  filter(LANDUSE == "DTMF") %>%  
  summarize(AREA_HA = sum(AREA_HA)) %>% 
  pull(AREA_HA)

 # landuses_def_cropped_pa_cat %>% 
 #  group_by(CATEGORY) %>% 
 #  reframe(AREA_HA = sum(AREA_HA)) %>% 
 #  mutate(PERC = AREA_HA/sum(AREA_HA)*100) %>% 
 #  filter(CATEGORY == "National Park / Game Reserve") %>% 
 #  pull(PERC)

def_cat <- landuses_def_cropped_pa_cat %>% 
  group_by(CATEGORY) %>% 
  reframe(AREA_HA = sum(AREA_HA)) %>% 
  left_join(UTMF_2000_cropped_pa_cat, by = "CATEGORY") %>%
  mutate(PERC = AREA_HA/UTMF*100) 
```

Over the 2000-2019 period, we found that, in the entire high forest zone (HFZ) of Ghana, `r round(sum(landuses_def_rg$AREA_HA)/10^6,1)` Mha (`r round(sum(landuses_def_rg$AREA_HA)/utmf2000 *100, 1)`) % of the initial `r round(utmf2000/10^6, 1)` Mha of undisturbed tropical moist forest (TMF) was lost, either degraded or deforested, leaving the country with about `r round(leftUTMF/10^6, 0)` Mha of undisturbed TMF and `r round(leftDTMF/10^6, 1)` Mha of degraded TMF in 2019, representing respectively `r round(leftUTMF/10^6/sum(gha_regions_2007basedon2019$REGION_HA/10^6)*100, 0)`% and `r round(leftDTMF/10^6/sum(gha_regions_2007basedon2019$REGION_HA/10^6)*100, 0)`% of the HFZ.

Over 2000-2019, `r round((sum(landuses_def_rg$AREA_HA)-sum(landuses_def_cropped_pa_rg$AREA_HA))/sum(landuses_def_rg$AREA_HA)*100, 0)`% of the total deforestation and degradation in the high forest zone happened off-reserve, where 75% of the undisturbed forest of 2000 was deforested and 11% degraded. Gazetted forests in the HFZ lost `r round(sum(landuses_def_cropped_pa_rg$AREA_HA)/utmf2000_cropped_pa *100, 1)` % of their 2000 undisturbed forest cover (`r round(sum(landuses_def_cropped_pa_rg$AREA_HA), 0)` ha of the initial `r round(utmf2000_cropped_pa/10^6, 1)`  Mha, with 14% deforested and 12% degraded).

`r round(def_cat %>% filter(CATEGORY == "Forest Reserve") %>% pull(AREA_HA)/sum(def_cat$AREA_HA)*100,0)`% of this loss happened in forest reserves. Only `r round(def_cat %>% filter(CATEGORY == "National Park / Game Reserve") %>% pull(AREA_HA),0)` ha were lost in the 5 national parks and game reserves of the HFZ (i.e. `r round(def_cat %>% filter(CATEGORY == "National Park / Game Reserve") %>% pull(PERC),0)` % of their 2000 undisturbed forest).

In 2019, `r round(leftUTMF_cropped_pa_2019_notadjusted/leftUTMF*100,0)`% of Ghana’s remaining `r round(leftUTMF/10^6, 0)` Mha of undisturbed tropical moist forest were located in protected forests (93% in forest reserves). [NB: the amount of UTMF remaining in protected areas cropped to the HFZ is very similar the area in uncropped protected areas].



#### Table S7: Tropical Moist Forest (TMF) area on- and off-reserve in the high forest zone in 2000 and 2019

```{r}
utmf_2000_2019_on_off_pa <- data.frame(
  FOREST = "UTMF", 
  LOCATION = c("TOTAL", "ON RESERVE", "OFF RESERVE"), 
  AREA_2000_HA = c(utmf2000, utmf2000_cropped_pa, 
                   utmf2000 - utmf2000_cropped_pa),
  AREA_2019_HA = c(leftUTMF, 
                   leftUTMF_cropped_pa_2019_notadjusted, 
                   leftUTMF 
                   - leftUTMF_cropped_pa_2019_notadjusted)) 
  
dtmf_2000_2019_on_off_pa <- data.frame(
  FOREST = "DTMF", 
  LOCATION = c("TOTAL", "ON RESERVE", "OFF RESERVE"), 
  AREA_2000_HA = c(landuses_hfz_2000 %>% filter(LANDUSE == "DTMF") %>% 
                     filter(SCENARIO == "MID_EXP") %>% pull(LANDUSE_HA_2000),
                   tmf_cropped_pa_2000 %>% mutate(AREA_HA = DTMF) %>% 
                     reframe(AREA_HA = sum(AREA_HA)) %>% pull(AREA_HA), 
                   landuses_hfz_2000 %>% filter(LANDUSE == "DTMF") %>% 
                     filter(SCENARIO == "MID_EXP") %>% pull(LANDUSE_HA_2000)
                   - tmf_cropped_pa_2000 %>% mutate(AREA_HA = DTMF) %>% 
                     reframe(AREA_HA = sum(AREA_HA)) %>% pull(AREA_HA)),
  AREA_2019_HA = c(leftDTMF, 
                   leftDTMF_cropped_pa_2019_notadjusted, 
                   leftDTMF
                   - leftDTMF_cropped_pa_2019_notadjusted))

forest_2000_2019_on_off_pa <- utmf_2000_2019_on_off_pa %>% 
  rbind(dtmf_2000_2019_on_off_pa) 

forest_2000_2019_on_off_pa

# Remarks:
# DTMF in 2019 accounts for all degraded forest (includes the 27Kha of 2000 => But these were not included in the land balance analysis which only considers loss of UTMF 2000)
# Areas in protected areas are only for cropped protected areas (matched with HFZ) and were not adjusted for the land use overlap (assumption is that this forest & food crop overlap was not occurring in PAs) while the total area in 2000 and 2019 were adjusted (subtracting the assumed overlap from food crop to the UTMF area)

```

#### fig.5: on- and off- forest reserve deforestation and degradation per land use

```{r,fig.height=5, fig.width=13}

# Calculate direct def and deg off-reserves (subtracting on- to Total)
data_plot_def <- landuses_def_rg %>% 
         group_by(DRIVER) %>% 
         reframe(TOTAL = sum(AREA_HA)) %>% 
  left_join(landuses_def_cropped_pa_rg %>% 
              group_by(DRIVER) %>% 
              reframe(ON_RESERVE = sum(AREA_HA)), 
            by = "DRIVER") %>%
  mutate(OFF_RESERVE = TOTAL - ON_RESERVE) %>% 
  pivot_longer(cols = c(2:4),
                 names_to = "LOCATIONS", 
               values_to = "AREA_HA")


data_plot_def_arranged <- data_plot_def %>% 
  mutate(LANDUSE = case_when(
                grepl("COCOA", DRIVER) ~ "COCOA",
                grepl("MINING", DRIVER) ~ "MINING",
                grepl("TREEPLANTATION", DRIVER) ~ "TREE PLANTATION",
                grepl("SETTLEMENT", DRIVER) ~ "SETTLEMENT",
                grepl("OTHER", DRIVER) ~ "REMAINING\n[INCL.FOOD CROP\nAND LOGGING]",
                TRUE ~ "NA")) %>% 
  mutate(DRIVER = case_when(
                grepl("COCOA_DEF", DRIVER) ~ "COCOA DIRECT DEFORESTATION",
                grepl("COCOA_DEG", DRIVER) ~ "COCOA DIRECT DEGRADATION",
                grepl("MINING", DRIVER) ~ "MINING DIRECT DEFORESTATION",
                grepl("TREEPLANTATION", DRIVER) ~ "TREE PLANTATION DIRECT DEFORESTATION",
                grepl("SETTLEMENT", DRIVER) ~ "SETTLEMENT DIRECT DEFORESTATION",
                grepl("OTHER_DEF", DRIVER) ~ "REMAINING DEF. & DEG.\n[INCL. FOOD CROP AND LOGGING]",
                grepl("OTHER_DEG", DRIVER) ~ "REMAINING DEF. & DEG.\n[INCL. FOOD CROP AND LOGGING]",
                TRUE ~ "NA")) %>% 
  group_by(DRIVER, LOCATIONS) %>% 
  reframe(AREA_HA = sum(AREA_HA), LANDUSE) %>% 
  unique() %>% 
      mutate(DRIVER = factor(DRIVER, levels = c("COCOA DIRECT DEFORESTATION", "COCOA DIRECT DEGRADATION",  "REMAINING DEF. & DEG.\n[INCL. FOOD CROP AND LOGGING]", "MINING DIRECT DEFORESTATION", "SETTLEMENT DIRECT DEFORESTATION","TREE PLANTATION DIRECT DEFORESTATION"))) %>% 
      arrange(DRIVER) %>% 
    mutate(LANDUSE = factor(LANDUSE, levels = c("COCOA", "REMAINING\n[INCL.FOOD CROP\nAND LOGGING]", "MINING", "SETTLEMENT","TREE PLANTATION"))) %>% 
          arrange(LANDUSE) 

# New facet label names 
locations_labels <- c("TOTAL HFZ", "OFF-RESERVE", "IN GAZETTED AREAS")
names(locations_labels) <- c("TOTAL", "OFF_RESERVE", "ON_RESERVE")


# Plot
ggplot(data = data_plot_def_arranged %>% 
         mutate(LOCATIONS = factor(LOCATIONS, levels= c("TOTAL", "OFF_RESERVE", "ON_RESERVE"))),
       aes(x = LANDUSE, y = AREA_HA, fill = DRIVER)) + 
    geom_bar(stat="identity", position = "stack") +
  facet_wrap(~LOCATIONS, labeller = labeller(LOCATIONS = locations_labels)) +
    coord_flip() +
    theme(axis.title.y = element_text(angle = 90, size=14),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 14),
        axis.text.x = element_text(size = 10), 
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 12),
        strip.background = element_rect(colour = "grey50"),
        legend.position = "right", 
        legend.text = element_text(size = 9), 
        legend.title = element_blank(),
        panel.spacing = unit(1.5, "lines"), 
        panel.background = element_rect(fill = "white", colour = "grey50"), 
        panel.grid.major.y = element_blank(), 
        panel.grid.major.x = element_line(colour = "grey90")) +
   scale_fill_manual(values = def_scale_direct) +
  labs(y = "Area deforested and/or degraded per land use (ha) 2000-2019", x = "") +
  scale_y_continuous(label = unit_format(unit = "K", scale = 1e-3, sep = ""))

```

#### fig.5 per region

```{r,fig.height=12, fig.width=7}

# Calculate direct def and deg off-reserves (subtracting on- to Total)
data_plot_def_rg <- landuses_def_rg %>% rename(TOTAL = AREA_HA) %>% 
  left_join(landuses_def_cropped_pa_rg %>% rename(ON_RESERVE = AREA_HA), by = c("DRIVER", "REGION")) %>%
  mutate(OFF_RESERVE = TOTAL - ON_RESERVE) %>% 
  pivot_longer(cols = c(3:5),
                 names_to = "LOCATIONS", 
               values_to = "AREA_HA")


data_plot_def_rg_arranged <- data_plot_def_rg %>% 
  mutate(LANDUSE = case_when(
                grepl("COCOA", DRIVER) ~ "COCOA",
                grepl("MINING", DRIVER) ~ "MINING",
                grepl("TREEPLANTATION", DRIVER) ~ "TREE PLANTATION",
                grepl("SETTLEMENT", DRIVER) ~ "SETTLEMENT",
                grepl("OTHER", DRIVER) ~ "REMAINING\n[INCL.FOOD CROP\nAND LOGGING]",
                TRUE ~ "NA")) %>% 
  mutate(DRIVER = case_when(
                grepl("COCOA_DEF", DRIVER) ~ "COCOA DIRECT DEFORESTATION",
                grepl("COCOA_DEG", DRIVER) ~ "COCOA DIRECT DEGRADATION",
                grepl("MINING", DRIVER) ~ "MINING DIRECT DEFORESTATION",
                grepl("TREEPLANTATION", DRIVER) ~ "TREE PLANTATION DIRECT DEFORESTATION",
                grepl("SETTLEMENT", DRIVER) ~ "SETTLEMENT DIRECT DEFORESTATION",
                grepl("OTHER_DEF", DRIVER) ~ "REMAINING DEF. & DEG.\n[INCL. FOOD CROP AND LOGGING]",
                grepl("OTHER_DEG", DRIVER) ~ "REMAINING DEF. & DEG.\n[INCL. FOOD CROP AND LOGGING]",
                TRUE ~ "NA")) %>% 
      mutate(DRIVER = factor(DRIVER, levels = c("COCOA DIRECT DEFORESTATION", "COCOA DIRECT DEGRADATION",  "REMAINING DEF. & DEG.\n[INCL. FOOD CROP AND LOGGING]", "MINING DIRECT DEFORESTATION", "SETTLEMENT DIRECT DEFORESTATION","TREE PLANTATION DIRECT DEFORESTATION"))) %>% 
      arrange(DRIVER) %>% 
    mutate(LANDUSE = factor(LANDUSE, levels = c("COCOA", "REMAINING\n[INCL.FOOD CROP\nAND LOGGING]", "MINING", "SETTLEMENT","TREE PLANTATION"))) %>% 
          arrange(LANDUSE) 

# Plot
ggplot(data = data_plot_def_rg_arranged %>% 
         mutate(LOCATIONS = factor(LOCATIONS, levels= c("TOTAL", "OFF_RESERVE", "ON_RESERVE"))),
       aes(x = LANDUSE, y = AREA_HA, fill = DRIVER)) + 
    geom_bar(stat="identity", position = "stack") +
   facet_grid(REGION ~ LOCATIONS) + # labeller = labeller(LOCATIONS = locations_labels)
    coord_flip() +
    theme(axis.title.y = element_text(angle = 90, size=11),
        axis.text.y = element_text(size = 10),
        axis.title.x = element_text(size = 10),
        axis.text.x = element_text(size = 9), 
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 10),
        strip.background = element_rect(colour = "grey50"),
        legend.position = "bottom", 
        legend.text = element_text(size = 7), 
        legend.title = element_blank(),
        panel.spacing = unit(1.5, "lines"), 
        panel.background = element_rect(fill = "white", colour = "grey50"), 
        panel.grid.major.y = element_blank(), 
        panel.grid.major.x = element_line(colour = "grey90")) +
   scale_fill_manual(values = def_scale_direct) +
  labs(y = "Area deforested and/or degraded per land use (ha) 2000-2019", x = "") +
  scale_y_continuous(label = unit_format(unit = "K", scale = 1e-3, sep = ""))

```

#### Table S9: on- and off- forest reserve deforestation and degradation per land use

```{r}
forest_2000_2019_on_off_pa_per_lu <- data_plot_def_arranged %>% 
  pivot_wider(names_from = "LOCATIONS",
              values_from = "AREA_HA") %>% 
  select(-LANDUSE) %>% 
  mutate(PERC_OFF = OFF_RESERVE/TOTAL*100, PERC_ON = ON_RESERVE/TOTAL*100) %>% 
  reframe(DRIVER, OFF_RESERVE, PERC_OFF, ON_RESERVE, PERC_ON, TOTAL)

forest_2000_2019_on_off_pa_per_lu
```


# SECTION 2.3. QUANTIFYING INDIRECT DEFORESTATION AND DEGRADATION THROUGH FOOD CROP DISPLACEMENT

# SECTION 2.3.1. Estimating the expansion of land uses and food crops in total and outside forests

## SECTION 2.3.1.a. Total expansion area of land uses

```{r}
# Cocoa expansion area per region
cocoa_EA <- cocoa_area_rg_2000 %>% 
  reframe(REGION, YEAR1_HA = AREA_HA, SCENARIO) %>% 
  left_join(cocoa_area_rg_2019 %>% 
              reframe(REGION, YEAR2_HA = LANDUSE_HA), by = "REGION") %>% 
  mutate(EA_HA = YEAR2_HA - YEAR1_HA, LANDUSE = "COCOA")

cocoa_EA_hfz <- cocoa_EA %>% 
  group_by(SCENARIO) %>% 
  summarise(YEAR1_HA = sum(YEAR1_HA), YEAR2_HA = sum(YEAR2_HA), EA_HA = sum(EA_HA))
```

```{r}
# Merge all Land use expansion areas (calculated on GEE and R)

landuses_EA_allrg <- landuses_EA_gee %>% 
  group_by(LANDUSE, SCENARIO) %>% 
  reframe(EA_HA = sum(EA_HA)) %>% 
  rbind(cocoa_EA %>% 
          group_by(LANDUSE, SCENARIO) %>% 
          reframe(EA_HA = sum(EA_HA))) %>% 
    rbind(crops_EA %>% 
          summarise(EA_HA = sum(EA_HA), LANDUSE, SCENARIO = "MID_EXP") %>% unique())

landuses_EA <- landuses_EA_gee %>% 
  group_by(REGION, LANDUSE, SCENARIO) %>% 
  reframe(EA_HA = sum(EA_HA)) %>% 
  rbind(cocoa_EA %>% 
          summarise(REGION, EA_HA, LANDUSE, SCENARIO)) %>% 
    rbind(crops_EA %>% 
          summarise(REGION, EA_HA, LANDUSE, SCENARIO = "MID_EXP")) %>% 
    rbind(crops_EA %>% 
          summarise(REGION, EA_HA, LANDUSE, SCENARIO = "MIN_EXP")) %>% 
    rbind(crops_EA %>% 
          summarise(REGION, EA_HA, LANDUSE, SCENARIO = "MAX_EXP")) 

# Adding the percentage of expansion ONLY FOR POSITIVE EXPANSION !
landuses_EA_perc <- landuses_EA %>% 
  mutate(EA_HA_raw = EA_HA) %>% 
  mutate(EA_HA = ifelse(EA_HA < 0, 0, EA_HA)) %>% 
  group_by(REGION, SCENARIO) %>% 
  mutate(TOT_EA = sum(EA_HA)) %>% 
  ungroup() %>% 
  mutate(PERC_EA = (EA_HA / TOT_EA) *100) 


```


## SECTION 2.3.1.b. Expansion area outside forests

### a) Per region

To attribute foodcrop deforestation to LULC_i expansion, we need to remove all the expansion area that occurred on land that was forest in 2000.

```{r}
# Total LU expansion areas OUTSIDE UTMF 2005
# = EA - Direct deforestation 

landuses_EA_outUTMF <- landuses_EA %>%  
  left_join(landuses_def_rg %>% 
              mutate(DRIVER = ifelse(grepl("COCOA", DRIVER), "COCOA_DD_HA", DRIVER)) %>% 
              group_by(REGION, DRIVER) %>% 
              mutate(AREA_HA = sum(AREA_HA)) %>% 
              unique() %>% 
              mutate(LANDUSE = case_when(
                grepl("COCOA", DRIVER) ~ "COCOA",
                grepl("MINING", DRIVER) ~ "MINING",
                grepl("TREEPLANTATION", DRIVER) ~ "TREEPLANTATION",
                grepl("SETTLEMENT", DRIVER) ~ "SETTLEMENT",
                grepl("OVERLAP", DRIVER) ~ "OVERLAP",
                grepl("OTHER", DRIVER) ~ "OTHER",
                TRUE ~ "NA"
              )), by =c("REGION", "LANDUSE")) %>%  
  group_by(REGION, LANDUSE, SCENARIO) %>% 
  mutate(EA_OUT_HA = EA_HA - AREA_HA) %>% 
   filter(LANDUSE != "CROP")

```


```{r}
# Table
# Expansion Areas in total in cocoa growing regions 
landuses_EA_outUTMF_tot <- landuses_EA_outUTMF %>% 
  group_by(DRIVER, SCENARIO) %>% 
  reframe(EA_OUT_HA = sum(EA_OUT_HA))

# landuses_EA_outUTMF_tot
```


### Summary table: total, outside- and on-forests expansion area per land use and per region

```{r}
# Plot EA aggregated over HFZ

landuses_EA_tot_out_on_UTMF_rg <- landuses_EA_outUTMF %>% 
  filter(SCENARIO == "MID_EXP") %>% 
         pivot_longer(cols = c(EA_HA, AREA_HA, EA_OUT_HA), 
               names_to = "METRIC", 
               values_to = "AREA_HA_MID") %>% 
  reframe(REGION, LANDUSE, METRIC, AREA_HA_MID) %>% 

  left_join(landuses_EA_outUTMF %>% 
  filter(SCENARIO == "MIN_EXP") %>% 
         pivot_longer(cols = c(EA_HA, AREA_HA, EA_OUT_HA), 
               names_to = "METRIC", 
               values_to = "AREA_HA_MIN") %>% 
  reframe(REGION, LANDUSE, METRIC, AREA_HA_MIN), by = c("REGION", "METRIC", "LANDUSE")) %>% 

  left_join(landuses_EA_outUTMF %>% 
  filter(SCENARIO == "MAX_EXP") %>% 
         pivot_longer(cols = c(EA_HA, AREA_HA, EA_OUT_HA), 
               names_to = "METRIC", 
               values_to = "AREA_HA_MAX") %>% 
  reframe(REGION, LANDUSE, METRIC, AREA_HA_MAX), by = c("REGION", "METRIC", "LANDUSE")) 
```


```{r}
# Expansion of cocoa in Western region percentage over total LU expansion

# landuses_EA_tot_out_on_UTMF_rg %>% 
#   filter(REGION == "WESTERN") %>% 
#   filter(METRIC == "EA_HA") %>% 
#   select(REGION, LANDUSE, AREA_HA_MID) %>% 
#   rbind(crops_EA %>% 
#           filter(REGION == "WESTERN") %>% 
#            select(REGION, LANDUSE, AREA_HA_MID = EA_HA)) %>% 
#   group_by(REGION) %>% 
#   mutate(TOT = sum(AREA_HA_MID)) %>% 
#   mutate(PERC = AREA_HA_MID / TOT)
  
```

#### fig.S5: Land use expansion over 2000-2019 per land use (aggregated at HFZ level)

```{r,fig.height=4, fig.width=7}

# Regional results aggregated at the HFZ level

ggplot(data = landuses_EA_tot_out_on_UTMF_rg %>% 
         group_by(LANDUSE, METRIC) %>% 
         reframe(AREA_HA_MID = sum(AREA_HA_MID), AREA_HA_MIN = sum(AREA_HA_MIN), AREA_HA_MAX = sum(AREA_HA_MAX)) %>% 
               mutate(METRIC = factor(METRIC, levels= c("EA_HA", "AREA_HA", "EA_OUT_HA"))),
       aes(x = METRIC, y = AREA_HA_MID, fill = METRIC))  + 
  geom_bar(stat="identity", width = 0.95) +
  geom_errorbar(aes(ymin = AREA_HA_MIN, ymax = AREA_HA_MAX), width = 0.2) +
  facet_grid(~LANDUSE, switch = "x") +
    theme(axis.title.y = element_text(angle = 90, size=15),
        axis.text.y = element_text(size = 14),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(), 
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 10),
        strip.background = element_blank(),
        legend.position = "right", 
        legend.text = element_text(size = 15), 
        legend.title = element_text(size = 15),
        panel.spacing = unit(0.2, "lines"), 
        panel.background = element_rect(fill = "grey95"), 
        panel.grid.major.y = element_line(colour = "white"),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank()) +
  scale_fill_manual(values = exp_scale, labels=c("Total expansion", "Expansion in forest", "Expansion outside forest")) + #  
  labs(y = "Area (ha)", x = "", fill = "Land use expansion\n2000-2019") +
  scale_y_continuous(label = unit_format(unit = "K", scale = 1e-3, sep = ""))

```

#### fig.S5: Land use expansion over 2000-2019 per land use and per region

```{r,fig.height=7, fig.width=7}
# Plot EA per region

ggplot(data = landuses_EA_tot_out_on_UTMF_rg %>% 
               mutate(METRIC = factor(METRIC, levels= c("EA_HA", "AREA_HA", "EA_OUT_HA"))), aes(x = METRIC, y = AREA_HA_MID, fill = METRIC))  + 
    geom_bar(stat="identity", width = 0.95) +
   geom_errorbar(aes(ymin = AREA_HA_MIN, ymax = AREA_HA_MAX), width = 0.2) +
  facet_grid(rows = vars(REGION), cols = vars(LANDUSE)) +
    theme(axis.title.y = element_text(angle = 90, size=15),
        axis.text.y = element_text(size = 11),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(), 
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 13),
        strip.text.y = element_text(size = 13),
        strip.background = element_blank(),
        legend.position = "below", 
        legend.text = element_text(size = 11), 
        legend.title = element_text(size = 11),
        panel.spacing = unit(0.4, "lines"), 
        panel.background = element_rect(fill = "grey95"), 
        panel.grid.major.y = element_line(colour = "white"),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank()) +
   scale_fill_manual(values = exp_scale, labels=c("Total expansion", "Expansion on forest", "Expansion outside forest")) +
  labs(y = "Area (ha)", x = "", fill = "2000-2019") +
  scale_y_continuous(label = unit_format(unit = "K", scale = 1e-3, sep = ""))

```


### b) On all HFZ

```{r,fig.height=4, fig.width=7}
# Total LU expansion areas OUTSIDE UTMF 2000
# = EA - Direct deforestation 

landuses_EA_outUTMF_hfz <- landuses_EA %>%  
  left_join(landuses_def_rg %>% 
              mutate(DRIVER = ifelse(grepl("COCOA", DRIVER), "COCOA_DD_HA", DRIVER)) %>% 
              group_by(REGION, DRIVER) %>% 
              mutate(AREA_HA = sum(AREA_HA)) %>% 
              unique() %>% 
              mutate(LANDUSE = case_when(
                grepl("COCOA", DRIVER) ~ "COCOA",
                grepl("MINING", DRIVER) ~ "MINING",
                grepl("TREEPLANTATION", DRIVER) ~ "TREEPLANTATION",
                grepl("SETTLEMENT", DRIVER) ~ "SETTLEMENT",
                grepl("OVERLAP", DRIVER) ~ "OVERLAP",
                grepl("OTHER", DRIVER) ~ "OTHER",
                TRUE ~ "NA"
              )), by =c("REGION", "LANDUSE")) %>%  
  select(-DRIVER) %>% 
  group_by(LANDUSE, SCENARIO) %>% 
  mutate(EA_HA = sum(EA_HA), AREA_HA = sum(AREA_HA)) %>% 
  mutate(EA_OUT_HA = EA_HA - AREA_HA, HFZ = "HFZ") %>% 
  select(-REGION) %>% 
filter(LANDUSE != "CROP") %>% 
  unique()

```



```{r}
#### fig. S5: Land use expansion over 2000-2019 per land use (calculated at HFZ level) = exactly the same as calculated at regional level

# ggplot(data = landuses_EA_outUTMF_hfz %>% 
#   pivot_longer(cols = c(EA_HA, AREA_HA, EA_OUT_HA), 
#                names_to = "METRIC", 
#                values_to = "AREA_HA") %>% 
#   pivot_wider( names_from = "SCENARIO", 
#                values_from = "AREA_HA") %>% 
#    mutate(METRIC = factor(METRIC, levels= c("EA_HA", "AREA_HA", "EA_OUT_HA"))), 
#   aes(x = METRIC, y = MID_EXP, fill = METRIC)) +
#   geom_bar(stat="identity", width = 0.95) +
#   geom_errorbar(aes(ymin = MIN_EXP, ymax = MAX_EXP), width = 0.2) +
#   facet_grid(~LANDUSE, switch = "x") +
#     theme(axis.title.y = element_text(angle = 90, size=15),
#         axis.text.y = element_text(size = 14),
#         axis.title.x = element_blank(),
#         axis.text.x = element_blank(), 
#         axis.ticks.x = element_blank(),
#         strip.text.x = element_text(size = 10),
#         strip.background = element_blank(),
#         legend.position = "right", 
#         legend.text = element_text(size = 15), 
#         legend.title = element_text(size = 15),
#         panel.spacing = unit(0.2, "lines"), 
#         panel.background = element_rect(fill = "grey95"), 
#         panel.grid.major.y = element_line(colour = "white"),
#         panel.grid.minor.y = element_blank(),
#         panel.grid.major.x = element_blank()) +
#   scale_fill_manual(values = exp_scale,  labels=c("Total expansion", "Expansion in forest", "Expansion outside forest")) +
#   labs(y = "Area (ha)", x = "", fill= "Land use expansion\n2000-2019") +
#   scale_y_continuous(label = unit_format(unit = "K", scale = 1e-3, sep = ""))

```

# SECTION 2.3.2. Estimating food crop displacement

## SECTION 2.3.2.a	Expansion of land uses on each land pool outside forests, including on food crops 

### a) Per region

### Area of converted cocoa per region

```{r}
# Where expansion of cocoa outside forest is negative = converted cocoa 

converted_cocoa_rg <- landuses_EA_outUTMF %>% 
  filter(LANDUSE == "COCOA") %>% 
  mutate(EA_OUT_HA = ifelse(EA_OUT_HA < 0, EA_OUT_HA, 0)) %>%
  group_by(REGION, SCENARIO) %>% 
  reframe(CONVERTED_CC_HA = sum(EA_OUT_HA)) 
```

### Total 'land pools' area 

```{r}
# For mining, tree plantations and settlement
# Denominateur: Calculate the total area of Other land + Food crop 2000 and Converted Cocoa 

land_for_expansion <- merged_landuses_rg_2000_step3 %>% 
  select(REGION, CROP_HA, OTHERLAND_HA, SCENARIO) %>% 
  left_join(converted_cocoa_rg, by = c("REGION", "SCENARIO")) %>% 
  group_by(REGION, SCENARIO) %>% 
  mutate(LAND_FOR_LUi_EXPANSION = CROP_HA + OTHERLAND_HA - CONVERTED_CC_HA)

```

```{r}
# For cocoa
# Denominateur: Calculate the total area of Other land + Food crop 2000 

land_for_cc_expansion <- merged_landuses_rg_2000_step3 %>% 
  select(REGION, SCENARIO, CROP_HA, OTHERLAND_HA) %>% 
  group_by(REGION, SCENARIO) %>% 
  mutate(LAND_FOR_CC_EXPANSION = CROP_HA + OTHERLAND_HA)

```

### Land balance: Distributing land use expansion according to the relative pool areas

```{r}
# For mining, tree plantations and settlement
# EA of all merged LUi on other land / Food crops / converted cocoa

ea_lui_off <- landuses_EA_outUTMF %>% 
  filter(LANDUSE != "COCOA") %>% 
  group_by(REGION, SCENARIO) %>% 
  reframe(EA_LUi_OFF_HA = sum(EA_OUT_HA)) %>% 
  left_join(land_for_expansion, by = c("REGION", "SCENARIO")) %>% 
  mutate(EA_LUi_ON_OL = EA_LUi_OFF_HA * (OTHERLAND_HA / LAND_FOR_LUi_EXPANSION)) %>% 
  mutate(EA_LUi_ON_FC = EA_LUi_OFF_HA * (CROP_HA / LAND_FOR_LUi_EXPANSION)) %>% 
  mutate(EA_LUi_ON_CC = EA_LUi_OFF_HA * (-CONVERTED_CC_HA / LAND_FOR_LUi_EXPANSION)) %>% 
  # Check caps = ok for all
  mutate(CHECK_CAP_OL = ifelse(EA_LUi_ON_OL > OTHERLAND_HA, "STOP", "OK")) %>% 
  mutate(CHECK_CAP_FC = ifelse(EA_LUi_ON_FC > CROP_HA, "STOP", "OK")) %>% 
  mutate(CHECK_CAP_CC = ifelse(EA_LUi_ON_CC > -CONVERTED_CC_HA , "STOP", "OK"))  

# EA of each LUi on other land / Food crops / converted cocoa

ea_per_lui_off <- landuses_EA_outUTMF %>% 
  filter(LANDUSE != "COCOA") %>% 
  rename(EA_LUi_OFF_HA = EA_OUT_HA) %>% 
  left_join(land_for_expansion, by = c("REGION", "SCENARIO")) %>% 
  mutate(EA_LUi_ON_OL = EA_LUi_OFF_HA * (OTHERLAND_HA / LAND_FOR_LUi_EXPANSION)) %>% 
  mutate(EA_LUi_ON_FC = EA_LUi_OFF_HA * (CROP_HA / LAND_FOR_LUi_EXPANSION)) %>% 
  mutate(EA_LUi_ON_CC = EA_LUi_OFF_HA * (-CONVERTED_CC_HA / LAND_FOR_LUi_EXPANSION)) 
# %>% select(REGION, LANDUSE, EA_LUi_OFF_HA, EA_LUi_ON_OL, EA_LUi_ON_FC, EA_LUi_ON_CC)

# ea_per_lui_off

# Check caps : ok
# ea_per_lui_off %>% 
#   group_by(REGION, SCENARIO) %>% 
#   reframe(OTHERLAND_HA, CROP_HA, CONVERTED_CC_HA, EA_LUi_ON_OL = sum(EA_LUi_ON_OL), EA_LUi_ON_FC = sum(EA_LUi_ON_FC), EA_LUi_ON_CC = sum(EA_LUi_ON_CC)) %>% 
#   mutate(CHECK_CAP_OL = ifelse(EA_LUi_ON_OL > OTHERLAND_HA, "STOP", "OK")) %>% 
#   mutate(CHECK_CAP_FC = ifelse(EA_LUi_ON_FC > CROP_HA, "STOP", "OK")) %>% 
#   mutate(CHECK_CAP_CC = ifelse(EA_LUi_ON_CC > -CONVERTED_CC_HA , "STOP", "OK"))  

# Check if same as all_Lui : ok
# ea_per_lui_off %>% 
#   select(c(REGION, LANDUSE, SCENARIO, EA_LUi_ON_OL, EA_LUi_ON_FC, EA_LUi_ON_CC)) %>% 
#   pivot_longer(cols = c(EA_LUi_ON_OL, EA_LUi_ON_FC, EA_LUi_ON_CC), 
#                names_to = "EXPANSION", 
#                values_to = "AREA_HA") %>% 
#   group_by(REGION, SCENARIO, EXPANSION) %>% 
#   reframe(AREA_HA = sum(AREA_HA)) %>% 
#   left_join(ea_lui_off %>% 
#               select(REGION, SCENARIO, EA_LUi_ON_OL, EA_LUi_ON_FC, EA_LUi_ON_CC) %>% 
#               pivot_longer(cols = c(EA_LUi_ON_OL, EA_LUi_ON_FC, EA_LUi_ON_CC), 
#               names_to = "EXPANSION", 
#               values_to = "ALL_AREA_HA"), 
#             by = c("REGION", "EXPANSION", "SCENARIO")) %>%
#   group_by(REGION, EXPANSION, SCENARIO) %>% 
#   mutate(CHECK_SAME = AREA_HA - ALL_AREA_HA)

```

```{r}
# For cocoa
# EA of cocoa on other land / Food crops 

ea_cc_off <- landuses_EA_outUTMF %>% 
  filter(LANDUSE == "COCOA") %>% 
  mutate(EA_OUT_HA = ifelse(EA_OUT_HA < 0, 0, EA_OUT_HA)) %>% 
  left_join(land_for_cc_expansion, by = c("REGION", "SCENARIO")) %>% 
  group_by(REGION, SCENARIO) %>% 
  mutate(EA_CC_ON_OL = EA_OUT_HA * (OTHERLAND_HA / LAND_FOR_CC_EXPANSION)) %>% 
  mutate(EA_CC_ON_FC = EA_OUT_HA * (CROP_HA / LAND_FOR_CC_EXPANSION)) %>%  
  # Check caps = ok for all
  mutate(CHECK_CAP_OL = ifelse(EA_CC_ON_OL > OTHERLAND_HA, "STOP", "OK")) %>% 
  mutate(CHECK_CAP_FC = ifelse(EA_CC_ON_FC > CROP_HA, "STOP", "OK")) 

```


### b) On all HFZ

### Area of converted cocoa

```{r}
# Where expansion of cocoa outside forest is negative = converted cocoa 

converted_cocoa_hfz <- landuses_EA_outUTMF_hfz %>% 
  filter(LANDUSE == "COCOA") %>% 
  mutate(EA_OUT_HA = ifelse(EA_OUT_HA < 0, EA_OUT_HA, 0)) %>% 
  group_by(LANDUSE, SCENARIO) %>% 
  reframe(HFZ, CONVERTED_CC_HA = sum(EA_OUT_HA)) %>% 
  select(-LANDUSE)

```

### Total 'land pools' area

```{r}
# For mining, tree plantations, settlement
# Denominateur: Calculate the total area of Other land + Food crop 2000 and Converted Cocoa 

land_for_expansion_hfz <- merged_landuses_hfz_2000_step3 %>% 
  select(HFZ, CROP_HA, OTHERLAND_HA, SCENARIO) %>% 
  left_join(converted_cocoa_hfz, by = c("HFZ", "SCENARIO")) %>% 
  group_by(HFZ, SCENARIO) %>% 
  mutate(LAND_FOR_LUi_EXPANSION = CROP_HA + OTHERLAND_HA - CONVERTED_CC_HA)

```

```{r}
# For cocoa
# Denominateur: Calculate the total area of Other land + Food crop 2000 

land_for_cc_expansion_hfz <- merged_landuses_hfz_2000_step3 %>%
  select(HFZ, CROP_HA, OTHERLAND_HA, SCENARIO) %>%
  group_by(HFZ, SCENARIO) %>% 
  mutate(LAND_FOR_CC_EXPANSION = CROP_HA + OTHERLAND_HA)

```

### Land balance: Distributing land use expansion according to the relative pool areas


```{r}
# For mining, tree plantations, settlement
# EA of all LUi on other land / Food crops / converted cocoa

ea_lui_off_hfz <- landuses_EA_outUTMF_hfz %>% 
  filter(LANDUSE != "COCOA") %>% 
  group_by(HFZ, SCENARIO) %>% 
  reframe(EA_LUi_OFF_HA = sum(EA_OUT_HA)) %>% 
  left_join(land_for_expansion_hfz, by = c("HFZ", "SCENARIO")) %>% 
  mutate(EA_LUi_ON_OL = EA_LUi_OFF_HA * (OTHERLAND_HA / LAND_FOR_LUi_EXPANSION)) %>% 
  mutate(EA_LUi_ON_FC = EA_LUi_OFF_HA * (CROP_HA / LAND_FOR_LUi_EXPANSION)) %>% 
  mutate(EA_LUi_ON_CC = EA_LUi_OFF_HA * (-CONVERTED_CC_HA / LAND_FOR_LUi_EXPANSION)) %>% 
  # Check caps = ok for all
  mutate(CHECK_CAP_OL = ifelse(EA_LUi_ON_OL > OTHERLAND_HA, "STOP", "OK")) %>% 
  mutate(CHECK_CAP_FC = ifelse(EA_LUi_ON_FC > CROP_HA, "STOP", "OK")) %>% 
  mutate(CHECK_CAP_CC = ifelse(EA_LUi_ON_CC > -CONVERTED_CC_HA , "STOP", "OK"))  


# EA of each LUi on other land / Food crops / converted cocoa

ea_per_lui_off_hfz <- landuses_EA_outUTMF_hfz %>% 
  filter(LANDUSE != "COCOA") %>% 
  rename(EA_LUi_OFF_HA = EA_OUT_HA) %>% 
  left_join(land_for_expansion_hfz, by = c("HFZ", "SCENARIO")) %>% 
  group_by(LANDUSE, SCENARIO) %>% 
  mutate(EA_LUi_ON_OL = EA_LUi_OFF_HA * (OTHERLAND_HA / LAND_FOR_LUi_EXPANSION)) %>% 
  mutate(EA_LUi_ON_FC = EA_LUi_OFF_HA * (CROP_HA / LAND_FOR_LUi_EXPANSION)) %>% 
  mutate(EA_LUi_ON_CC = EA_LUi_OFF_HA * (-CONVERTED_CC_HA / LAND_FOR_LUi_EXPANSION)) 
# %>% select(REGION, LANDUSE, EA_LUi_OFF_HA, EA_LUi_ON_OL, EA_LUi_ON_FC, EA_LUi_ON_CC)

# ea_per_lui_off_hfz

# Check caps : ok 
# ea_per_lui_off_hfz %>% 
#   group_by(HFZ, SCENARIO) %>% 
#   reframe(OTHERLAND_HA, CROP_HA, CONVERTED_CC_HA, EA_LUi_ON_OL = sum(EA_LUi_ON_OL), EA_LUi_ON_FC = sum(EA_LUi_ON_FC), EA_LUi_ON_CC = sum(EA_LUi_ON_CC)) %>% 
#   mutate(CHECK_CAP_OL = ifelse(EA_LUi_ON_OL > OTHERLAND_HA, "STOP", "OK")) %>% 
#   mutate(CHECK_CAP_FC = ifelse(EA_LUi_ON_FC > CROP_HA, "STOP", "OK")) %>% 
#   mutate(CHECK_CAP_CC = ifelse(EA_LUi_ON_CC > -CONVERTED_CC_HA , "STOP", "OK"))  

```

```{r}
# For cocoa
# EA of cocoa on other land / Food crops 

ea_cc_off_hfz <- landuses_EA_outUTMF_hfz %>%
  filter(LANDUSE == "COCOA") %>%
  group_by(LANDUSE, SCENARIO) %>% 
  mutate(EA_OUT_HA = ifelse(EA_OUT_HA < 0, 0, EA_OUT_HA)) %>%
  left_join(land_for_cc_expansion_hfz, by = c("HFZ", "SCENARIO")) %>%
  group_by(HFZ, SCENARIO) %>% 
  mutate(EA_CC_ON_OL = EA_OUT_HA * (OTHERLAND_HA / LAND_FOR_CC_EXPANSION)) %>%
  mutate(EA_CC_ON_FC = EA_OUT_HA * (CROP_HA / LAND_FOR_CC_EXPANSION)) %>%
  # Check caps = ok for all
  mutate(CHECK_CAP_OL = ifelse(EA_CC_ON_OL > OTHERLAND_HA, "STOP", "OK")) %>%
  mutate(CHECK_CAP_FC = ifelse(EA_CC_ON_FC > CROP_HA, "STOP", "OK"))

```


## SECTION 2.3.2.b.	Amount of food crop displacement and expansion

### a) Per region

```{r}
# Calculate total expansion + displacement of food crops

crops_dea_ha <- crops_EA %>% 
  left_join(ea_lui_off %>% 
              select(REGION, SCENARIO, EA_LUi_ON_FC), by = "REGION") %>% 
  left_join(ea_cc_off %>% 
              select(REGION, SCENARIO, EA_CC_ON_FC), by = c("REGION", "SCENARIO")) %>% 
  group_by(REGION, SCENARIO) %>% 
  mutate(DEA_FC_HA = EA_HA + EA_LUi_ON_FC + EA_CC_ON_FC)
```

```{r}
# Other deforestation & degradation per region to cap the DEA of food crops on forests 

other_defdeg_rg <- landuses_def_rg %>% 
              mutate(DRIVER = ifelse(grepl("COCOA", DRIVER), "COCOA_DD_HA", DRIVER)) %>%
              group_by(REGION, DRIVER) %>%
              mutate(AREA_HA = sum(AREA_HA)) %>%
              unique() %>%
              mutate(LANDUSE = case_when(
                grepl("COCOA", DRIVER) ~ "COCOA",
                grepl("MINING", DRIVER) ~ "MINING",
                grepl("TREEPLANTATION", DRIVER) ~ "TREEPLANTATION",
                grepl("SETTLEMENT", DRIVER) ~ "SETTLEMENT",
                grepl("OVERLAP", DRIVER) ~ "OVERLAP",
                grepl("OTHER", DRIVER) ~ "OTHER",
                TRUE ~ "NA"
              )) %>% 
  filter(LANDUSE == "OTHER") %>% 
  group_by(REGION) %>% 
  reframe(OTHER_DEFDEG_HA = sum(AREA_HA))
```


### b) On all HFZ

```{r}
# Calculate total expansion + displacement of food crops

crops_dea_ha_hfz <- crops_EA %>% 
  mutate(HFZ = "HFZ") %>% 
  group_by(HFZ) %>% 
  reframe(EA_HA = sum(EA_HA)) %>% 
  left_join(ea_lui_off_hfz %>% 
              select(HFZ, SCENARIO, EA_LUi_ON_FC), by = c("HFZ")) %>% 
  left_join(ea_cc_off_hfz %>% 
              select(HFZ, SCENARIO, EA_CC_ON_FC), by = c("HFZ", "SCENARIO")) %>% 
  mutate(DEA_FC_HA = EA_HA + EA_LUi_ON_FC + EA_CC_ON_FC)
```

```{r}
# Other deforestation & degradation per region to cap the DEA of food crops on forests 

other_defdeg_hfz <- landuses_def_rg %>% 
              mutate(DRIVER = ifelse(grepl("COCOA", DRIVER), "COCOA_DD_HA", DRIVER)) %>%
              group_by(REGION, DRIVER) %>%
              mutate(AREA_HA = sum(AREA_HA), HFZ = "HFZ") %>%
              unique() %>%
              mutate(LANDUSE = case_when(
                grepl("COCOA", DRIVER) ~ "COCOA",
                grepl("MINING", DRIVER) ~ "MINING",
                grepl("TREEPLANTATION", DRIVER) ~ "TREEPLANTATION",
                grepl("SETTLEMENT", DRIVER) ~ "SETTLEMENT",
                grepl("OVERLAP", DRIVER) ~ "OVERLAP",
                grepl("OTHER", DRIVER) ~ "OTHER",
                TRUE ~ "NA"
              )) %>% 
  filter(LANDUSE == "OTHER") %>% 
  group_by(HFZ) %>% 
  reframe(OTHER_DEFDEG_HA = sum(AREA_HA))
```


## SECTION 2.3.2.c.	Land balance: Distribution of food crop displacement and expansion

### a) Per region

```{r}
# Denominateur : calculate the total land for expansion of food crops, by subtracting all the land on which LUi and cocoa expanded 

converted_forest_rg <- landuses_def_rg %>% 
              mutate(DRIVER = ifelse(grepl("COCOA", DRIVER), "COCOA_DD_HA", DRIVER)) %>%
              group_by(REGION, DRIVER) %>%
              mutate(AREA_HA = sum(AREA_HA)) %>%
              unique() %>%
              mutate(LANDUSE = case_when(
                grepl("COCOA", DRIVER) ~ "COCOA",
                grepl("MINING", DRIVER) ~ "MINING",
                grepl("TREEPLANTATION", DRIVER) ~ "TREEPLANTATION",
                grepl("SETTLEMENT", DRIVER) ~ "SETTLEMENT",
                grepl("OVERLAP", DRIVER) ~ "OVERLAP",
                grepl("OTHER", DRIVER) ~ "OTHER",
                TRUE ~ "NA"
              )) %>% 
  filter(LANDUSE != "OTHER") %>% 
  group_by(REGION) %>% 
  reframe(CONVERTED_HA = sum(AREA_HA))

# Other land = Region area minus all land uses in 2000 (including water, regrowth and DTMF, UTMF and all LU)
# Left Forest = Only UTMF in 2000 (not DTMF) because "Converted forest" was only calculated on UTMF for other LU
# ... Thus makes more sense to only monitor FC def & deg on UTMF, even though we may slightly over estimate UTMF defdeg 
# ... Because we will direct FC DEA on UTMF while it may have replace DTMF in reality. 
# ... But as there is only 24Kha of DTMF in 2000 compared to 2.8 Mha of UTMF, we consider that it is too low to worry.
left_forest <- landuses_rg_2000 %>% 
  subset(LANDUSE == "UTMF") %>% 
  group_by(REGION, SCENARIO) %>% 
  reframe(FOREST_HA_2000 = sum(LANDUSE_HA_2000)) %>% 
  left_join(converted_forest_rg, by = "REGION") %>%
  group_by(REGION, SCENARIO) %>%
  mutate(FOREST_LEFT_HA = FOREST_HA_2000 - CONVERTED_HA)

land_for_fc_expansion <- ea_lui_off %>% 
  select(-CHECK_CAP_OL, -CHECK_CAP_FC, - CHECK_CAP_CC) %>% 
  left_join(ea_cc_off %>% 
              select(REGION, EA_CC_ON_OL, EA_CC_ON_FC, SCENARIO), by = c("REGION", "SCENARIO")) %>% 
  group_by(REGION, SCENARIO) %>% 
  mutate(OTHERLAND_LEFT_HA = OTHERLAND_HA - EA_LUi_ON_OL - EA_CC_ON_OL, 
         CONVERTED_CC_LEFT_HA = -CONVERTED_CC_HA - EA_LUi_ON_CC) %>% 
  left_join(left_forest, by = c("REGION", "SCENARIO")) %>%  
  group_by(REGION, SCENARIO) %>% 
  mutate(LAND_FOR_FC_EXPANSION = OTHERLAND_LEFT_HA + CONVERTED_CC_LEFT_HA + FOREST_LEFT_HA) 

```

```{r}
# Calculate the expansion & displacement of food crops (DEA_FC) on forest / other land / converted cocoa

crops_dea_notcapped <- crops_dea_ha %>% 
  select(-EA_LUi_ON_FC, -EA_CC_ON_FC) %>% 
  left_join(land_for_fc_expansion, by = c("REGION", "SCENARIO")) %>% 
  group_by(REGION, SCENARIO) %>% 
  mutate(DEA_FC_ON_F = DEA_FC_HA * (FOREST_LEFT_HA / LAND_FOR_FC_EXPANSION), # Includes direct deforestation due to FC and "indirect" due to LU&CC displacing FC
         DEA_FC_ON_OL = DEA_FC_HA * (OTHERLAND_LEFT_HA / LAND_FOR_FC_EXPANSION),
         DEA_FC_ON_CC = DEA_FC_HA * (CONVERTED_CC_LEFT_HA / LAND_FOR_FC_EXPANSION)) %>% 
  ungroup() %>% 
  # Check caps = ok for all EXCEPT for other def & def in Western and central regions
  left_join(other_defdeg_rg, by = "REGION") %>% 
  group_by(REGION, SCENARIO) %>% 
  mutate(CHECK_CAP_OTHERDEFDEG = ifelse(DEA_FC_ON_F > OTHER_DEFDEG_HA, "STOP", "OK"),
         CHECK_CAP_OL = ifelse(DEA_FC_ON_OL > OTHERLAND_LEFT_HA, "STOP", "OK"),
         CHECK_CAP_CC = ifelse(DEA_FC_ON_CC > CONVERTED_CC_LEFT_HA , "STOP", "OK"))  %>% 
  ungroup()
# crops_dea_notcapped
```

```{r}
# Capped expansion and displacement of food crops 

crops_dea <- crops_dea_notcapped %>% 
  group_by(REGION, SCENARIO) %>% 
  mutate(DEA_FC_ON_F_CAPPED = ifelse(DEA_FC_ON_F > OTHER_DEFDEG_HA, OTHER_DEFDEG_HA, DEA_FC_ON_F)) %>% 
  mutate(DIFF_FC = DEA_FC_ON_F - DEA_FC_ON_F_CAPPED) %>% 
  mutate(DEA_FC_ON_OL_ADJUSTED = DEA_FC_ON_OL + DIFF_FC) %>% 
  mutate(CHECK_NEW_CAP_OTHERDEFDEG = ifelse(DEA_FC_ON_F_CAPPED > OTHER_DEFDEG_HA, "STOP", "OK")) %>% 
  mutate(CHECK_NEW_CAP_OL = ifelse(DEA_FC_ON_OL_ADJUSTED > OTHERLAND_LEFT_HA, "STOP", "OK")) %>% 
  select(-DEA_FC_ON_F, -DEA_FC_ON_OL_ADJUSTED)

```
 

```{r}
# Table

landuse_crops_EA_tot <-  crops_dea %>%  
  group_by(SCENARIO) %>% 
  summarise_if(is.numeric, sum, na.rm=TRUE)
# landuse_crops_EA_tot
```


### b) On all HFZ

```{r}
# Denominateur : calculate the total land for expansion of food crops, by subtracting all the land on which LUi and cocoa expanded 

converted_forest_hfz <- landuses_def_rg %>% 
              mutate(DRIVER = ifelse(grepl("COCOA", DRIVER), "COCOA_DD_HA", DRIVER)) %>%
              group_by(REGION, DRIVER) %>%
              mutate(HFZ = "HFZ", AREA_HA = sum(AREA_HA)) %>%
              unique() %>%
              mutate(LANDUSE = case_when(
                grepl("COCOA", DRIVER) ~ "COCOA",
                grepl("MINING", DRIVER) ~ "MINING",
                grepl("TREEPLANTATION", DRIVER) ~ "TREEPLANTATION",
                grepl("SETTLEMENT", DRIVER) ~ "SETTLEMENT",
                grepl("OVERLAP", DRIVER) ~ "OVERLAP",
                grepl("OTHER", DRIVER) ~ "OTHER",
                TRUE ~ "NA"
              )) %>% 
  filter(LANDUSE != "OTHER") %>% 
  group_by(HFZ) %>% 
  reframe(CONVERTED_HA = sum(AREA_HA))

left_forest_hfz <- landuses_rg_2000 %>% 
  subset(LANDUSE == "UTMF" | LANDUSE == "DTMF") %>% 
  mutate(HFZ = "HFZ") %>% 
  group_by(HFZ, SCENARIO) %>% 
  reframe(FOREST_HA_2000 = sum(LANDUSE_HA_2000)) %>% 
  left_join(converted_forest_hfz, by = "HFZ") %>% 
  group_by(HFZ, SCENARIO) %>%
  mutate(FOREST_LEFT_HA = FOREST_HA_2000 - CONVERTED_HA)

land_for_fc_expansion_hfz <- ea_lui_off_hfz %>% 
  select(-CHECK_CAP_OL, -CHECK_CAP_FC, - CHECK_CAP_CC) %>% 
  left_join(ea_cc_off_hfz %>% 
              select(HFZ, SCENARIO, EA_CC_ON_OL, EA_CC_ON_FC), by = c("HFZ", "SCENARIO")) %>% 
  group_by(HFZ, SCENARIO) %>%
  mutate(OTHERLAND_LEFT_HA = OTHERLAND_HA - EA_LUi_ON_OL - EA_CC_ON_OL) %>% 
  mutate(CONVERTED_CC_LEFT_HA = -CONVERTED_CC_HA - EA_LUi_ON_CC) %>% 
  left_join(left_forest_hfz, by = c("HFZ", "SCENARIO")) %>%  
  group_by(HFZ, SCENARIO) %>%
  mutate(LAND_FOR_FC_EXPANSION = OTHERLAND_LEFT_HA + CONVERTED_CC_LEFT_HA + FOREST_LEFT_HA) 

```

```{r}
# Calculate the expansion & displacement of food crops (DEA_FC) on forest / other land / converted cocoa

crops_dea_notcapped_hfz <- crops_dea_ha_hfz %>% 
  select(-EA_LUi_ON_FC, -EA_CC_ON_FC) %>% 
  left_join(land_for_fc_expansion_hfz, by = c("HFZ", "SCENARIO")) %>% 
  group_by(HFZ, SCENARIO) %>%
  mutate(DEA_FC_ON_F = DEA_FC_HA * (FOREST_LEFT_HA / LAND_FOR_FC_EXPANSION)) %>% 
  mutate(DEA_FC_ON_OL = DEA_FC_HA * (OTHERLAND_LEFT_HA / LAND_FOR_FC_EXPANSION)) %>% 
  mutate(DEA_FC_ON_CC = DEA_FC_HA * (CONVERTED_CC_LEFT_HA / LAND_FOR_FC_EXPANSION)) %>% 
  # Check caps = ok 
  left_join(other_defdeg_hfz, by = "HFZ") %>% 
  mutate(CHECK_CAP_OTHERDEFDEG = ifelse(DEA_FC_ON_F > OTHER_DEFDEG_HA, "STOP", "OK")) %>% 
  mutate(CHECK_CAP_OL = ifelse(DEA_FC_ON_OL > OTHERLAND_LEFT_HA, "STOP", "OK")) %>% 
  mutate(CHECK_CAP_CC = ifelse(DEA_FC_ON_CC > CONVERTED_CC_LEFT_HA , "STOP", "OK"))  
```

```{r}
# Capped expansion and displacement of food crops 

crops_dea_hfz <- crops_dea_notcapped_hfz %>% 
  group_by(SCENARIO) %>%
  mutate(DEA_FC_ON_F_CAPPED = ifelse(DEA_FC_ON_F > OTHER_DEFDEG_HA, OTHER_DEFDEG_HA, DEA_FC_ON_F)) %>% 
  mutate(DIFF_FC = DEA_FC_ON_F - DEA_FC_ON_F_CAPPED) %>% 
  mutate(DEA_FC_ON_OL_ADJUSTED = DEA_FC_ON_OL + DIFF_FC) %>% 
  mutate(CHECK_NEW_CAP_OTHERDEFDEG = ifelse(DEA_FC_ON_F_CAPPED > OTHER_DEFDEG_HA, "STOP", "OK")) %>% 
  mutate(CHECK_NEW_CAP_OL = ifelse(DEA_FC_ON_OL_ADJUSTED > OTHERLAND_LEFT_HA, "STOP", "OK")) %>% 
  select(-DEA_FC_ON_F, -DEA_FC_ON_OL_ADJUSTED)

```

# SECTION 2.3.3. Attributing the deforestation and degradation of land used for food crops in 2019 to the expansion of other land uses  

### a) Per region

```{r}
# Denominateur : calculate the total EA of LUi off F + EA of CC off F + EA FC
# And distribute the DEF of FC to the three LUI/CC/FC

shared_fc_def_rg <- crops_dea %>% 
  select(REGION, SCENARIO, DEA_FC_ON_F_CAPPED, DEA_FC_HA, EA_LUi_ON_FC, EA_CC_ON_FC, EA_FC = EA_HA) %>% 
  group_by(REGION, SCENARIO) %>%
  mutate(SHARE_LUi = DEA_FC_ON_F_CAPPED * (EA_LUi_ON_FC / DEA_FC_HA)) %>% 
  mutate(SHARE_COCOA = DEA_FC_ON_F_CAPPED * (EA_CC_ON_FC / DEA_FC_HA)) %>% 
  mutate(SHARE_FOODCROP = DEA_FC_ON_F_CAPPED * (EA_FC / DEA_FC_HA)) %>% 
  left_join(ea_lui_off %>% select(REGION, SCENARIO, EA_LUi_OFF_HA), by = c("REGION", "SCENARIO")) %>% 
  left_join(ea_cc_off %>% select(REGION, SCENARIO, EA_OUT_HA), by = c("REGION", "SCENARIO")) %>% 
# CHECK if food crop deforestation allocated to a LUi or CC is higher than its OUTSIDE FOREST expansion: Ok all.
  mutate(CHECK_CAP_LUi = ifelse(SHARE_LUi > EA_LUi_OFF_HA, "STOP", "OK")) %>%
  mutate(CHECK_CAP_CC = ifelse(SHARE_COCOA > EA_OUT_HA, "STOP", "OK")) %>%
  mutate(CHECK_CAP_FC = ifelse(SHARE_FOODCROP > EA_FC, "STOP", "OK")) 

# shared_fc_def_rg 

shared_fc_def_rg_long <- shared_fc_def_rg %>% 
  select(REGION, SCENARIO, DEA_FC_ON_F_CAPPED, SHARE_LUi, SHARE_COCOA, SHARE_FOODCROP) %>% 
  pivot_longer(cols = c(4:6), 
               names_to = "METRIC", 
               values_to = "AREA_HA") %>% 
  mutate(PERC = AREA_HA/DEA_FC_ON_F_CAPPED*100)
```

```{r}
# Distribute the DEF of FC to the each LU separately MINING/SETTLEMENT/TREEPLANTATIONS/CC/FC

shared_fc_def_per_lui_rg <- crops_dea %>% 
  select(REGION, SCENARIO, DEA_FC_ON_F_CAPPED, DEA_FC_HA, EA_LUi_ON_FC, EA_CC_ON_FC, EA_FC = EA_HA) %>% 
  left_join(ea_per_lui_off %>% 
              select(REGION, SCENARIO, LANDUSE, EA_LUi_ON_FC) %>% 
              pivot_wider(names_from = LANDUSE, 
                          values_from = EA_LUi_ON_FC), by = c("REGION", "SCENARIO")) %>% 
  rename(EA_MINING_ON_FC = MINING, EA_TREEPLANTATION_ON_FC = TREEPLANTATION, EA_SETTLEMENTS_ON_FC = SETTLEMENT) %>% 
  group_by(REGION, SCENARIO) %>% 
  mutate(SHARE_MINING = DEA_FC_ON_F_CAPPED * (EA_MINING_ON_FC / DEA_FC_HA)) %>% 
  mutate(SHARE_TREEPLANTATION = DEA_FC_ON_F_CAPPED * (EA_TREEPLANTATION_ON_FC / DEA_FC_HA)) %>% 
  mutate(SHARE_SETTLEMENTS = DEA_FC_ON_F_CAPPED * (EA_SETTLEMENTS_ON_FC / DEA_FC_HA)) %>% 
  mutate(SHARE_COCOA = DEA_FC_ON_F_CAPPED * (EA_CC_ON_FC / DEA_FC_HA)) %>% 
  mutate(SHARE_FOODCROP = DEA_FC_ON_F_CAPPED * (EA_FC / DEA_FC_HA)) 


# # CHECK if food crop deforestation allocated to each LUi is higher than their OUTSIDE FOREST expansion: Ok all.
# shared_fc_def_per_lui_rg %>% 
#   select(REGION, SCENARIO, SHARE_MINING, SHARE_TREEPLANTATION, SHARE_SETTLEMENTS) %>% 
#   pivot_longer(cols = c(3:5),
#                names_to = "LANDUSE", 
#                values_to = "AREA_HA") %>% 
#      mutate(LANDUSE = case_when(
#                 grepl("MINING", LANDUSE) ~ "MINING",
#                 grepl("TREEPLANTATION", LANDUSE) ~ "TREEPLANTATION",
#                 grepl("SETTLEMENT", LANDUSE) ~ "SETTLEMENT",
#                 TRUE ~ "NA"
#               )) %>% 
#   left_join(ea_per_lui_off %>% select(REGION, SCENARIO, LANDUSE, EA_LUi_OFF_HA), by = c("REGION", "LANDUSE", "SCENARIO")) %>% 
#   mutate(CHECK_CAP_LUi = ifelse(AREA_HA > EA_LUi_OFF_HA, "STOP", "OK"))


shared_fc_def_per_lui_rg_long <- shared_fc_def_per_lui_rg %>% 
  select(REGION, SCENARIO, DEA_FC_ON_F_CAPPED, SHARE_COCOA, SHARE_MINING, SHARE_TREEPLANTATION, SHARE_SETTLEMENTS, SHARE_FOODCROP) %>% 
  pivot_longer(cols = c(4:8), 
               names_to = "METRIC", 
               values_to = "AREA_HA") %>% 
  mutate(PERC = AREA_HA/DEA_FC_ON_F_CAPPED*100)
```


### b) On all HFZ

Doing the analysis directly over all merged cocoa regions rather than separately and summed up as above

```{r}
# Denominateur : calculate the total EA of LUi off F + EA of CC off F + EA FC
# And distribute the DEF of FC to the three LUI/CC/FC

shared_fc_def_hfz <- crops_dea_hfz %>% 
  select(HFZ, SCENARIO, DEA_FC_ON_F_CAPPED, DEA_FC_HA, EA_LUi_ON_FC, EA_CC_ON_FC, EA_FC = EA_HA) %>% 
  group_by(HFZ, SCENARIO) %>% 
  mutate(SHARE_LUi = DEA_FC_ON_F_CAPPED * (EA_LUi_ON_FC / DEA_FC_HA)) %>% 
  mutate(SHARE_COCOA = DEA_FC_ON_F_CAPPED * (EA_CC_ON_FC / DEA_FC_HA)) %>% 
  mutate(SHARE_FOODCROP = DEA_FC_ON_F_CAPPED * (EA_FC / DEA_FC_HA)) %>% 
  left_join(ea_lui_off_hfz %>% select(HFZ, SCENARIO, EA_LUi_OFF_HA), by = c("HFZ", "SCENARIO")) %>% 
  left_join(ea_cc_off_hfz %>% select(HFZ, SCENARIO, EA_OUT_HA), by = c("HFZ", "SCENARIO")) %>% 
# CHECK if food crop deforestation allocated to a LUi or CC is higher than its OUTSIDE FOREST expansion: Ok all.
  mutate(CHECK_CAP_LUi = ifelse(SHARE_LUi > EA_LUi_OFF_HA, "STOP", "OK")) %>%
  mutate(CHECK_CAP_CC = ifelse(SHARE_COCOA > EA_OUT_HA, "STOP", "OK")) %>%
  mutate(CHECK_CAP_FC = ifelse(SHARE_FOODCROP > EA_FC, "STOP", "OK")) 

# shared_fc_def_hfz 

shared_fc_def_hfz_long <- shared_fc_def_hfz %>% 
  select(HFZ, SCENARIO, DEA_FC_ON_F_CAPPED, SHARE_LUi, SHARE_COCOA, SHARE_FOODCROP) %>% 
  pivot_longer(cols = c(4:6), 
               names_to = "METRIC", 
               values_to = "AREA_HA") %>% 
  mutate(PERC = AREA_HA/DEA_FC_ON_F_CAPPED*100)
```

```{r}
# Distribute the DEF of FC to the each LU separately MINING/SETTLEMENT/TREEPLANTATIONS/CC/FC

shared_fc_def_per_lui_hfz <- crops_dea_hfz %>% 
  select(HFZ, SCENARIO, DEA_FC_ON_F_CAPPED, DEA_FC_HA, EA_LUi_ON_FC, EA_CC_ON_FC, EA_FC = EA_HA) %>% 
  left_join(ea_per_lui_off_hfz %>% 
              select(HFZ, SCENARIO, LANDUSE, EA_LUi_ON_FC) %>% 
              pivot_wider(names_from = LANDUSE, 
                          values_from = EA_LUi_ON_FC), by = c("HFZ", "SCENARIO")) %>% 
  rename(EA_MINING_ON_FC = MINING, EA_TREEPLANTATION_ON_FC = TREEPLANTATION, EA_SETTLEMENTS_ON_FC = SETTLEMENT) %>% 
  group_by(HFZ, SCENARIO) %>% 
  mutate(SHARE_MINING = DEA_FC_ON_F_CAPPED * (EA_MINING_ON_FC / DEA_FC_HA)) %>% 
  mutate(SHARE_TREEPLANTATION = DEA_FC_ON_F_CAPPED * (EA_TREEPLANTATION_ON_FC / DEA_FC_HA)) %>% 
  mutate(SHARE_SETTLEMENTS = DEA_FC_ON_F_CAPPED * (EA_SETTLEMENTS_ON_FC / DEA_FC_HA)) %>% 
  mutate(SHARE_COCOA = DEA_FC_ON_F_CAPPED * (EA_CC_ON_FC / DEA_FC_HA)) %>% 
  mutate(SHARE_FOODCROP = DEA_FC_ON_F_CAPPED * (EA_FC / DEA_FC_HA)) 


# CHECK if food crop deforestation allocated to each LUi is higher than their OUTSIDE FOREST expansion: Ok all.
# shared_fc_def_per_lui_hfz %>% 
#   select(HFZ, SCENARIO, SHARE_MINING, SHARE_TREEPLANTATION, SHARE_SETTLEMENTS) %>% 
#   pivot_longer(cols = c(3:5),
#                names_to = "LANDUSE", 
#                values_to = "AREA_HA") %>% 
#      mutate(LANDUSE = case_when(
#                 grepl("MINING", LANDUSE) ~ "MINING",
#                 grepl("TREEPLANTATION", LANDUSE) ~ "TREEPLANTATION",
#                 grepl("SETTLEMENT", LANDUSE) ~ "SETTLEMENT",
#                 TRUE ~ "NA"
#               )) %>% 
#   left_join(ea_per_lui_off_hfz %>% select(HFZ, SCENARIO, LANDUSE, EA_LUi_OFF_HA), by = c("HFZ", "LANDUSE", "SCENARIO")) %>% 
#   mutate(CHECK_CAP_LUi = ifelse(AREA_HA > EA_LUi_OFF_HA, "STOP", "OK"))


shared_fc_def_per_lui_hfz_long <- shared_fc_def_per_lui_hfz %>% 
  select(HFZ, SCENARIO, DEA_FC_ON_F_CAPPED, SHARE_COCOA, SHARE_MINING, SHARE_TREEPLANTATION, SHARE_SETTLEMENTS, SHARE_FOODCROP) %>% 
  pivot_longer(cols = c(4:8), 
               names_to = "METRIC", 
               values_to = "AREA_HA") %>% 
  mutate(PERC = AREA_HA/DEA_FC_ON_F_CAPPED*100)
```


# SECTION 3.1. DIRECT AND INDIRECT DEFORESTATION AND DEGRADATION

### a) Per region


```{r}
tempa <- landuses_def_rg %>% 
              mutate(DRIVER = ifelse(grepl("OTHER", DRIVER), "OTHER_DFDG_HA", DRIVER)) %>%
              group_by(REGION, DRIVER) %>%
              mutate(AREA_HA = sum(AREA_HA)) %>%
              unique() %>%
              mutate(LANDUSE = case_when(
                grepl("COCOA", DRIVER) ~ "COCOA",
                grepl("MINING", DRIVER) ~ "MINING",
                grepl("TREEPLANTATION", DRIVER) ~ "TREE PLANTATION",
                grepl("SETTLEMENT", DRIVER) ~ "SETTLEMENT",
                grepl("OVERLAP", DRIVER) ~ "OVERLAP",
                grepl("OTHER", DRIVER) ~ "OTHER",
                TRUE ~ "NA")) %>% 
              mutate(TYPE = case_when(
                grepl("DEG", DRIVER) ~ "DIRECT DEGRADATION",
                grepl("DEF", DRIVER) ~ "DIRECT DEFORESTATION",
                grepl("DD", DRIVER) ~ "DIRECT DEFORESTATION",
                grepl("DFDG", DRIVER) ~ "DEFORESTATION & DEGRADATION",
                TRUE ~ "NA")) %>% 
  mutate(SCENARIO = "MID_EXP") %>% 
    select(REGION, LANDUSE, TYPE, SCENARIO, AREA_HA)

tempb <- shared_fc_def_per_lui_rg_long %>%
    mutate(LANDUSE = case_when(
                grepl("COCOA", METRIC) ~ "COCOA",
                grepl("MINING", METRIC) ~ "MINING",
                grepl("TREEPLANTATION", METRIC) ~ "TREE PLANTATION",
                grepl("SETTLEMENT", METRIC) ~ "SETTLEMENT",
                grepl("FOODCROP", METRIC) ~ "FOOD CROP",
                TRUE ~ "NA")) %>% 
              mutate(TYPE = "INDIRECT DEFORESTATION & DEGRADATION") %>% 
  mutate(TYPE = case_when(
                grepl("FOOD CROP", LANDUSE) ~ "DEFORESTATION & DEGRADATION",
                TRUE ~ TYPE)) %>% 
  select(REGION, LANDUSE, TYPE, SCENARIO, AREA_HA)

tempc <- tempa %>% rbind(tempb)

# Substract all the indirect deforestation & foodcrop direct def&deg to "other DD"
indirectDD_rg <- tempc %>% 
  filter(TYPE == "INDIRECT DEFORESTATION & DEGRADATION") %>% 
  group_by(REGION, SCENARIO, TYPE) %>% 
  mutate(AREA_HA = sum(AREA_HA)) %>% 
  mutate(LANDUSE = "ALL_LUI_CC") %>% 
  select(-DRIVER) %>% 
  unique()

fcDD_rg <- tempc %>% 
  filter(LANDUSE == "FOOD CROP") %>% 
  select(-DRIVER)

fc_direct_ind_rg <- indirectDD_rg %>% 
  rbind(fcDD_rg) %>% 
  group_by(REGION, SCENARIO) %>% 
  reframe(FC_totDD_AREA_HA = sum(AREA_HA)) 

timberDD_rg <- tempc %>% 
  filter(LANDUSE == "OTHER") %>% 
  select(-SCENARIO) %>% 
  left_join(fc_direct_ind_rg, by = c("REGION")) %>% 
  group_by(REGION, SCENARIO) %>% 
  mutate(TIMBER_HA = AREA_HA - FC_totDD_AREA_HA) %>% 
  reframe(REGION, LANDUSE = "LOGGING\nAND OTHER", TYPE, SCENARIO, AREA_HA = TIMBER_HA)

tempd <- tempc %>% 
  filter(LANDUSE != "OTHER") %>% 
  rbind(timberDD_rg)
```

#### fig.4B: 2000-2019 direct and indirect deforestation and degradation per region and per land use in hectares

```{r, fig.width = 13, fig.height = 7}

ggplot(data = tempd %>% pivot_wider(names_from = SCENARIO, values_from = AREA_HA) %>% 
         mutate(LANDUSE_TYPE = paste0(LANDUSE, " ", TYPE)), 
       aes(x = LANDUSE, y = MID_EXP, fill = LANDUSE_TYPE)) + 
    geom_bar(stat="identity", position = "stack") +
  geom_errorbar(aes(ymin = MAX_EXP, ymax = MIN_EXP), width = 0.2) +
  facet_wrap(~REGION) +
    coord_flip() +
    theme(axis.title.y = element_text(angle = 90, size=14),
        axis.text.y = element_text(size = 10),
        axis.title.x = element_text(size = 12),
        axis.text.x = element_text(size = 10), 
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 10),
        strip.background = element_rect(colour = "grey50"),
        legend.position = "right", 
        legend.text = element_text(size = 8), 
        legend.title = element_text(size = 10),
        panel.spacing = unit(1.5, "lines"), 
        panel.background = element_rect(fill = "white", colour = "grey50"), 
        panel.grid.major.y = element_blank(), 
        panel.grid.major.x = element_line(colour = "grey90")) +
   scale_fill_manual(values = def_scale) +
  labs(y = "Area deforested and/or degraded per land use (ha) 2000-2019", x = "", fill = "Direct and indirect deforestation and degradation\nper land use (ha) 2000-2019") +
  scale_y_continuous(label = unit_format(unit = "K", scale = 1e-3, sep = ""))

```

#### fig.4B: aggregated at the HFZ level

```{r,fig.height=5, fig.width=7}

ggplot(data = tempd %>% 
         group_by(LANDUSE, TYPE, SCENARIO) %>% 
         reframe(AREA_HA = sum(AREA_HA)) %>% 
         pivot_wider(names_from = SCENARIO, values_from = AREA_HA) %>% 
         mutate(LANDUSE_TYPE = paste0(LANDUSE, " ", TYPE)) %>% 
         mutate(TITLE = "TOTAL AGGREGATED\nAT HFZ LEVEL"), 
       aes(x = LANDUSE, y = MID_EXP, fill = LANDUSE_TYPE)) +
    geom_bar(stat="identity", position = "stack") +
  geom_errorbar(aes(ymin = MAX_EXP, ymax = MIN_EXP), width = 0.2) +
  facet_grid(. ~ TITLE) +
   coord_flip() +
      theme(axis.title.y = element_text(angle = 90, size=16),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_blank(),
        axis.text.x = element_text(size = 16), 
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 16),
        strip.background = element_rect(colour = "grey50"),
        legend.position = "right", 
        legend.text = element_blank(), 
        legend.title = element_blank(),
        panel.spacing = unit(1.5, "lines"), 
        panel.background = element_rect(fill = "white", colour = "grey50"), 
        panel.grid.major.y = element_blank(), 
        panel.grid.major.x = element_line(colour = "grey90")) +
   scale_fill_manual(values = def_scale) +
  labs(y = "Area (ha)", x = "") +
  scale_y_continuous(label = unit_format(unit = "K", scale = 1e-3, sep = ""))

```


#### fig.4C + fig.S4: 2000-2019 direct and indirect deforestation and degradation per region and per land use in percentage

```{r map-forestcover-3combi, fig.width = 7, fig.height = 5}

# To get figures S4, change y = MID_EXP to MIN_EXP or MAX_EXP in the ggplot

# Adding direct def & deg to Max and Min scenarios
def_rg_per_lu_dir_indir <- tempd %>% 
  rbind(tempa %>% mutate(SCENARIO = "MIN_EXP") %>% filter(LANDUSE != "OTHER")) %>%  
  rbind(tempa %>% mutate(SCENARIO = "MAX_EXP")  %>% filter(LANDUSE != "OTHER"))

# Total deforestation & degradation for each scenario 
## for HFZ
tot_def_hfz_agg <- def_rg_per_lu_dir_indir %>% 
  group_by(SCENARIO) %>% 
  reframe(TOT_AREA_HA = sum(AREA_HA))

## per region
tot_def_rg_agg <- def_rg_per_lu_dir_indir %>%
  group_by(REGION, SCENARIO) %>% 
  reframe(TOT_AREA_HA = sum(AREA_HA))

# Calculate percentage of each LU in the HFZ
sum_tot_def_hfz <- def_rg_per_lu_dir_indir %>% 
  left_join(tot_def_hfz_agg, by = "SCENARIO") %>% 
  group_by(LANDUSE, TYPE, SCENARIO) %>% 
  reframe(AREA_HA = sum(AREA_HA), TOT_AREA_HA) %>% 
  mutate(PERC = AREA_HA/TOT_AREA_HA *100, REGION = "HFZ") %>% 
  unique()

sum_tot_def_rg_hfz <- def_rg_per_lu_dir_indir %>% 
  left_join(tot_def_rg_agg, by = c("REGION", "SCENARIO")) %>% 
  group_by(REGION, LANDUSE, TYPE, SCENARIO) %>% 
  mutate(PERC = AREA_HA/TOT_AREA_HA*100) %>% 
  select(-DRIVER) %>% 
  rbind(sum_tot_def_hfz)


# Plot percentage of total direct & indirect def & deg per LU 
sum_tot_def_rg_hfz %>% 
         pivot_wider(names_from = SCENARIO, values_from = PERC) %>% 
         mutate(LANDUSE_TYPE = paste0(LANDUSE, " ", TYPE)) %>% 
         mutate(REGION = factor(REGION, 
                                levels= c("ASHANTI", "BRONG AHAFO", "CENTRAL", "EASTERN", "VOLTA", "WESTERN", "HFZ"))) %>% 
         arrange(REGION) %>% 
  #filter(REGION == "HFZ") %>% 
  mutate(LANDUSE_TYPE = factor(LANDUSE_TYPE, 
                                      levels = c(
                                        "TREE PLANTATION DIRECT DEFORESTATION",
                                        "SETTLEMENT DIRECT DEFORESTATION", 
                                        "MINING DIRECT DEFORESTATION",
                                        "COCOA DIRECT DEGRADATION", 
                                        "COCOA DIRECT DEFORESTATION",
                                        "LOGGING\nAND OTHER DEFORESTATION & DEGRADATION",
                                        "TREE PLANTATION INDIRECT DEFORESTATION & DEGRADATION",
                                        "SETTLEMENT INDIRECT DEFORESTATION & DEGRADATION",
                                        "MINING INDIRECT DEFORESTATION & DEGRADATION", 
                                        "COCOA INDIRECT DEFORESTATION & DEGRADATION", 
                                        "FOOD CROP DEFORESTATION & DEGRADATION"))) %>% 
        arrange(LANDUSE_TYPE) %>% 

ggplot(aes(x = REGION, y = MID_EXP, fill = LANDUSE_TYPE)) + 
    geom_bar(stat = "identity", position = "stack") + 
    coord_flip() +
      theme(axis.title.y = element_text(angle = 90, size=16),
        axis.text.y = element_text(size = 17),
        axis.title.x = element_text(size = 17),
        axis.text.x = element_text(size = 17), 
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 16),
        strip.background = element_rect(colour = "grey50"),
        legend.position = "none", 
        # legend.text = element_blank(), 
        # legend.title =  element_blank(),
        panel.spacing = unit(1.5, "lines"), 
        panel.background = element_rect(fill = "white", colour = "grey50"), 
        panel.grid.major.y = element_blank(), 
        panel.grid.major.x = element_line(colour = "grey90")) +

     scale_fill_manual(values = def_scale_levels) + #def_scale_levels
  labs(y = "Area deforested and/or degraded\nper land use (%)", x = "")

```


#### Table S8: 2000-2019 direct and indirect deforestation and forest degradation per region and per land use

```{r}

tot_def_rg_wide_perc <- sum_tot_def_rg_hfz %>% 
  select(-AREA_HA, -TOT_AREA_HA) %>% 
  mutate(PERC = round(PERC, 2)) %>% 
         pivot_wider(names_from = SCENARIO, values_from = PERC) %>% 
  rename(PERC_MID_EXP = MID_EXP, PERC_MAX_EXP = MAX_EXP, PERC_MIN_EXP = MIN_EXP) %>% 
  left_join(sum_tot_def_rg_hfz %>% 
              mutate(AREA_HA = round(AREA_HA, 2)) %>% 
  select(-PERC, -TOT_AREA_HA) %>% 
         pivot_wider(names_from = SCENARIO, values_from = AREA_HA) %>% 
  reframe(AREA_MID_EXP = MID_EXP, AREA_MAX_EXP = MAX_EXP, AREA_MIN_EXP = MIN_EXP), by = c("REGION", "LANDUSE", "TYPE"))

tot_def_rg_wide_perc  %>% 
  #mutate(AREA_HA = format(AREA_HA, digits = 2, big.mark = ",", decimal.mark = ".")) %>% 
  #mutate(PERC = format(PERC, digits = 2, big.mark = ",", decimal.mark = ".")) %>% 
  #dcast(., COCOA_MAP ~ FOREST_MAP) %>%
  kbl(caption = "Total direct & indirect deforestation (2000 - 2019) in ha", align = "c") %>% #col.names = c("Cocoa maps", names(.)[-1])
  kable_styling(full_width = F) %>% 
  #kable_paper() %>% 
  column_spec(1, bold = T, border_right = T) #add_header_above(c(" " = 1, "Forest maps" = 3))  
```

```{r}
#Table per landuse (combined direct & indirect)
tot_def_rg_perc_summed <- sum_tot_def_rg_hfz %>% 
  select(-AREA_HA, -TOT_AREA_HA) %>% 
  mutate(PERC = round(PERC, 2)) %>% 
  group_by(LANDUSE, REGION, SCENARIO) %>% 
  reframe(PERC = sum(PERC)) %>% 
  left_join(sum_tot_def_rg_hfz %>% 
              mutate(AREA_HA = round(AREA_HA, 2)) %>% 
  select(-PERC, -TOT_AREA_HA) %>% 
      group_by(LANDUSE, REGION, SCENARIO) %>% 
  reframe(AREA_HA = sum(AREA_HA)), by = c("REGION", "LANDUSE", "SCENARIO"))

tot_def_rg_perc_summed_wide <- tot_def_rg_perc_summed %>% 
  pivot_wider(names_from = 'SCENARIO', 
              values_from= c("PERC", "AREA_HA"))

```

```{r}
#Table summing up the total def&deg of food crops (direct & indirect)
sum_tot_def_rg_hfz %>% 
  filter(TYPE == "DEFORESTATION & DEGRADATION" | TYPE == "INDIRECT DEFORESTATION & DEGRADATION") %>% 
  filter(LANDUSE != "LOGGING\nAND OTHER") %>% 
  subset(REGION != "HFZ") %>% 
  group_by(REGION, SCENARIO) %>% 
  reframe(PERC = sum(PERC), AREA_HA = sum(AREA_HA))
  
```

```{r}
#Table summing up the indirect FC def and deg
sum_tot_def_rg_hfz %>% 
  filter(TYPE == "INDIRECT DEFORESTATION & DEGRADATION") %>% 
  filter(LANDUSE != "LOGGING\nAND OTHER") %>% 
  group_by(REGION, SCENARIO) %>% 
  reframe(PERC = sum(PERC), AREA_HA = sum(AREA_HA))
```



### b) On all HFZ

#### fig.S6A: 2000-2019 direct and indirect deforestation and degradation and per land use in hectares 

```{r,fig.height=5, fig.width=7}
tempa_hfz <- landuses_def_rg %>% 
              mutate(DRIVER = ifelse(grepl("OTHER", DRIVER), "OTHER_DFDG_HA", DRIVER)) %>%
              group_by(DRIVER) %>%
              mutate(AREA_HA = sum(AREA_HA)) %>%
              unique() %>%
              mutate(LANDUSE = case_when(
                grepl("COCOA", DRIVER) ~ "COCOA",
                grepl("MINING", DRIVER) ~ "MINING",
                grepl("TREEPLANTATION", DRIVER) ~ "TREE PLANTATION",
                grepl("SETTLEMENT", DRIVER) ~ "SETTLEMENT",
                grepl("OVERLAP", DRIVER) ~ "OVERLAP",
                grepl("OTHER", DRIVER) ~ "OTHER",
                TRUE ~ "NA")) %>% 
              mutate(TYPE = case_when(
                grepl("DEG", DRIVER) ~ "DIRECT DEGRADATION",
                grepl("DEF", DRIVER) ~ "DIRECT DEFORESTATION",
                grepl("DD", DRIVER) ~ "DIRECT DEFORESTATION",
                grepl("DFDG", DRIVER) ~ "DEFORESTATION & DEGRADATION",
                TRUE ~ "NA")) %>% 
  mutate(SCENARIO = "MID_EXP") %>% 
    select(LANDUSE, TYPE, SCENARIO, AREA_HA) %>% 
  unique()

tempb_hfz <- shared_fc_def_per_lui_hfz_long %>%
    mutate(LANDUSE = case_when(
                grepl("COCOA", METRIC) ~ "COCOA",
                grepl("MINING", METRIC) ~ "MINING",
                grepl("TREEPLANTATION", METRIC) ~ "TREE PLANTATION",
                grepl("SETTLEMENT", METRIC) ~ "SETTLEMENT",
                grepl("FOODCROP", METRIC) ~ "FOOD CROP",
                TRUE ~ "NA")) %>% 
              mutate(TYPE = "INDIRECT DEFORESTATION & DEGRADATION") %>% 
  mutate(TYPE = case_when(
                grepl("FOOD CROP", LANDUSE) ~ "DEFORESTATION & DEGRADATION",
                TRUE ~ TYPE)) %>% 
  group_by(LANDUSE, TYPE, SCENARIO) %>%
              mutate(AREA_HA = sum(AREA_HA)) %>%
  select(LANDUSE, TYPE, SCENARIO, AREA_HA) %>% 
   unique() 

tempc_hfz <- tempa_hfz  %>% rbind(tempb_hfz)

# Substract all the indirect deforestation & foodcrop direct def&def to "other DD"
indirectDD_hfz <- tempc_hfz %>% 
  filter(TYPE == "INDIRECT DEFORESTATION & DEGRADATION") %>% 
  group_by(SCENARIO, TYPE) %>% 
  mutate(AREA_HA = sum(AREA_HA)) %>% 
  mutate(LANDUSE = "ALL_LUI_CC") %>% 
  select(-DRIVER) %>% 
  unique()

fcDD_hfz <- tempc_hfz %>% 
  filter(LANDUSE == "FOOD CROP") %>% 
  select(-DRIVER)

fc_direct_ind_hfz <- indirectDD_hfz %>% 
  rbind(fcDD_hfz) %>% 
  group_by(SCENARIO) %>% 
  reframe(FC_totDD_AREA_HA = sum(AREA_HA)) 

timberDD_hfz <- tempc_hfz %>% 
  filter(LANDUSE == "OTHER") %>% 
  select(-SCENARIO) %>% 
  cbind(fc_direct_ind_hfz) %>% 
  group_by(SCENARIO) %>% 
  mutate(TIMBER_HA = AREA_HA - FC_totDD_AREA_HA) %>% 
  reframe(LANDUSE = "LOGGING\nAND OTHER", TYPE, SCENARIO, AREA_HA = TIMBER_HA)

tempd_hfz <- tempc_hfz %>% 
  filter(LANDUSE != "OTHER") %>% 
  rbind(timberDD_hfz)


# Plot
ggplot(data = tempd_hfz %>% 
         pivot_wider(names_from = SCENARIO, values_from = AREA_HA) %>% 
         mutate(LANDUSE_TYPE = paste0(LANDUSE, " ", TYPE)),
       aes(x = LANDUSE, y = MID_EXP, fill = LANDUSE_TYPE)) + 
    geom_bar(stat="identity", position = "stack") +
  geom_errorbar(aes(ymin = MAX_EXP, ymax = MIN_EXP), width = 0.2) +
     #facet_grid(. ~ TITLE) +
   coord_flip() +
      theme(axis.title.y = element_text(angle = 90, size=16),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 16),
        axis.text.x = element_text(size = 16), 
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 16),
        strip.background = element_rect(colour = "grey50"),
        legend.position = "none", 
        legend.text = element_blank(), 
        legend.title = element_blank(),
        panel.spacing = unit(1.5, "lines"), 
        panel.background = element_rect(fill = "white", colour = "grey50"), 
        panel.grid.major.y = element_blank(), 
        panel.grid.major.x = element_line(colour = "grey90")) +
   scale_fill_manual(values = def_scale) +
  labs(y = "Area deforested and/or degraded\nper land use (ha) 2000-2019", x = "") +
  scale_y_continuous(label = unit_format(unit = "K", scale = 1e-3, sep = ""))

```

#### fig.S6B: 2000-2019 direct and indirect deforestation and degradation per land use in percentage 

```{r,fig.height=4, fig.width=6}
# Plot percentages

# Total deforestation & degradation for each scenario 
# ... (needs to add direct def & deg to Max and Min scenarios)
# for HFZ
def_hfz_agg <- tempd_hfz %>% 
  mutate(SCENARIO =  case_when(
    SCENARIO == "MID_EXP" ~ "MID", 
    SCENARIO == "MIN_EXP" ~ "MIN", 
    SCENARIO == "MAX_EXP" ~ "MAX", 
                TRUE ~ SCENARIO)) %>% 
  rbind(tempa_hfz %>% mutate(SCENARIO = "MIN") %>% filter(LANDUSE != "OTHER")) %>%
  rbind(tempa_hfz %>% mutate(SCENARIO = "MAX")  %>% filter(LANDUSE != "OTHER")) 

tot_def_hfz_agg_2 <- def_hfz_agg %>% 
  group_by(SCENARIO) %>% 
  reframe(TOT_AREA_HA = sum(AREA_HA))

# Calculate percentage of each LU in the HFZ
sum_tot_def_at_hfz <- def_hfz_agg %>% 
  left_join(tot_def_hfz_agg_2, by = "SCENARIO") %>% 
  group_by(LANDUSE, TYPE, SCENARIO) %>% 
  # reframe(AREA_HA = sum(AREA_HA), TOT_AREA_HA) %>% 
  reframe(PERC = AREA_HA/TOT_AREA_HA *100) %>% 
  unique()



# Plot percentage of total direct & indirect def & deg per LU 
sum_tot_def_at_hfz %>% 
  mutate(LANDUSE_TYPE = paste0(LANDUSE, " ", TYPE)) %>% 
  mutate(LANDUSE_TYPE = factor(LANDUSE_TYPE, 
                                      levels = c(
                                        "TREE PLANTATION DIRECT DEFORESTATION",
                                        "SETTLEMENT DIRECT DEFORESTATION", 
                                        "MINING DIRECT DEFORESTATION",
                                        "COCOA DIRECT DEGRADATION", 
                                        "COCOA DIRECT DEFORESTATION",
                                        "LOGGING\nAND OTHER DEFORESTATION & DEGRADATION",
                                        "TREE PLANTATION INDIRECT DEFORESTATION & DEGRADATION",
                                        "SETTLEMENT INDIRECT DEFORESTATION & DEGRADATION",
                                        "MINING INDIRECT DEFORESTATION & DEGRADATION", 
                                        "COCOA INDIRECT DEFORESTATION & DEGRADATION", 
                                        "FOOD CROP DEFORESTATION & DEGRADATION"))) %>% 
        arrange(LANDUSE_TYPE) %>% 

ggplot(aes(x = SCENARIO, y = PERC, fill = LANDUSE_TYPE)) + 
    geom_bar(stat = "identity", position = "stack") + 
    coord_flip() +
      theme(axis.title.y = element_text(angle = 90, size=16),
        axis.text.y = element_text(size = 17),
        axis.title.x = element_text(size = 17),
        axis.text.x = element_text(size = 17), 
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 16),
        strip.background = element_rect(colour = "grey50"),
        legend.position = "none", 
        panel.spacing = unit(1.5, "lines"), 
        panel.background = element_rect(fill = "white", colour = "grey50"), 
        panel.grid.major.y = element_blank(), 
        panel.grid.major.x = element_line(colour = "grey90")) +

     scale_fill_manual(values = def_scale_levels) +
  labs(y = "Area deforested and/or degraded\nper land use (%)", x = "")

```


# SECTION 3.4. PROJECTING TRENDS IN THE FUTURE

#### fig.6A: Percentage of each region (restricted to HFZ) occupied by each land use in 2019-2020

```{r,fig.height=5, fig.width=13}

left_area_rg <- merged_landuses_rg_2019 %>% 
  left_join(crops_rg_2019, by = "REGION") %>% 
  mutate(LEFT_HA =  REGION_HA - LANDUSE_HA - CROP_HA) %>% 
  mutate(LEFT_HA = ifelse(LEFT_HA <0, 0, LEFT_HA)) %>% 
  mutate(LANDUSE_HA = LEFT_HA, LANDUSE = "LEFT") %>% 
  select(-REGION_HA, -LEFT_HA, -CROP_HA)

all_landuses_rg_2019 <- landuses_rg_2019 %>% 
  rbind(crops_rg_2019 %>% rename(LANDUSE_HA = CROP_HA) %>% mutate(LANDUSE = "FOOD CROP")) %>% 
  rbind(left_area_rg) %>% 
  left_join(merged_landuses_rg_2019 %>% select(REGION, REGION_HA), by = "REGION") %>% 
  mutate(PERC_OCCUP = LANDUSE_HA/REGION_HA*100)

all_landuses_2019 <- landuses_rg_2019 %>% 
  rbind(crops_rg_2019 %>% rename(LANDUSE_HA = CROP_HA) %>% mutate(LANDUSE = "FOOD CROP")) %>% 
  rbind(left_area_rg) %>% 
  left_join(merged_landuses_rg_2019 %>% select(REGION, REGION_HA), by = "REGION") %>% 
  group_by(LANDUSE) %>% 
  reframe(REGION = "HFZ", LANDUSE_HA = sum(LANDUSE_HA), REGION_HA = sum(REGION_HA)) %>% 
  mutate(PERC_OCCUP = LANDUSE_HA/REGION_HA*100) 

all_landuses_rg_tot_2019 <- all_landuses_rg_2019 %>% 
  rbind(all_landuses_2019)


# Plot percentages
ggplot(data = all_landuses_rg_tot_2019  %>% 
         mutate(LANDUSE = case_when(
                grepl("OTHERTMF", LANDUSE) ~ "OTHER TROPICAL MOIST FOREST (TMF)",
                grepl("DTMF", LANDUSE) ~ "DEGRADED TMF",
                grepl("UTMF", LANDUSE) ~ "UNDISTURBED TMF",
                grepl("LEFT", LANDUSE) ~ "OTHER LAND USE/LAND COVER",
                grepl("TREEPLANTATION", LANDUSE) ~ "TREE PLANTATION",
                grepl("SETTLEMENTS", LANDUSE) ~ "SETTLEMENT",
                TRUE ~ LANDUSE))  %>% 
         
         mutate(REGION = factor(REGION, levels= c("ASHANTI","BRONG AHAFO","CENTRAL", "EASTERN","VOLTA","WESTERN","HFZ"))) %>% 
         mutate(LANDUSE = factor(LANDUSE, levels= c("OTHER LAND USE/LAND COVER","WATER","REGROWTH", "OTHER TROPICAL MOIST FOREST (TMF)", "DEGRADED TMF","UNDISTURBED TMF","SETTLEMENT","TREE PLANTATION", "MINING", "FOOD CROP", "COCOA"))) %>% 
         mutate(TITLE = "LAND USE/COVER 2019")) + 
  geom_bar(aes(x = REGION, y = PERC_OCCUP, fill = LANDUSE), stat = "identity", position = "stack") +
  facet_grid(. ~ TITLE) +  
  coord_flip() +
      theme(axis.title.y = element_text(angle = 90, size=16),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_blank(),
        axis.text.x = element_text(size = 16), 
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 16),
        strip.background = element_rect(colour = "grey50"),
        legend.position = "right", 
        legend.text = element_text(size = 14), 
        legend.title = element_blank(),
        panel.spacing = unit(1.5, "lines"), 
        panel.background = element_rect(fill = "white", colour = "grey50"), 
        panel.grid.major.y = element_blank(), 
        panel.grid.major.x = element_line(colour = "grey90")) +
  scale_x_discrete(limits = c("ASHANTI","BRONG AHAFO","CENTRAL", "EASTERN","VOLTA","WESTERN","HFZ")) +
  scale_fill_manual(values = c(
    "#DDDDDD", #left
    "#88CCEE", #water
    "#999933", #regrowth
    "#8FBC8F", #othertmf
    "#00c88c", #dtmf
    "#008962", #utmf
    "#332288", #settlements
    "#CC6677", #tree plantations
    "#882255", #mining
    "#DDCC77", #foodcrops
    "#A0522D" #cocoa
    )) +
  labs(y = "Expansion area (%)", x = "", fill = "Occupation of\nthe region area\nper LU (%)") +
  scale_y_continuous(label = unit_format(unit = "%", sep = ""))

```


#### fig.6B: Trend over 2020-2030 of the expansion area of each land use 

```{r,fig.height=5, fig.width=13}
# Overall annual expansion of each land use per region : 2000-2019 divided by 19
# Left over of forest per region 
# "Left" land per region 

def_scale_trend0 <- c(
  "#EB5E00", #cocoa
  "#9BDBFF", #"food crops"
  "#FFD700", #mining
  "#646464", #settlement
  "#00FFFF", #tree plantations
  "#008962", #utmf 
  "#00c88c", #dtmf
  "#8FBC8F", #other tmf
  "#DDDDDD" #left
)


def_scale_trend <- c(
  "#A0522D", #cocoa
  "#DDCC77", #"food crops"
  "#882255", #mining
  "#332288", #settlement
  "#CC6677", #tree plantations
  "#008962", #utmf 
  "#00c88c", #dtmf
  "#8FBC8F", #other tmf
  "#DDDDDD" #left
)


temp_exp_10y_rg <- landuses_EA %>% 
  mutate(EA_HA_ANNUAL = EA_HA/19) %>% 
  mutate(AREA_HA = EA_HA_ANNUAL*11) %>% 
  mutate(TYPE = "ESTIMATED 2020-2030 EXPANSION AREA") %>% 
  reframe(REGION, LANDUSE, TYPE, SCENARIO, AREA_HA) %>% 
  mutate(LANDUSE = case_when(
                grepl("CROP", LANDUSE) ~ "FOOD CROP",
                TRUE ~ LANDUSE))  
  
temp_exp_10y_tot_rg <- temp_exp_10y_rg %>% 
  group_by(REGION, SCENARIO) %>% 
  reframe(AREA_HA = sum(AREA_HA), TYPE) %>% 
  mutate(LANDUSE = "TOTAL LAND USES") %>% 
  unique()

exp_10y_tot_rg <- temp_exp_10y_rg %>% rbind(temp_exp_10y_tot_rg)
  
left_landforexp_rg <- all_landuses_rg_tot_2019 %>% 
  filter(LANDUSE == "OTHERTMF" | LANDUSE == "DTMF" | LANDUSE == "UTMF" | LANDUSE == "LEFT") %>% 
  mutate(TYPE = case_when(
                grepl("OTHERTMF", LANDUSE) ~ "FOREST",
                grepl("DTMF", LANDUSE) ~ "FOREST",
                grepl("UTMF", LANDUSE) ~ "FOREST",
                grepl("LEFT", LANDUSE) ~ "LEFT",
                TRUE ~ "NA")) %>% 
  reframe(REGION, LANDUSE, TYPE, SCENARIO = "MID_EXP", AREA_HA = LANDUSE_HA) %>% 
  rbind(temp_exp_10y_rg) %>% 
  mutate(LANDUSE = factor(LANDUSE, levels= c("COCOA","FOOD CROP","MINING", "SETTLEMENT","TREEPLANTATION", "UTMF", "DTMF", "OTHERTMF", "LEFT"))) %>% 
          arrange(LANDUSE) %>% 
  filter(REGION != "HFZ") %>% 
    mutate(LANDUSE = case_when(
                grepl("OTHERTMF", LANDUSE) ~ "OTHER TMF 2019",
                grepl("DTMF", LANDUSE) ~ "DEGRADED TMF 2019",
                grepl("UTMF", LANDUSE) ~ "UNDISTURBED TMF 2019",
                grepl("LEFT", LANDUSE) ~ "OTHER LAND USE/LAND COVER",
                grepl("COCOA", LANDUSE) ~ "COCOA ESTIMATED EXPANSION 2020-2030",
                grepl("FOOD CROP", LANDUSE) ~ "FOOD CROP ESTIMATED EXPANSION 2020-2030",
                grepl("MINING", LANDUSE) ~ "MINING ESTIMATED EXPANSION 2020-2030",
                grepl("TREEPLANTATION", LANDUSE) ~ "TREE PLANTATION ESTIMATED EXPANSION 2020-2030",
                grepl("SETTLEMENT", LANDUSE) ~ "SETTLEMENT ESTIMATED EXPANSION 2020-2030",
                TRUE ~ "NA"))  %>% 
  mutate(LANDUSE = factor(LANDUSE, levels= c("COCOA ESTIMATED EXPANSION 2020-2030","FOOD CROP ESTIMATED EXPANSION 2020-2030","MINING ESTIMATED EXPANSION 2020-2030", "SETTLEMENT ESTIMATED EXPANSION 2020-2030","TREE PLANTATION ESTIMATED EXPANSION 2020-2030", "UNDISTURBED TMF 2019", "DEGRADED TMF 2019", "OTHER TMF 2019", "OTHER LAND USE/LAND COVER"))) %>% 
          arrange(LANDUSE) %>% 
  mutate(TYPE = factor(TYPE, levels= c("LEFT", "FOREST", "ESTIMATED 2020-2030 EXPANSION AREA"))) %>% 
          arrange(TYPE) 
  
# Plot
ggplot(data = left_landforexp_rg %>% 
         pivot_wider(names_from = SCENARIO, values_from = AREA_HA) %>% 
         mutate(TITLE = "TREND 2020-2030"),
       aes(x = REGION, y = MID_EXP, fill = LANDUSE, group = TYPE)) + 
    geom_bar(stat="identity", position = position_dodge(width=0.9)) +
  geom_errorbar(aes(ymin = MAX_EXP, ymax = MIN_EXP), position = position_dodge(width=0.9), width = 0.1) +
   facet_grid(. ~ TITLE) +
    coord_flip() +
      theme(axis.title.y = element_text(angle = 90, size=16),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_blank(),
        axis.text.x = element_text(size = 20), 
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 20),
        strip.background = element_rect(colour = "grey50"),
        legend.position = "right", 
        legend.text = element_text(size = 12), 
        legend.title = element_blank(),
        panel.spacing = unit(1.5, "lines"), 
        panel.background = element_rect(fill = "white", colour = "grey50"), 
        panel.grid.major.y = element_blank(), 
        panel.grid.major.x = element_line(colour = "grey90")) +
   scale_fill_manual(values = def_scale_trend) +
  labs(y = "Area (ha)", x = "", fill = "") +
  scale_y_continuous(label = unit_format(unit = "K", scale = 1e-3, sep = ""))
```

















