Fieldwork - data and approach

[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.

H1, effect of JAs: NO PAs, NOT ZDC
H2, effect of PAs: NO JAs, NOT ZDC
H3, effect of ZDCs: NO JAs, NO PAs, ARE in ZDC areas
H4, synergistic effect of ZDCs + JAs, NO PAs, ARE in areas with JAs and ZDCs
H5, synergistic effect of ZDCs + JAs + PAs, ARE in areas with JAs, PAs and ZDCs

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 cooperative (for which we have trade data) whether it matches the above hypotheses or not.

The analysis here is done on the cooperatives for which we have trade data, but it was also done at the department level (e.g. if we want to consider other cooperatives or locations) (script : survey_site_selection_per_dep_CIV.Rmd).

In this script, some indicators are calculated at the cooperative level, and others at the department level:

Protected areas (PAs) : cooperative level
Zero Deforestation commitment (ZDC) : cooperative level (also available at the department level)
Certification : cooperative level
Jurisdictional areas (JAs) (i.e. CFI landscapes) : department level
Additional information on deforestation and cocoa driven deforestation (not required initially in the above hypotheses, but included in the excel as it could be useful) : department level

Loading data

Setting thresholds

cut_off_buffer_pa <- 25000 #Buffer for protected areas

PROTECTED AREA: For a cooperative to be considered as close to a protected area, i.e. higher risk of having farms of members located in protected areas (and/or higher risk of buying cocoa coming from protected areas), it should be located within a buffer of 2.510^{4}m around protected areas (TO DISCUSS)

CIV departments

# 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) # the cooperative dots do not appear on the map when using the same crs as civ departments (EPSG 32630) - I don't know why. So I rather change the crs of civ departments into EPSG 4326 (WGS84)

# 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)

Cooperatives

# Cocoa accountability map - to get the coordinates of the cooperatives
cam <- s3read_using(
  object = "cote_divoire/logistics/mighty/accountability_map_v2/cam_2021-03-02_geocoded.csv",
  bucket = "trase-storage",
  opts = c("check_region" = T),
  FUN = read_delim,
  delim = ";") 

cam_coord <- cam[which(!is.na(cam$latitude) & !is.na(cam$longitude)),]
cam_coord <- cam_coord %>% 
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>% 
  dplyr::select(coopacronym, coopname) 

# Problem with duplicates (coop with same acronym & denomination but different location)
coop_name_duplicates <- cam_coord %>% 
 get_dupes(coopacronym, coopname) # 20 duplicates - coordinates are different, but keeping them might generate issues (e.g. in group_by), so they are removed

cam_coord_dupesout <- cam_coord %>% 
 distinct(coopacronym, coopname, .keep_all = TRUE)

number_all_coops <- nrow(cam_coord_dupesout)

# Keep cooperatives for which we have trade data (source : Trase)
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)

#  Joining the latitude / Longitude that we have from the CAM table 
volume_fob_exporter_country_provider_coord <- inner_join(volume_fob_exporter_country_provider, cam_coord_dupesout, by = c("COOPACRONYM_AANAME" = "coopacronym", "COOPNAME_AADENOMINATION" = "coopname"))

# Keeping only information on cooperatives and traders
coop_vol_exporter_coord <- volume_fob_exporter_country_provider_coord %>% 
  dplyr::select(
    YEAR, 
    GEOCODE, 
    COOPACRONYM_AANAME, 
    COOPNAME_AADENOMINATION, 
    UTZ, 
    RA, 
    FAIRTRADE, 
    COCOALIFE, 
    SUSTAINABLE_ORIGINS, 
    CERTIFIED, 
    EXPORTER_CLEAN, 
    geometry) %>% 
  distinct(COOPACRONYM_AANAME, COOPNAME_AADENOMINATION, EXPORTER_CLEAN, .keep_all = TRUE) 

number_coops_with_trade_info <- nrow(distinct(coop_vol_exporter_coord, COOPACRONYM_AANAME, COOPNAME_AADENOMINATION, .keep_all = FALSE))

