[Email Rachael Garrett 18/05/2021] "Our fieldwork plan is more or less to sample across these gradients in regions where cocoa is produced: - Control areas (NO Jurisdictional areas (JAs), NO protected ares (PAs), NOT in areas where Cargill or Barry Callebut control a majority of the market (ZDC) - this is because Cargill and BC seem to have the strongest commitments, but this will be adjusted with more precision as we analyze the commitment data and market shares.
Each of these, except the PAs areas also have an a and b for certified and non-certified farms (RA/Utz and/or Fair Trade), likely clustered at the village level or wherever cooperatives typically aggregate. I am not entirely sure what to do with the PA areas - I would really like to remove them as case studies - in that case they could also be considered a spillover area instead. I’m not sure how many cocoa farms we would find in those areas…"
The purpose of this script is to produce an csv file indicating for each CIV department whether it matches the above hypotheses or not. The analysis is also done on the cooperatives for which we have trade data (script : survey_site_selection_per_coop_CIV.Rmd).
cut_off_prod_perc <- 15 #percentage of cocoa beans purchased compared to the department's production
cut_off_pa_area <- 15 #percentage of protected area coverage in the department
ZDC: For a department to be considered as part of a ZDC, the share of the volume purchased by Barry Callebaut and Cargill should be a least 15% of the department’s total produced volume (TO DISCUSS).
PROTECTED AREA: Departments with a percentage of protected area coverage higher than 15 % are considered as having a higher risk of having farms located in protected areas (and/or higher risk of cocoa being purchased from protected areas) (TO DISCUSS).
# CIV departments ----------------------------
civ_departments <- s3read_using(
object = "cote_divoire/spatial/BOUNDARIES/DEPARTEMENT/OUT/CIV_DEPARTEMENTS.geojson",
FUN = read_sf,
bucket = "trase-storage",
opts = c("check_region" = T)) %>%
dplyr::select(LVL_4_CODE, LVL_4_NAME, LVL_3_CODE, LVL_3_NAME, geometry) %>%
st_transform(4326)
# Correcting an invalid polygon (self-intersection)
civ_departments <- civ_departments %>%
as("Spatial") %>%
rgeos::gBuffer(byid=TRUE, width=0) %>%
st_as_sf()
# Adding department area in km2
civ_departments <- civ_departments %>%
mutate(DEP_AREA_KM2 = drop_units(units::set_units(st_area(civ_departments), value = km^2)))
tm_shape(civ_departments) +
tm_polygons(col = "gray90", border.col= "white") +
tm_text(text = "LVL_4_NAME", size = 0.5) +
tmap_options(check.and.fix = TRUE) +
tm_layout(frame = FALSE)
Data per department on:
Relative cdd in ha/ha of forest can be added if needed.
# COCOA PRODUCTION 2019 ----------------------------
cocoa_production_pg <- s3read_using(
object = "cote_divoire/cocoa/production_estimates/UPDATED_MARCH21_jrc_cocoa_production_pg_2019.csv",
bucket = "trase-storage",
opts = c("check_region" = T),
FUN = read_delim,
delim = ";") %>%
rename(DEP_COCOA_PROD_KG_2019 = jrc_production_2019_kg)
#icco_2019_tonnes <- 2154000 #in tonnes
# COCOA AREA 2019 (JRC map) ----------------------------
cocoa_area_ha_pg <- read_delim(
get_object("cote_divoire/spatial/LULC/OUT/UPDATED_JAN21_jrc_area_under_cocoa_ha.csv",
bucket = "trase-storage",
check_region = T),
delim = ","
) %>%
dplyr::select(GEOCODE = LVL_4_CODE, COCOA_AREA_HA_2019 = sum)
# COCOA DRIVEN DEFORESTATION 2000 - 2015 (FAO / BNETD) ----------------------------
cdd <-
read_delim(
get_object("cote_divoire/cocoa/cocoa_driven_deforestation/OUT/past_cocoa_deforestation.csv",
bucket = "trase-storage",
check_region = T),
delim = ","
) %>%
dplyr::select(LVL_4_CODE, CDD_HA = sum) %>%
left_join(cocoa_production_pg, by = "LVL_4_CODE") %>%
#Adding relative CDD (compared to production)
mutate(CDD_HA_TONNE = (CDD_HA) / (DEP_COCOA_PROD_KG_2019 / 1000),
CDD_HA_TONNE = ifelse(is.na(CDD_HA_TONNE), 0, CDD_HA_TONNE)) %>%
#Annualized deforestation (15 year-deforestation period divided by 15 to get ha/year)
mutate(CDD_HA_TONNE_YEAR = (CDD_HA / 15) / (DEP_COCOA_PROD_KG_2019 / 1000),
CDD_HA_TONNE_YEAR = ifelse(is.na(CDD_HA_TONNE_YEAR), 0, CDD_HA_TONNE_YEAR)) %>%
left_join(cocoa_area_ha_pg, by = c("LVL_4_CODE" = "GEOCODE"))
cdd_sf <- civ_departments %>% # turned into sf object
left_join(cdd, by = "LVL_4_CODE") %>%
#Adding relative CDD (compared to department area)
mutate(CDD_HA_DEP_AREA = CDD_HA / DEP_AREA_KM2)
# Plot Cocoa Area 2019
tm_shape(civ_departments) +
tm_polygons(col = "gray90") +
tmap_options(check.and.fix = TRUE) +
tm_shape(cdd_sf) +
tm_polygons(title = "Cocoa Area in 2019 (ha) (Source: JRC)", col = "COCOA_AREA_HA_2019") +
tm_scale_bar() +
tm_layout(legend.outside = T, frame = FALSE)
# Plot Cocoa Production 2019
tm_shape(civ_departments) +
tm_polygons(col = "gray90") +
tmap_options(check.and.fix = TRUE) +
tm_shape(cdd_sf) +
tm_polygons(title = "Cocoa Production in 2019 (kg) \n (Source: Trase based on JRC)", col = "DEP_COCOA_PROD_KG_2019") +
tm_scale_bar() +
tm_layout(legend.outside = T, frame = FALSE)
# Plot Cocoa Driven Deforestation 2000 - 2015
# CDD
tm_shape(civ_departments) +
tm_polygons(col = "gray90") +
tmap_options(check.and.fix = TRUE) +
tm_shape(cdd_sf) +
tm_polygons(title = "CDD 2000 - 2015 (ha) \n (Source: Trase based on FAO/BNETD)", col = "CDD_HA") +
tm_scale_bar() +
tm_layout(legend.outside = T, frame = FALSE)
# Plot CDD relative to cocoa production
tm_shape(civ_departments) +
tm_polygons(col = "gray90") +
tmap_options(check.and.fix = TRUE) +
tm_shape(cdd_sf) +
tm_polygons(title = "Relative CDD 2000 - 2015 (ha/ton) \n (Source: Trase based on FAO/BNETD)", col = "CDD_HA_TONNE") +
tm_scale_bar() +
tm_layout(legend.outside = T, frame = FALSE)
# Plot CDD relative to department's surface
tm_shape(civ_departments) +
tm_polygons(col = "gray90") +
tmap_options(check.and.fix = TRUE) +
tm_shape(cdd_sf) +
tm_polygons(title = "Relative CDD 2000 - 2015 (ha/km2) \n as compared to departement's area \n (Source: Trase based on FAO/BNETD)", col = "CDD_HA_DEP_AREA") +
tm_scale_bar() +
tm_layout(legend.outside = T, frame = FALSE)
# TERRITORIAL DEFORESTATION ----------------------------
# Deforestation per department in 2019 (total deforestation, not only cocoa-related)
territorial_deforestation <- read_delim(
get_object("cote_divoire/cocoa/indicators/out/territorial_deforestation_q2_2021.csv",
bucket = "trase-storage",
check_region = T
),
delim = ";"
) %>%
transmute(GEOCODE, TERRITORIAL_DEFORESTATION_HA_2019 = TERRITORIAL_DEFORESTATION_HA)
territorial_defo_sf <- civ_departments %>%
left_join(territorial_deforestation, by = c("LVL_4_CODE" = "GEOCODE")) %>%
#Adding relative territorial deforestation (compared to department area)
mutate(TERR_DEFO_HA_DEP_AREA = TERRITORIAL_DEFORESTATION_HA_2019 / DEP_AREA_KM2)
# Plot Territorial deforestation
tm_shape(civ_departments) +
tm_polygons(col = "gray90") +
tmap_options(check.and.fix = TRUE) +
tm_shape(territorial_defo_sf) +
tm_polygons(title = "Total deforestation in 2019 (ha)", col = "TERRITORIAL_DEFORESTATION_HA_2019") +
tm_scale_bar() +
tm_layout(legend.outside = T, frame = FALSE)
# Plot Territorial deforestation relative to department's surface
tm_shape(civ_departments) +
tm_polygons(col = "gray90") +
tmap_options(check.and.fix = TRUE) +
tm_shape(territorial_defo_sf) +
tm_polygons(title = "Relative total deforestation 2019 (ha/km2) \n as compared to departement's area", col = "TERR_DEFO_HA_DEP_AREA") +
tm_scale_bar() +
tm_layout(legend.outside = T, frame = FALSE)
# PROTECTED AREAS (ALL CATEGORIES) ---------------------------------------
# We use the BNETD map - for the comparison of various maps, see the script per cooperatives (survey_site_selection_per_coop_CIV.Rmd).
civ_protected_areas <- s3read_using(
object = "cote_divoire/spatial/BOUNDARIES/PROTECTED_AREAS/BNETD/OUT/BNETD_PROTECTED_AREAS_CIV.gpkg",
FUN = read_sf,
bucket = "trase-storage",
opts = c("check_region" = T)) %>%
st_transform(4326)
# Plot results
tm_shape(civ_departments) +
tm_polygons(col = "gray90") +
tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas) +
tm_polygons(border.col = NULL, col = "forestgreen") +
tm_scale_bar() +
tm_layout(legend.outside = T, frame = FALSE)
As we are working at the department level, an initial idea was to look at which departments intersect with protected areas, but as many do, to get a finer filtering, we can :
(TO DISCUSS) Decision is to keep option 1 but if option 2 is relevant, it was already done and can be added.
# Intersection between departments and protected areas
pa_pg <- st_intersection(civ_departments, civ_protected_areas)
# Calculate the total surface of protected area polygons in each department
pa_pg_area <- pa_pg %>%
mutate(PA_AREA_HA = drop_units(units::set_units(st_area(geometry), value = hectare))) %>%
group_by(LVL_4_CODE) %>%
summarize(PA_AREA_HA = sum(PA_AREA_HA))
# Calculate the area of each department and calculate the percentage that the protected area surface represents compared to the department area
civ_departments_area <- civ_departments %>%
transmute(LVL_4_CODE, DEP_AREA_HA= drop_units(units::set_units(st_area(geometry), value = hectare))) %>%
left_join(st_set_geometry(pa_pg_area, NULL), by = "LVL_4_CODE") %>%
mutate(DEP_PERC_OF_PA = PA_AREA_HA / DEP_AREA_HA * 100)
# Threshold to be discussed. 1 = Close to a protected area, 0 = far from a protected area
civ_departments_pa_area <- civ_departments_area %>%
mutate(
DEP_HIGH_PA_COVER = case_when(
DEP_PERC_OF_PA >= cut_off_pa_area ~ "1",
DEP_PERC_OF_PA < cut_off_pa_area ~ "0",
TRUE ~ "0"
))
# Plot results
tm_shape(civ_departments) +
tm_polygons(col = "gray90") +
tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas) +
tm_polygons(border.col = NULL, col = "forestgreen") +
tm_shape(civ_departments_pa_area) +
tm_polygons(col = "DEP_HIGH_PA_COVER", title = "Protected area coverage (%) \n 1 = High, 0 = Low", alpha = 0.4) +
tm_scale_bar() +
tm_layout(legend.outside = T, frame = FALSE)
A threshold of 15 % is fixed to segregate departments considered as having a low protected area coverage (< 15 % of the department area) and those with high protected area coverage (> 15 %). The analysis at the cooperative level shows which cooperatives are close to protected areas.
The CFI landscapes in Cote d’Ivoire are 5 priority regions covering the largest remaining forests. They are Guémon, Cavally, Nawa, San-Pedro and La Mé.
# CFI LANDSCAPES (JURISDICTIONAL APPROACH) ---------------------------------------
cfi_landscapes_civ <- civ_departments %>%
filter(LVL_3_NAME %in% c("CAVALLY", "SAN-PEDRO", "NAWA", "ME", "GUEMON")) %>%
dplyr::select(LVL_3_NAME, LVL_3_CODE)
# Plot CFI landscape
tm_shape(civ_departments) +
tm_polygons(col = "gray90") +
tmap_options(check.and.fix = TRUE) +
tm_shape(cfi_landscapes_civ) +
tm_polygons(title ="CFI landscapes", col = "LVL_3_NAME") +
tm_scale_bar() +
tm_layout(legend.outside = T, frame = FALSE) +
tm_shape(civ_protected_areas) +
tm_polygons(border.col = NULL, col = "forestgreen", alpha = 0.5)
Assumption: knowing that Cargill & BC have more stringent ZDC commitments (according to ETH), it is assumed that departments from which they source their cocoa (threshold to be defined) are under a ZDC.
To be discussed: we could consider various ways of defining whether a department is part of a ZDC or not:
(TO DISCUSS) Decision is to keep option 1 but if option 2 is relevant, it was already done and can be added.
# ZDC: VOLUME SOURCED BY CARGILL & BARRY CALLEBAUT PER DEPARTMENT 2019 ------------
# all volumes exported per exporter, per destination country, per cocoa seller
volume_fob_exporter_country_provider <-
s3readRDS(
"cote_divoire/cocoa/sei-pcs/OUT/volume_fob_exporteroriginal_country_producer_2019.rds",
bucket = "trase-storage"
) %>%
rename(GEOCODE = GEOCODE_PROVIDER) %>%
filter(! GEOCODE == "UNKNOWN") %>%
dplyr::select(-GEOCODE_COCOA)
volume_fob_exporter_country_provider <- volume_fob_exporter_country_provider %>%
mutate(COOP_SOURCED_DIRECT_VOLUME_BEQ_MEAN = map_dbl(COOP_SOURCED_DIRECT_VOLUME_BEQ_distribution,
~ mean(.x$COOP_SOURCED_DIRECT_VOLUME_BEQ, na.rm = T))) # calculate the mean of the estimates of the volume sourced by trader from the coops
#Note on the above line: there is uncertainty in the volumes sourced per cooperative, which we simulate using Monte Carlo sampling of the available data. The 1000 estimates are contained in COOP_SOURCED_DIRECT_VOLUME_BEQ_distribution.
#We reflect this uncertainty by plotting the confidence intervals (e.g. in the bar charts), drawn from the resulting distribution, but when plotting maps, we extract the mean of the distribution (i.e. the mean estimate of sourcing per COOP).
# Total sourced by trader, per GEOCODE (pg)
exp_pg_beq_trader <- volume_fob_exporter_country_provider %>%
filter(EXPORTER_CLEAN %in% c("CARGILL", "BARRY CALLEBAUT")) %>%
group_by(EXPORTER_CLEAN, GEOCODE) %>%
summarise(BEQ_TRADER = sum(COOP_SOURCED_DIRECT_VOLUME_BEQ_MEAN))
exp_pg_per_trader <- exp_pg_beq_trader %>%
split(exp_pg_beq_trader$EXPORTER_CLEAN) #split the table into two (one with Cargill, one with Barry), to keep only one line per department
# Option 1: Calculate the percentage that the sourcing of BARRY CALLEBAUT & CARGILL represents compared to the Department's production
exp_pg_per_trader_joined <- civ_departments %>%
left_join(cdd, by = "LVL_4_CODE") %>%
left_join(
transmute(exp_pg_per_trader$`BARRY CALLEBAUT`, GEOCODE, TOTAL_BEQ_BARRY = BEQ_TRADER),
by = c("LVL_4_CODE" = "GEOCODE")) %>%
left_join(
transmute(exp_pg_per_trader$`CARGILL`, GEOCODE, TOTAL_BEQ_CARGILL = BEQ_TRADER),
by = c("LVL_4_CODE" = "GEOCODE"))
exp_pg_per_trader_joined$TOTAL_BEQ_BARRY[is.na(exp_pg_per_trader_joined$TOTAL_BEQ_BARRY)] <- 0
exp_pg_per_trader_joined$TOTAL_BEQ_CARGILL[is.na(exp_pg_per_trader_joined$TOTAL_BEQ_CARGILL)] <- 0
exp_pg_per_trader_joined <- exp_pg_per_trader_joined %>%
mutate(TOTAL_BEQ_BARRY_CARG = TOTAL_BEQ_BARRY + TOTAL_BEQ_CARGILL) %>%
mutate(PERC_BARRY_CARG_PROD = TOTAL_BEQ_BARRY_CARG / DEP_COCOA_PROD_KG_2019 * 100)
# Threshold to be discussed (section 1). 1 = YES (part of a ZDC), 0 = NO (not part of ZDC)
exp_pg_per_trader_joined <- exp_pg_per_trader_joined %>%
mutate(
DEP_IN_ZDC = case_when(
PERC_BARRY_CARG_PROD >= cut_off_prod_perc ~ "1",
PERC_BARRY_CARG_PROD < cut_off_prod_perc ~ "0",
TRUE ~ "0"
))
Here the threshold applied on the production percentage is 15% (with 1 for percentages higher or equal to 15%, and 0 for percentages below (and for departments that are not supplying BC & C, and/or are not producing cocoa)), which means that departments from which Cargill and Barry Callebaut purchase at least 15% of the total production are considered as part of a ZDC.
# PLOT ZDC INDICATORS PER GEOCODE ----------------------------
# Plot percentage of the departments' production sourced by BC & C
tm_shape(civ_departments) +
tm_polygons(col = "gray90") +
tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas) +
tm_polygons(border.col = NULL, col = "forestgreen") +
tm_shape(exp_pg_per_trader_joined) +
tm_polygons(col = "PERC_BARRY_CARG_PROD", title = "Share of production directly sourced \n by Barry Callebaut & Cargill (%)", alpha = 0.4) +
tm_scale_bar() +
tm_layout(legend.outside = T, frame = FALSE)
# Plot percentage of the departments' production sourced by BC & C (classes)
tm_shape(civ_departments) +
tm_polygons(col = "gray90") +
tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas) +
tm_polygons(border.col = NULL, col = "forestgreen") +
tm_shape(exp_pg_per_trader_joined) +
tm_polygons(col = "DEP_IN_ZDC", title = "Departments part of a ZDC \n (i.e. supplying Barry Callebaut and/or Cargill \n above a defined threshold)", alpha = 0.4) +
tm_scale_bar() +
tm_layout(legend.outside = T, frame = FALSE)
# MERGE INDICATORS PER GEOCODE ----------------------------
indicators_pg <- exp_pg_per_trader_joined %>%
left_join(territorial_deforestation, by = c("LVL_4_CODE" = "GEOCODE")) %>%
left_join(st_set_geometry(
transmute(civ_departments_pa_area, LVL_4_CODE, DEP_PERC_OF_PA, DEP_HIGH_PA_COVER), NULL),
by = "LVL_4_CODE") %>%
mutate(DEP_IN_CFI_LANDSCAPE = ifelse(
LVL_3_NAME %in% c("CAVALLY", "SAN-PEDRO", "NAWA", "ME", "GUEMON"), "1", "0")) %>% #1 = YES, 0 = NO
filter(COCOA_AREA_HA_2019 != 0) # removing departments with 0 cocoa area in 2019
# ALL HYPOTHESES
indicators_pg_hypotheses <- indicators_pg %>%
mutate(HYPOTHESES = case_when(
DEP_IN_CFI_LANDSCAPE == "0" & DEP_HIGH_PA_COVER == "0" & DEP_IN_ZDC == "0" ~ "CONTROL",
DEP_IN_CFI_LANDSCAPE == "1" & DEP_HIGH_PA_COVER == "0" & DEP_IN_ZDC == "0" ~ "H1",
DEP_IN_CFI_LANDSCAPE == "0" & DEP_HIGH_PA_COVER == "1" & DEP_IN_ZDC == "0" ~ "H2",
DEP_IN_CFI_LANDSCAPE == "0" & DEP_HIGH_PA_COVER == "0" & DEP_IN_ZDC == "1" ~ "H3",
DEP_IN_CFI_LANDSCAPE == "1" & DEP_HIGH_PA_COVER == "0" & DEP_IN_ZDC == "1" ~ "H4",
DEP_IN_CFI_LANDSCAPE == "1" & DEP_HIGH_PA_COVER == "1" & DEP_IN_ZDC == "1" ~ "H5",
TRUE ~ as.character(NA)
)
)
# Control areas : NO Jurisdictional areas (JAs), NO protected ares (PAs), NOT in areas where Cargill or Barry Callebut control a majority of the market (ZDC)
ncontrol <- indicators_pg_hypotheses %>%
filter(HYPOTHESES == "CONTROL") %>%
nrow()
nH1 <- indicators_pg_hypotheses %>%
filter(HYPOTHESES == "H1") %>%
nrow()
nH2 <- indicators_pg_hypotheses %>%
filter(HYPOTHESES == "H2") %>%
nrow()
nH3 <- indicators_pg_hypotheses %>%
filter(HYPOTHESES == "H3") %>%
nrow()
nH4 <- indicators_pg_hypotheses %>%
filter(HYPOTHESES == "H4") %>%
nrow()
nH5 <- indicators_pg_hypotheses %>%
filter(HYPOTHESES == "H5") %>%
nrow()
nMissing <- indicators_pg_hypotheses %>%
filter(is.na(HYPOTHESES)) %>%
nrow()
# Plot departments according to hypotheses
tm_shape(civ_departments) +
tm_polygons(col = "white") +
tmap_options(check.and.fix = TRUE) +
tm_shape(indicators_pg_hypotheses) +
tm_polygons(col = "HYPOTHESES", title = "Departments matching \n selection hypotheses", palette = c("blue", "red", "purple", "orange", "green", "yellow")) +
tm_scale_bar() +
tm_layout(legend.outside = T, frame = FALSE)
This gives us:
28 departments as Control (NO Jurisdictional areas (JAs), NO protected ares (PAs), NOT in areas where Cargill or Barry Callebut control a majority of the market (ZDC)
3 departments for H1, effect of JAs: NO PAs, NOT ZDC
17 departments for H2, effect of PAs: NO JAs, NOT ZDC
10 departments for H3, effect of ZDCs: NO JAs, NO PAs, ARE in ZDC areas
2 departments for H4, synergistic effect of ZDCs + JAs, NO PAs, ARE in areas with JAs and ZDCs
6 departments for H5, synergistic effect of ZDCs + JAs + PAs, ARE in areas with JAs, PAs and ZDCs
14 departments do not match any of the hypotheses (Missing), those are areas with JAs and PAs and NOT ZDCs as well as areas with NO JAs, with PAs and ZDCs.
Departments with no cocoa area in 2019 were removed. For deforestation, the information is in the table and can be filtered or added to the above selection.
s3write_using(indicators_pg_hypotheses,
object = "cote_divoire/cocoa/indicators/out/indicators_dep_HHsurvey_hypotheses.csv",
bucket = "trase-storage",
FUN = write_delim, delim = ";",
opts = c("check_region" = T)
)
indicators_dep_HHsurvey_hypotheses <- s3read_using(
object = "cote_divoire/cocoa/indicators/out/indicators_dep_HHsurvey_hypotheses.csv",
bucket = "trase-storage",
FUN = read_delim,
delim = ";",
opts = c("check_region" = T)
)