# Turn it into a sf object
coop_vol_exporter_coord_sf <- coop_vol_exporter_coord %>% 
  st_as_sf(crs = 4326) 
#st_crs(coop_vol_exporter_coord_sf)

# List of cooperatives 
coop_coord_sf <- coop_vol_exporter_coord_sf %>% 
  distinct(COOPACRONYM_AANAME, COOPNAME_AADENOMINATION, geometry) 

# Plot total CAM cooperatives and coops for which we have Trade information 
tm_shape(civ_departments) + 
  tm_polygons(col = "gray90", border.col= "white") + 
  tmap_options(check.and.fix = TRUE) +
  tm_layout(legend.outside = T, frame = FALSE) +
tm_shape(cam_coord_dupesout) +
  tm_dots(col = "coral1", size = 0.1) +
tm_shape(coop_coord_sf) +
  tm_dots(col = "steelblue", size = 0.1) +
  tm_add_legend(title = "Cooperatives \n (Source: Cocoa Accountability Map)", labels = c("Coops without trade info", "Coops with trade info"), col = c("coral1", "steelblue"), border.col = NULL)

The Cocoa Accountability Map (CAM) has 4253 cooperatives. After joining it with the trase excel file (volume_fob_exporter_country_provider) to keep only cooperatives for which we have trade data, we are left with 593 cooperatives.

Cocoa production, area & cocoa driven deforestation

Data per department on:

Cocoa production in 2019 in kg from Trase (DEP_COCOA_PROD_KG_2019),
Cocoa area in 2019 from JRC map (COCOA_AREA_HA_2019)
Cocoa driven deforestation in ha: areas deforested between 2000 and 2015 where cocoa farms are detected in 2019 (CDD_HA)
Relative cocoa driven deforestation in ha/tonnes: cocoa driven deforestation divided by the production in 2019 (CDD_HA_TONNE)
Annualized relative cocoa driven deforestation in ha/tonnes.year: relative ccd divided by 15 years (CDD_HA_TONNE_YEAR)
Relative cococa driven deforestation compared to the department area in ha/km2 (CDD_HA_DEP_AREA)

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)

Deforestation

# 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

# PROTECTED AREAS (ALL CATEGORIES)  ---------------------------------------

# All categories WDPA 
civ_protected_areas_wdpa <- s3read_using(
  object = "cote_divoire/spatial/BOUNDARIES/PROTECTED_AREAS/WDPA/OUT/WDPA_PROTECTED_AREAS_CIV_AUG21.gpkg",
  FUN = read_sf,
  bucket = "trase-storage",
  opts = c("check_region" = T)) %>% 
  st_transform(4326) # making sure that the CRS is the same as the one used for civ_departments and coops

civ_protected_areas_wdpa$IUCN_CAT[civ_protected_areas_wdpa$IUCN_CAT == "Not Reported"] <- NA
civ_protected_areas_wdpa$IUCN_CAT[civ_protected_areas_wdpa$IUCN_CAT == "Not Applicable"] <- NA

# Forêts classées SST GEOPORTAIL
civ_forets_classees_sst <- s3read_using(
  object = "cote_divoire/spatial/BOUNDARIES/PROTECTED_AREAS/SST_GEOPORTAIL/OUT/SST_FORETS_CLASSEES_CIV_AUG21.gpkg",
  FUN = read_sf,
  bucket = "trase-storage",
  opts = c("check_region" = T)) %>% 
  st_transform(4326)

# Forêts classées, Parcs and Réserves BNETD
civ_forets_classees_bnetd <- s3read_using(
  object = "cote_divoire/spatial/BOUNDARIES/PROTECTED_AREAS/BNETD/OUT/BNETD_FORETS_CLASSEES_CIV.gpkg",
  FUN = read_sf,
  bucket = "trase-storage",
  opts = c("check_region" = T)) %>% 
  st_transform(4326) %>% 
  transmute(NAME = NOM, CATEGORY = NATURE, IUCNCAT = NA, SUP_HA, YEAR = NA) %>% 
  mutate(CATEGORY = "Forêt classée")

civ_parcs_reserves_bnetd <- s3read_using(
  object = "cote_divoire/spatial/BOUNDARIES/PROTECTED_AREAS/BNETD/OUT/BNETD_PARCS_RESERVES_CIV.gpkg",
  FUN = read_sf,
  bucket = "trase-storage",
  opts = c("check_region" = T)) %>% 
  st_transform(4326) %>% 
  transmute(NAME = AREANAME, CATEGORY = DESIGNATE, IUCNCAT, SUP_HA = SIZE, YEAR)

civ_protected_areas <- rbind(civ_parcs_reserves_bnetd, civ_forets_classees_bnetd)
#write_sf(civ_protected_areas, dsn = "BNETD_PROTECTED_AREAS_CIV.gpkg")

# Plot protected areas WDPA and Forêts classées SST
tm_shape(civ_departments) + 
  tm_polygons(col = "gray90") + 
  tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas_wdpa) + 
  tm_polygons(border.col = NULL, col = "forestgreen") +
  tm_scale_bar() + 
  tm_layout(legend.outside = T, frame = FALSE) +
tm_shape(civ_forets_classees_sst) +
  tm_polygons(border.col = NULL, col = "coral1", alpha = 0.6) +
  tm_add_legend(labels = c("All protected areas (WDPA)", "Forets classees (SST)"), col = c("forestgreen", "coral1"), border.col = NULL)

# Plot protected areas BNETD
tm_shape(civ_departments) + 
  tm_polygons(col = "gray90") + 
  tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas_wdpa) + 
  tm_polygons(border.col = NULL, col = "gray20") +
  tm_scale_bar() + 
  tm_layout(legend.outside = T, frame = FALSE) +
tm_shape(civ_forets_classees_bnetd) +
  tm_polygons(border.col = NULL, col = "coral1", alpha = 0.6) +
tm_shape(civ_parcs_reserves_bnetd) +
  tm_polygons(border.col = NULL, col = "forestgreen", alpha = 0.6) +
  tm_add_legend(labels = c("All protected areas (WDPA)", "Forets classees (BNETD)", "Parcs & Reserves (BNETD)"), col = c("gray20", "coral1", "forestgreen"), border.col = NULL) 

# Plot Forêts classées SST and BNETD
tm_shape(civ_departments) + 
  tm_polygons(col = "gray90") + 
  tmap_options(check.and.fix = TRUE) +
tm_shape(civ_forets_classees_bnetd) + 
  tm_polygons(border.col = NULL, col = "forestgreen") +
  tm_scale_bar() + 
  tm_layout(legend.outside = T, frame = FALSE) +
tm_shape(civ_forets_classees_sst) +
  tm_polygons(border.col = NULL, col = "coral1", alpha = 0.6) +
  tm_add_legend(labels = c("Forets classees (BNETD)", "Forets classees (SST)"), col = c("forestgreen", "coral1"), border.col = NULL)

# Plot protected areas WDPA & reported IUCN Categories
tm_shape(civ_departments) + 
  tm_polygons(col = "gray90") + 
  tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas_wdpa) + 
  tm_polygons(border.col = NULL, col = "IUCN_CAT", palette = c("forestgreen", "green", "steelblue"), title = "IUCN Category (WDPA)") +
  tm_text(text = "IUCN_CAT", size = 0.5) + 
  tm_scale_bar() + 
  tm_layout(legend.outside = T, frame = FALSE)

# Plot protected areas BNETD & reported Categories
tm_shape(civ_departments) + 
  tm_polygons(col = "gray95") + 
  tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas) + 
  tm_polygons(border.col = NULL, col = "CATEGORY", title = "Protected area category (BNETD)") +
  tm_scale_bar() + 
  tm_layout(legend.outside = T, frame = FALSE)

We have three data sources :

WDPA : all categories (parcs, reserves, forêts classées). This map identifies more areas as protected compared to BNETD. It might be because the map was not updated after the revision of the code forestier in 2019. Most protected areas do not have a reported IUCN category - cf. map above and table below.
SST (Geoportail of the Surveillance Spatiale des Terres) : Forêts classées only (could not find a shp for parcs & reserves on this website). It does not match exactly the boundaries of the WDPA protected areas.
BNETD (sent by Nitidae, but not sure if this map takes in account the revision of the code forestier) : all categories in two shapefiles, one with forêts classées and one with Parcs & Réserves. Variations can be observed, with areas indicated as protected in WDPA shp that are not present in BNETD. As compared to SST, some areas are indicated as forêts classées by SST and not by BNETD and vice-versa.

(TO DISCUSS) Decision is to use the BNETD map as it seems more up to date, but it would be good to have a confirmation.

# CIV Protected areas IUCN categories
iucn_cat_civ <- as.data.frame(civ_protected_areas_wdpa) %>% 
  get_dupes(IUCN_CAT) %>% 
  distinct(IUCN_CAT, dupe_count, .keep_all = FALSE) %>% 
  rename("Count" = "dupe_count")
iucn_cat_civ

##   IUCN_CAT Count
## 1       Ia     3
## 2       II     9
## 3       IV     3
## 4     <NA>   228

Mighty Earth have in the past used the following: The average distance a cocoa bean travels before it reaches a cooperative is in the range of 25-46 km. Using this figure, cooperatives were categorized into small (0-800 members), medium (800-1000 members) and large (1000+ members). The average distance a cocoa bean travels was then applied to these categories, with 25 km, 35 km and 45 km assigned to small, medium and large cooperatives respectively.

(TO DISCUSS) To identify cooperatives that are close to protected area, i.e. higher risk of having farms of members located in protected areas (and/or higher risk of buying cocoa coming from protected areas), a buffer of 2.510^{4} m is applied, otherwise almost the entire country falls into the buffer. It could be further decreased as many coops fall within the buffer. An other option would be to categorize cooperatives according to the estimated number of farmers and apply the corresponding bean travel distance.

(TO DISCUSS) All categories of protected areas are taken in account - but it could be interesting to distinguish the various IUCN categories and/or to look at the forested areas rather than protected areas.

# Apply a buffer around all protected areas
civ_protected_areas_buffer <- st_buffer(civ_protected_areas, dist = cut_off_buffer_pa) 

# Select cooperatives located in and around protected areas (within buffer)
selection <- st_intersects(coop_coord_sf, civ_protected_areas_buffer)
coop_close_to_pas <- coop_coord_sf[lengths(selection) > 0, ] # lengths() > 0  converts the sparse geometry binary predicate (sgbp) into a logical vector with TRUE for length > 0 (i.e. intersecting a protected area)
rm(selection)
coop_close_to_pas <- coop_close_to_pas %>% 
  mutate(COOP_IN_BUFFER_PA = "1")

# Plot cooperatives located in and around protected areas (within buffer)
tm_shape(civ_departments) + 
  tm_polygons(col = "gray90") + 
  tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas_buffer) + 
  tm_polygons(border.col = NULL, col = "red", alpha = 0.4) +
  tm_scale_bar() + 
  tm_layout(legend.outside = T, frame = FALSE) + 
tm_shape(civ_protected_areas) + 
  tm_polygons(border.col = NULL, col = "forestgreen") +
  tm_scale_bar() + 
  tm_layout(legend.outside = T, frame = FALSE) + 
tm_shape(coop_coord_sf) +
  tm_dots(col = "steelblue", size = 0.1) +
tm_shape(coop_close_to_pas) +
  tm_dots(col = "navy", size = 0.1) +
  tm_add_legend(title = "Proximity to Protected Areas ", labels = c("Coops inside/close to PAs", "Coops far from PAs", "Buffer"), col = c("navy", "steelblue", "red"), border.col = NULL)

CFI Landscapes

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)

# Select cooperatives located in the CFI landscapes
selection <- st_intersects(coop_coord_sf, cfi_landscapes_civ)
coop_in_cfi <- coop_coord_sf[lengths(selection) > 0, ] # lengths() > 0  converts the sparse geometry binary predicate (sgbp) into a logical vector with TRUE for length > 0 (i.e. intersecting a protected area)
rm(selection)
coop_in_cfi <- coop_in_cfi %>% 
  mutate(COOP_IN_CFI_LANDSCAPE = "1")

# Plot cooperatives located within CFI landscapes
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) + 
tm_shape(coop_coord_sf) +
  tm_dots(col = "steelblue", size = 0.1) +
tm_shape(coop_in_cfi) +
  tm_dots(col = "navy", size = 0.1)

Zero Deforestation Commitments (ZDC)

(TO DISCUSS) Assumption is that Cargill & BC have more stringent ZDC commitments, it is therefore assumed that departments from which they source their cocoa are under a ZDC.

As the data on estimated volumes produced by cooperative should be considered with error margins, it was decided to rather only show cooperatives supplying Cargill & BC, with 1 if Cargill or BC is a known buyer, 0 if not, instead of considering volumes.

# ZDC: SOURCING BY CARGILL & BARRY CALLEBAUT PER COOPERATIVE 2019 ------------

# Do Barry Callebaut and/or Cargill source from the cooperative? 1 = Yes or 0 = No 
coop_exporter_BC_C_coord_sf <- coop_vol_exporter_coord_sf %>%
  group_by(COOPACRONYM_AANAME, COOPNAME_AADENOMINATION) %>%
  # Barry Callebaut and/or Cargill are part of the buyers of the cooperative?
  mutate(COOP_IN_ZDC = case_when(
      EXPORTER_CLEAN %in% c("BARRY CALLEBAUT", "CARGILL") ~ "1", 
      TRUE ~ "0"
  )) %>% 
  # Remove duplicates (lines are merged to have only one line per cooperative)
  distinct(COOPACRONYM_AANAME, COOPNAME_AADENOMINATION, .keep_all = TRUE) %>% 
  # The Exporter column might lead to misinterpretation if it's kept in the table as some cooperatives sell to both Cargill and BC, but only one is shown in exporter clean
  dplyr::select(-EXPORTER_CLEAN) 


# Plot coops in and out of ZDCs 
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(coop_exporter_BC_C_coord_sf) +
  tm_dots(col = "COOP_IN_ZDC", title = "Cooperatives part of ZDC \n (i.e. supplying Barry C. and/or Cargill)", size = 0.1, alpha = 0.7, palette = c("coral1", "navy")) +
  tm_scale_bar() +
  tm_layout(legend.outside = T, frame = FALSE)

Merge indicators per cooperative into a single table

# MERGE INDICATORS PER COOPERATIVE ----------------------------

indicators_coop <- coop_vol_exporter_coord_sf  %>%
  dplyr::select(-YEAR) %>% 
  left_join(st_set_geometry(transmute(coop_exporter_BC_C_coord_sf, 
                                      COOPACRONYM_AANAME, COOPNAME_AADENOMINATION,
                                      COOP_IN_ZDC), NULL),
              by = c("COOPACRONYM_AANAME", "COOPNAME_AADENOMINATION")) %>% 
  left_join(st_set_geometry(coop_close_to_pas, NULL), 
              by = c("COOPACRONYM_AANAME", "COOPNAME_AADENOMINATION")) %>%
  left_join(st_set_geometry(coop_in_cfi, NULL), 
            by = c("COOPACRONYM_AANAME", "COOPNAME_AADENOMINATION")) %>% 
  left_join(cdd, by = c("GEOCODE" = "LVL_4_CODE")) %>%
  left_join(territorial_deforestation, by = c("GEOCODE")) %>% 
  dplyr::select(-EXPORTER_CLEAN) #Might lead to misinterpretation if it's kept in the table as some cooperatives sell to other exporters, but only one is shown in exporter clean.

indicators_coop$COOP_IN_CFI_LANDSCAPE[is.na(indicators_coop$COOP_IN_CFI_LANDSCAPE)] <- "0"
indicators_coop$COOP_IN_BUFFER_PA[is.na(indicators_coop$COOP_IN_BUFFER_PA)] <- "0"

Selection of cooperatives based on hypotheses

# ALL HYPOTHESES
indicators_coop_hypotheses <- indicators_coop %>% 
  mutate(HYPOTHESES = case_when(
    COOP_IN_CFI_LANDSCAPE == "0" & COOP_IN_BUFFER_PA == "0" & COOP_IN_ZDC == "0" ~ "CONTROL",
    COOP_IN_CFI_LANDSCAPE == "1" & COOP_IN_BUFFER_PA == "0" & COOP_IN_ZDC == "0" ~ "H1",
    COOP_IN_CFI_LANDSCAPE == "0" & COOP_IN_BUFFER_PA == "1" & COOP_IN_ZDC == "0" ~ "H2",
    COOP_IN_CFI_LANDSCAPE == "0" & COOP_IN_BUFFER_PA == "0" & COOP_IN_ZDC == "1" ~ "H3",
    COOP_IN_CFI_LANDSCAPE == "1" & COOP_IN_BUFFER_PA == "0" & COOP_IN_ZDC == "1" ~ "H4",
    COOP_IN_CFI_LANDSCAPE == "1" & COOP_IN_BUFFER_PA == "1" & COOP_IN_ZDC == "1" ~ "H5",
    TRUE ~ as.character(NA)
  )
    )

ncontrol <- indicators_coop_hypotheses %>% 
  filter(HYPOTHESES == "CONTROL") %>% 
  nrow()

nH1 <- indicators_coop_hypotheses %>% 
  filter(HYPOTHESES == "H1") %>% 
  nrow()

nH2 <- indicators_coop_hypotheses %>% 
  filter(HYPOTHESES == "H2") %>% 
  nrow()

nH3 <- indicators_coop_hypotheses %>% 
  filter(HYPOTHESES == "H3") %>% 
  nrow()

nH4 <- indicators_coop_hypotheses %>% 
  filter(HYPOTHESES == "H4") %>% 
  nrow()

nH5 <- indicators_coop_hypotheses %>% 
  filter(HYPOTHESES == "H5") %>% 
  nrow()



# Plot coops according to hypotheses 
tm_shape(civ_departments) + 
  tm_polygons(col = "gray95") + 
  tmap_options(check.and.fix = TRUE) +
tm_shape(civ_protected_areas) + 
  tm_polygons(border.col = NULL, col = "forestgreen") +
tm_shape(indicators_coop_hypotheses) +
  tm_dots(col = "HYPOTHESES", title = "Cooperatives matching \n selection hypotheses", size = 0.1, palette = c("blue", "red", "purple", "orange", "green", "yellow")) +
  tm_scale_bar() +
  tm_layout(legend.outside = T, frame = FALSE)

This gives us:

18 coops 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)
8 coops for H1, effect of JAs: NO PAs, NOT ZDC
193 coops for H2, effect of PAs: NO JAs, NOT ZDC
18 coops for H3, effect of ZDCs: NO JAs, NO PAs, ARE in ZDC areas
2 coops for H4, synergistic effect of ZDCs + JAs, NO PAs, ARE in areas with JAs and ZDCs
128 coops for H5, synergistic effect of ZDCs + JAs + PAs, ARE in areas with JAs, PAs and ZDCs

For certification, the information is in the table and can be filtered (columns: UTZ; RA; FAIRTRADE; COCOALIFE; SUSTAINABLE_ORIGINS; CERTIFIED), as well as information on deforestation.

Save output file in S3

s3write_using(indicators_coop_hypotheses,
  object = "cote_divoire/cocoa/indicators/out/indicators_coop_HHsurvey_hypotheses.csv",
  bucket = "trase-storage",
  FUN = write_delim, delim = ";",
  opts = c("check_region" = T)
)

Upload file from S3

indicators_coop_HHsurvey_hypotheses <- s3read_using(
  object = "cote_divoire/cocoa/indicators/out/indicators_coop_HHsurvey_hypotheses.csv",
  bucket = "trase-storage",
  FUN = read_delim, 
  delim = ";",
  opts = c("check_region" = T)
)

Sustain-cocoa : household survey site selection (per cooperative)

Cecile Renier

17/08/2021