library(stringr)

# Get the current page number from the file name
current_page <- as.numeric(str_extract(knitr::current_input(), "\\d+"))

# Set the total number of pages
total_pages <- 9

# Generate the URLs for the previous and next pages
previous_page <- ifelse(current_page > 1, paste0("visual_", current_page - 1, "-darfur_violence-code_included.html"), NA)
next_page <- ifelse(current_page < total_pages, paste0("visual_", current_page + 1, "-darfur_violence-code_included.html"), NA)


# Note: All code, insights shared, and methodological notes in this document and other documents in this portfolio ARE NOT open source and MUST be cited.

# For this page, I create 3 sets of faceted pairs of choropleth maps using the tmap package, just like I did for the last visual (Visual 7). Just as in Visual 7, second-level administrative units are colored based on the average number of UN troops (derived from figures in the Geo-PKO dataset) stationed there in the given year for that map and recently closed down UN bases (also taken from the Geo-PKO dataset) are plotted with blue dots. The ACLED dataset is used to plot red bubbles depicting the location of attacks on civilians and fatality numbers per attack during the given year. 

# This set of faceted map pairs differs from Visual 7 only in that UN bases withdrawn from South Darfur in 2017, 2018, and 2019 are plotted for the 3 pairs of maps comparing the year prior to UN base closures with the year of UN base closures, and 50 mile radius buffers are shown emanating from these closed South Darfur bases. The bases closed down in South Darfur during these 3 years can be verified from the same sources listed in Visual 7 (also shown below). Some bases listed in the Geo-PKO dataset as active in a particular year had already been withdrawn in a prior year according to these sources; as a result, I made some edits to the Geo-PKO dataset to correct these inaccuracies. Once again, the buffers plotted on the map are used to calculate and compare the number of attacks on civilians and number of fatalities from such attacks, but for the locations in South Darfur where UN bases were withdrawn rather than the West Darfur base withdrawals.

# Please refer to the detailed comments on the Visual 7 page for a deeper explanation of why I use the code for Visual 8 the way I do below.

# See the following sources for UN base closure dates:

# https://oios.un.org/file/8866/download?token=F6tc9ZI5
# https://www.ecoi.net/en/file/local/2052225/S_2021_470_E.pdf

library(rgdal)
library(tidyverse)
library(tmap)
library(sp)
library(sf)

setwd("C:/Users/rsb84/Desktop/RB/COLUMBIA/QMSS/COURSES/Spring_2021/Data Visualization/End_project")


darfur_internal<-raster::shapefile("C:/Users/rsb84/Desktop/RB/COLUMBIA/QMSS/COURSES/Spring_2021/Data Visualization/End_project/darfur_internal.shp")

# Convert to a simple feature spatial object
darfur_internal_sf <- sf::st_as_sf(darfur_internal)

original_crs <- st_crs(darfur_internal) #Original CRS is EPSG:4326, i.e., the World Geodetic System 1984 (WGS84), a geographic CRS
CRS.new <- st_crs("EPSG:20135") #This is a projected CRS well-suited for Darfur

darfur_internal_sf <- st_transform(darfur_internal_sf, crs = CRS.new) # Transform the CRS to EPSG:20135

ACLED_data <- readxl::read_excel("ACLED-DARFUR-VAC-2008-2021 (After Course Ended-for 2023 Portfolio).xlsx", 
     col_types = c("date", "numeric", "text", 
         "text", "text", "text", "text", "text", 
         "text", "text", "text", "text", "text", 
         "numeric", "numeric", "text", 
         "text", "numeric", "numeric"))


acled_short=ACLED_data[c("event_date","year","sub_event_type","actor1","inter1","admin1","admin2","location","latitude","longitude","attacks","fatalities")]


acled_2016 = subset(acled_short, year == 2016)
acled_2017 = subset(acled_short, year == 2017)
acled_2018 = subset(acled_short, year == 2018)
acled_2019 = subset(acled_short, year == 2019)
acled_2020 = subset(acled_short, year == 2020)
acled_2021 = subset(acled_short, year == 2021)


# Make objects special feature objects and set CRS to the original CRS of the darfur map shapefile - 
acled_2016_sf <- st_as_sf(acled_2016, coords = c("longitude", "latitude"), crs = original_crs)
acled_2017_sf <- st_as_sf(acled_2017, coords = c("longitude", "latitude"), crs = original_crs)
acled_2018_sf <- st_as_sf(acled_2018, coords = c("longitude", "latitude"), crs = original_crs)
acled_2019_sf <- st_as_sf(acled_2019, coords = c("longitude", "latitude"), crs = original_crs)
acled_2020_sf <- st_as_sf(acled_2020, coords = c("longitude", "latitude"), crs = original_crs)
acled_2021_sf <- st_as_sf(acled_2021, coords = c("longitude", "latitude"), crs = original_crs)

# Transform the CRS to EPSG:20135
acled_2016_sf <- st_transform(acled_2016_sf, crs=CRS.new)
acled_2017_sf <- st_transform(acled_2017_sf, crs=CRS.new)
acled_2018_sf <- st_transform(acled_2018_sf, crs=CRS.new)
acled_2019_sf <- st_transform(acled_2019_sf, crs=CRS.new)
acled_2020_sf <- st_transform(acled_2020_sf, crs=CRS.new)
acled_2021_sf <- st_transform(acled_2021_sf, crs=CRS.new)

geo_pko <- readxl::read_excel("Geo-PKO-v-2-1-2023_portfolio-verified_bases_only-LATEST.xlsx")


# 2016

geo_pko_for_maps_16 = filter(geo_pko, year == 2016)
geo_pko_df_16 = geo_pko_for_maps_16[c("year","location","latitude", "longitude")]
geo_pko_16_avg.troops_per_base = geo_pko_for_maps_16 %>% group_by(location) %>% summarise(avg.troops = mean(no.troops, na.rm=TRUE))
geo_pko_df_16 = inner_join(geo_pko_16_avg.troops_per_base, geo_pko_df_16, by="location")
geo_pko_df_16 = distinct(geo_pko_df_16)
active_bases_2016_points = geo_pko_df_16

# 2017

geo_pko_for_maps_17 = filter(geo_pko, year == 2017)
geo_pko_df_17 = geo_pko_for_maps_17[c("year","location","latitude", "longitude")]
geo_pko_17_avg.troops_per_base = geo_pko_for_maps_17 %>% group_by(location) %>% summarise(avg.troops = mean(no.troops, na.rm=TRUE))
geo_pko_df_17 = inner_join(geo_pko_17_avg.troops_per_base, geo_pko_df_17, by="location")
geo_pko_df_17 = distinct(geo_pko_df_17)
active_bases_2017_points = geo_pko_df_17

# 2018

geo_pko_for_maps_18 = filter(geo_pko, year == 2018)
geo_pko_df_18 = geo_pko_for_maps_18[c("year","location","latitude", "longitude")]
geo_pko_18_avg.troops_per_base = geo_pko_for_maps_18 %>% group_by(location) %>% summarise(avg.troops = mean(no.troops, na.rm=TRUE))
geo_pko_df_18 = inner_join(geo_pko_18_avg.troops_per_base, geo_pko_df_18, by="location")
geo_pko_df_18 = distinct(geo_pko_df_18)
active_bases_2018_points = geo_pko_df_18

# 2019

geo_pko_for_maps_19 = filter(geo_pko, year == 2019)
geo_pko_df_19 = geo_pko_for_maps_19[c("year","location","latitude", "longitude")]
geo_pko_19_avg.troops_per_base = geo_pko_for_maps_19 %>% group_by(location) %>% summarise(avg.troops = mean(no.troops, na.rm=TRUE))
geo_pko_df_19 = inner_join(geo_pko_19_avg.troops_per_base, geo_pko_df_19, by="location")
geo_pko_df_19 = distinct(geo_pko_df_19)
active_bases_2019_points = geo_pko_df_19


# Convert active UN bases for each given year to simple feature spatial objects
active_bases_2016_points_sf <- sf::st_as_sf(active_bases_2016_points, coords = c("longitude", "latitude"), crs = original_crs)

active_bases_2017_points_sf <- sf::st_as_sf(active_bases_2017_points, coords = c("longitude", "latitude"), crs = original_crs)

active_bases_2018_points_sf <- sf::st_as_sf(active_bases_2018_points, coords = c("longitude", "latitude"), crs = original_crs)

active_bases_2019_points_sf <- sf::st_as_sf(active_bases_2019_points, coords = c("longitude", "latitude"), crs = original_crs)


# Transform the CRS to EPSG:20135
active_bases_2016_points_sf <- st_transform(active_bases_2016_points_sf, crs=CRS.new)
active_bases_2017_points_sf <- st_transform(active_bases_2017_points_sf, crs=CRS.new)
active_bases_2018_points_sf <- st_transform(active_bases_2018_points_sf, crs=CRS.new)
active_bases_2019_points_sf <- st_transform(active_bases_2019_points_sf, crs=CRS.new)



# It is critical to replace all NA values for the variable `avg.troops` to 0
active_bases_2016_points_sf$avg.troops[is.na(active_bases_2016_points_sf$avg.troops)] <- 0
active_bases_2017_points_sf$avg.troops[is.na(active_bases_2017_points_sf$avg.troops)] <- 0
active_bases_2018_points_sf$avg.troops[is.na(active_bases_2018_points_sf$avg.troops)] <- 0
active_bases_2019_points_sf$avg.troops[is.na(active_bases_2019_points_sf$avg.troops)] <- 0



# Ensure all avg_troop variable values for the active_bases point objects are read as numbers
active_bases_2016_points_sf$avg.troops <- as.numeric(active_bases_2016_points_sf$avg.troops)
active_bases_2017_points_sf$avg.troops <- as.numeric(active_bases_2017_points_sf$avg.troops)
active_bases_2018_points_sf$avg.troops <- as.numeric(active_bases_2018_points_sf$avg.troops)
active_bases_2019_points_sf$avg.troops <- as.numeric(active_bases_2019_points_sf$avg.troops)



# Now, for active_bases_2016_points_sf - active_bases_2020_points_sf, filter these objects such that only the UN bases listed occur where the average number of troops is not 0. This ensures only active UN bases having UN troops present remain, which will then be plotted in the maps further below:

active_bases_2016_points_sf <- active_bases_2016_points_sf[active_bases_2016_points_sf$avg.troops != 0, ]
active_bases_2017_points_sf <- active_bases_2017_points_sf[active_bases_2017_points_sf$avg.troops != 0, ]
active_bases_2018_points_sf <- active_bases_2018_points_sf[active_bases_2018_points_sf$avg.troops != 0, ]
active_bases_2019_points_sf <- active_bases_2019_points_sf[active_bases_2019_points_sf$avg.troops != 0, ]


# At this point, I make "contains" spatial joins with R. See my comments in the code of Visual 6 for more understanding of the following code:
acled_2016_joined_sf <- sf::st_join(acled_2016_sf, join = st_within, left=FALSE, darfur_internal_sf)
acled_2017_joined_sf <- sf::st_join(acled_2017_sf, join = st_within, left=FALSE, darfur_internal_sf)
acled_2018_joined_sf <- sf::st_join(acled_2018_sf, join = st_within, left=FALSE, darfur_internal_sf)
acled_2019_joined_sf <- sf::st_join(acled_2019_sf, join = st_within, left=FALSE, darfur_internal_sf)
acled_2020_joined_sf <- sf::st_join(acled_2020_sf, join = st_within, left=FALSE, darfur_internal_sf)
acled_2021_joined_sf <- sf::st_join(acled_2021_sf, join = st_within, left=FALSE, darfur_internal_sf)



acled_2016_joined_sf$fatalities <- as.numeric(acled_2016_joined_sf$fatalities)
acled_2017_joined_sf$fatalities <- as.numeric(acled_2017_joined_sf$fatalities)
acled_2018_joined_sf$fatalities <- as.numeric(acled_2018_joined_sf$fatalities)
acled_2019_joined_sf$fatalities <- as.numeric(acled_2019_joined_sf$fatalities)
acled_2020_joined_sf$fatalities <- as.numeric(acled_2020_joined_sf$fatalities)
acled_2021_joined_sf$fatalities <- as.numeric(acled_2021_joined_sf$fatalities)


# Spatial joining data from the active_bases point objects to the darfur_internal shapefile. Refer back to my comments in the identical code in Visual 6 if you have any questions about the code:
darfur_map_contains_2016_bases <- sf::st_join(darfur_internal_sf, join = st_contains, left=TRUE, active_bases_2016_points_sf)
darfur_map_contains_2017_bases <- sf::st_join(darfur_internal_sf, join = st_contains, left=TRUE, active_bases_2017_points_sf)
darfur_map_contains_2018_bases <- sf::st_join(darfur_internal_sf, join = st_contains, left=TRUE, active_bases_2018_points_sf)
darfur_map_contains_2019_bases <- sf::st_join(darfur_internal_sf, join = st_contains, left=TRUE, active_bases_2019_points_sf)


# It is critical to replace all NA values for the variable `avg.troops` with 0 so that the choropleth map can be colored based on values:
darfur_map_contains_2016_bases$avg.troops[is.na(darfur_map_contains_2016_bases$avg.troops)] <- 0
darfur_map_contains_2017_bases$avg.troops[is.na(darfur_map_contains_2017_bases$avg.troops)] <- 0
darfur_map_contains_2018_bases$avg.troops[is.na(darfur_map_contains_2018_bases$avg.troops)] <- 0
darfur_map_contains_2019_bases$avg.troops[is.na(darfur_map_contains_2019_bases$avg.troops)] <- 0


# Next, for the variable `year`, any observations that have the value NA instead of 2016 for the 2016 map or 2020 for the 2020 map will be changed to the correct year.
darfur_map_contains_2016_bases$year[is.na(darfur_map_contains_2016_bases$year)] <- 2016
darfur_map_contains_2017_bases$year[is.na(darfur_map_contains_2017_bases$year)] <- 2017
darfur_map_contains_2018_bases$year[is.na(darfur_map_contains_2018_bases$year)] <- 2018
darfur_map_contains_2019_bases$year[is.na(darfur_map_contains_2019_bases$year)] <- 2019


darfur_map_contains_2016_bases$avg.troops <- as.numeric(darfur_map_contains_2016_bases$avg.troops)
darfur_map_contains_2017_bases$avg.troops <- as.numeric(darfur_map_contains_2017_bases$avg.troops)
darfur_map_contains_2018_bases$avg.troops <- as.numeric(darfur_map_contains_2018_bases$avg.troops)
darfur_map_contains_2019_bases$avg.troops <- as.numeric(darfur_map_contains_2019_bases$avg.troops)


setwd("C:/Users/rsb84/Desktop/RB/COLUMBIA/QMSS/COURSES/Spring_2021/Data Visualization/End_project")

# Create a bounding box in WGS84
my_box <- rgeos::bbox2SP(n = 16.12054, s = 8.64132, w = 21.61368, e = 29.58313, proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

# Reproject the bounding box to EPSG:20135
my_box <- spTransform(my_box, CRS("+init=epsg:20135"))


fifty_miles_in_meters <- 50 * 1609.344

# UN Bases in South Darfur Closed in 2017:
closed_bases_south_2017 <- readxl::read_excel("closed_bases_south_2017-portfolio_2023.xlsx", 
     col_types = c("text", "text", 
         "numeric", "numeric"))
closed_bases_south_2017_sf=st_as_sf(closed_bases_south_2017, coords = c("longitude", "latitude"), crs = original_crs)

# UN Bases in South Darfur Closed in 2018:
closed_bases_south_2018 <- readxl::read_excel("closed_bases_south_2018-portfolio_2023.xlsx", 
     col_types = c("text", "text", 
         "numeric", "numeric"))
closed_bases_south_2018_sf=st_as_sf(closed_bases_south_2018, coords = c("longitude", "latitude"), crs = original_crs)

# UN Bases in South Darfur Closed in 2019:
closed_bases_south_2019 <- readxl::read_excel("closed_bases_south_2019-portfolio_2023.xlsx", 
     col_types = c("text", "text", 
         "numeric", "numeric"))
closed_bases_south_2019_sf=st_as_sf(closed_bases_south_2019, coords = c("longitude", "latitude"), crs = original_crs)

# Transform the following data to be plotted into a Projected CRS (i.e., EPSG:20135) that uses units of measurement of meters rather than degrees

closed_bases_south_2017_sf <- st_transform(closed_bases_south_2017_sf, CRS.new)
closed_bases_south_2018_sf <- st_transform(closed_bases_south_2018_sf, CRS.new)
closed_bases_south_2019_sf <- st_transform(closed_bases_south_2019_sf, CRS.new)


# Create buffer and union buffer objects around UN bases closed in South Darfur in 2017, 2018, and 2019:

# 2017
closed_bases_south_2017_buffers <- st_buffer(x=closed_bases_south_2017_sf, dist =  fifty_miles_in_meters)
union_buffers_2017 <- st_union(closed_bases_south_2017_buffers)

# 2018
closed_bases_south_2018_buffers <- st_buffer(x=closed_bases_south_2018_sf, dist =  fifty_miles_in_meters)
union_buffers_2018 <- st_union(closed_bases_south_2018_buffers)

# 2019
closed_bases_south_2019_buffers <- st_buffer(x=closed_bases_south_2019_sf, dist =  fifty_miles_in_meters)
union_buffers_2019 <- st_union(closed_bases_south_2019_buffers)


troops_breaks <- c(0, 200, 400, 600, 800, 1000)
fatalities_breaks <- c(0, 10, 20, 30, 40, 50)
fatalities_labels <- c("0 to 10", "10 to 20", "20 to 30", "30 to 40", "40 to 50")


# Define thick black borders for first order administrative units
adm1_borders_col = "black"  # Color for first order administrative borders
adm1_borders_lwd = 2       # Line width for first order administrative borders

# Currently, the darfur_map_contains_2016_bases, darfur_map_contains_2017_bases, darfur_map_contains_2018_bases, and darfur_map_contains_2019_bases sf map objects will display the geometries of their second order administrative districts. We want to create a second map of darfur that contains the names and geometries of the first order administrative districts; we can do this using the code below. This map will be used to provide thick black borders of North, South, East, West, and Central Darfur (which are Darfur's first order administrative districts) in the faceted tmap maps below even though each faceted maps will have their second order administrative districts be colored based on the average number of UN troops stationed within those respective second order administrative districts. Since all we are doing with this first order administrative district map is providing thick black boundaries, we can just use one of the second order administrative district maps in the code below

darfur_map_1st_order_adm <- darfur_map_contains_2016_bases %>%
  group_by(ADM1_EN) %>% # column with names of first order administrative districts
  summarize(geometry = st_union(geometry)) %>% #column with geometries of the second order administrative districts
  ungroup()


# You should put tmap_mode('plot') (for static maps) or tmap_mode('view') (for interactive maps) before creating the map with tm_shape() or any other tmap function. This sets the mode for all subsequent tmap functions until the mode is changed.

tmap_mode('plot')


# Layout with Title
layout1 <- tm_layout(
          title = "Attacks on Civilians, Before vs. After the 2017\n Edd el Fursan and Tulus Base Closures",
          title.position = c("left","top"),
          title.size=20, title.color="steelblue", fontface = 7,
          legend.title.size = 1,
          legend.title.color = "steelblue",
          legend.text.size = 0.6,
          legend.text.color = "steelblue",
          #legend.position = c("right","center"),
          legend.position = c(0.75,0.25),#c(horizontal, vertical)
          legend.bg.color = "white",
          legend.bg.alpha = 1,
          bg.color="white",
          frame=FALSE)

layout2 <- tm_layout(
          title = "Attacks on Civilians, Before vs. After the 2018\n Buram and Graida Base Closures",
          title.position = c("left","top"),
          title.size=20, title.color="steelblue", fontface = 7,
          legend.title.size = 1,
          legend.title.color = "steelblue",
          legend.text.size = 0.6,
          legend.text.color = "steelblue",
          #legend.position = c("right","center"),
          legend.position = c(0.75,0.25),#c(horizontal, vertical)
          legend.bg.color = "white",
          legend.bg.alpha = 1,
          bg.color="white",
          frame=FALSE)

layout3 <- tm_layout(
          title = "Attacks on Civilians, Before vs. After the 2019\n Nyala Base Closure",
          title.position = c("left","top"),
          title.size=20, title.color="steelblue", fontface = 7,
          legend.title.size = 1,
          legend.title.color = "steelblue",
          legend.text.size = 0.6,
          legend.text.color = "steelblue",
          #legend.position = c("right","center"),
          legend.position = c(0.76,0.25),#c(horizontal, vertical)
          legend.bg.color = "white",
          legend.bg.alpha = 1,
          bg.color="white",
          frame=FALSE)

# Define the layout object without title
layout_no_title <- tm_layout(
          title = "",  
          title.size=20, title.color="steelblue", fontface = 7,
          legend.title.size = 1,
          legend.title.color = "steelblue",
          legend.text.size = 0.6,
          legend.text.color = "steelblue",
          #legend.position = c("right","center"),
          #legend.position = c(0.79,0.25),#c(horizontal, vertical)
          legend.bg.color = "white",
          legend.bg.alpha = 1,
          bg.color="white",
          frame=FALSE)

# 2016 vs. 2017 - First Set of Maps for Comparison

t_2016.1 <- tm_shape(darfur_map_contains_2016_bases, bbox = my_box) + 
tm_fill("avg.troops", title="Average Number\nof UN Troops", palette="Greens") + 
tm_borders() +
tm_shape(darfur_map_1st_order_adm, bbox = my_box) +
tm_borders(col=adm1_borders_col, lwd=adm1_borders_lwd) +
tm_shape(closed_bases_south_2017_sf) + 
tm_dots(size=0.5, col="steelblue") +
tm_shape(acled_2016_joined_sf) +
tm_bubbles(size="fatalities", scale = 2.7, col="red", alpha=0.3, jitter=0.05, title.size="Fatalities", breaks=c(0,10,20,30,40,50), sizes.legend=c((16/17)*5,(16/17)*10,(16/17)*15,(16/17)*20,(16/17)*25), sizes.legend.labels=c("10","20","30","40","50")) +
tm_shape(union_buffers_2017, projection = CRS.new) +
#tm_fill(col = "blue", alpha = 0.2) +
tm_borders(col = "blue", alpha = 1, lty = 1, lwd = 1.8) +
layout1 +
tm_add_legend(title="") +
tm_credits(text="2016", position=c("center", "bottom"), size=1.5)



t_2017.1 <- tm_shape(darfur_map_contains_2017_bases, bbox = my_box) + 
tm_fill("avg.troops", title="Average Number\nof UN Troops", palette="Greens") + 
tm_borders() +
tm_shape(darfur_map_1st_order_adm, bbox = my_box) +
tm_borders(col=adm1_borders_col, lwd=adm1_borders_lwd) +
tm_shape(closed_bases_south_2017_sf) + 
tm_dots(size=0.5, col="steelblue") +
tm_shape(acled_2017_joined_sf) +
tm_bubbles(size="fatalities", scale=1.25, col="red", alpha=0.3, jitter=0.05, title.size="Fatalities", breaks=c(0,10,20,30,40,50), sizes.legend=c(5,10,15,20,25), sizes.legend.labels=c("10","20","30","40","50")) + 
tm_shape(union_buffers_2017, projection = CRS.new) +
#tm_fill(col = "blue", alpha = 0.2) +
tm_borders(col = "blue", alpha = 1, lty = 1, lwd = 1.8) +
layout_no_title +
tm_add_legend(title="") +
tm_credits(text="2017", position=c("center", "bottom"), size=1.5) +
tm_legend(show = FALSE)

# 2017 vs. 2018 - Second Set of Maps for Comparison

t_2017.2 <- tm_shape(darfur_map_contains_2017_bases, bbox = my_box) + 
tm_fill("avg.troops", title="Average Number\nof UN Troops", palette="Greens") + 
tm_borders() +
tm_shape(darfur_map_1st_order_adm, bbox = my_box) +
tm_borders(col=adm1_borders_col, lwd=adm1_borders_lwd) +
tm_shape(closed_bases_south_2018_sf) + 
tm_dots(size=0.5, col="steelblue") +
tm_shape(acled_2017_joined_sf)+
tm_bubbles(size="fatalities", scale=1.25, col="red", alpha=0.3, jitter=0.05, title.size="Fatalities", breaks=c(0,10,20,30,40,50), sizes.legend=c((37/32)*5,(37/32)*10,(37/32)*15,(37/32)*20,(37/32)*25), sizes.legend.labels=c("10","20","30","40","50")) +   
tm_shape(union_buffers_2018, projection = CRS.new) +
#tm_fill(col = "blue", alpha = 0.2) +
tm_borders(col = "blue", alpha = 1, lty = 1, lwd = 1.8) +
layout2 +
tm_add_legend(title="") +
tm_credits(text="2017", position=c("center", "bottom"), size=1.5)


t_2018.2 <- tm_shape(darfur_map_contains_2018_bases, bbox =my_box) + 
tm_fill("avg.troops", title="Average Number\nof UN Troops", palette="Greens") + 
tm_borders() +
tm_shape(darfur_map_1st_order_adm, bbox = my_box) +
tm_borders(col=adm1_borders_col, lwd=adm1_borders_lwd) +
tm_shape(closed_bases_south_2018_sf) + 
tm_dots(size=0.5, col="steelblue") +
tm_shape(acled_2018_joined_sf) +
tm_bubbles(size="fatalities", scale=2, col="red", alpha=0.3, jitter=0.05, title.size="Fatalities", breaks=c(0,10,20,30,40,50), sizes.legend=c(10,20,30,40,50), sizes.legend.labels=c("10","20","30","40","50")) + 
tm_shape(union_buffers_2018, projection = CRS.new) +  
#tm_fill(col = "blue", alpha = 0.2) +
tm_borders(col = "blue", alpha = 1, lty = 1, lwd = 1.8) +
layout_no_title +
tm_add_legend(title="") +
tm_credits(text="2018", position=c("center", "bottom"), size=1.5) + 
tm_legend(show = FALSE)

# 2018 v. 2019 - Third Set of Maps for Comparison

t_2018.3 <- tm_shape(darfur_map_contains_2018_bases, bbox =my_box) + 
tm_fill("avg.troops", title="Average Number\nof UN Troops", palette="Greens") + 
tm_borders() +
tm_shape(darfur_map_1st_order_adm, bbox = my_box) +
tm_borders(col=adm1_borders_col, lwd=adm1_borders_lwd) +
tm_shape(closed_bases_south_2019_sf) + 
tm_dots(size=0.5, col="steelblue") +
tm_shape(acled_2018_joined_sf)+
tm_bubbles(size="fatalities", scale=2, col="red", alpha=0.3, jitter=0.05, title.size="Fatalities", breaks=c(0,10,20,30,40,50), sizes.legend=c((27/60)*10,(27/60)*20,(27/60)*30,(27/60)*40,(27/60)*50), sizes.legend.labels=c("10","20","30","40","50")) + 
tm_shape(union_buffers_2019, projection = CRS.new) +
#tm_fill(col = "blue", alpha = 0.2) +
tm_borders(col = "blue", alpha = 1, lty = 1, lwd = 1.8) +
layout3 +
tm_add_legend(title="") +
tm_credits(text="2018", position=c("center", "bottom"), size=1.5)


t_2019.3 <- tm_shape(darfur_map_contains_2019_bases, bbox = my_box) + 
tm_fill("avg.troops", title="Average Number\nof UN Troops", palette="Greens") + 
tm_borders() +
tm_shape(darfur_map_1st_order_adm, bbox = my_box) +
tm_borders(col=adm1_borders_col, lwd=adm1_borders_lwd) +
tm_shape(closed_bases_south_2019_sf) + 
tm_dots(size=0.5, col="steelblue") +
tm_shape(acled_2019_joined_sf) +
tm_bubbles(size="fatalities", scale=2, col="red", alpha=0.3, jitter=0.05, title.size="Fatalities", breaks=c(0,10,20,30,40,50), sizes.legend=c(5,10,15,20,25), sizes.legend.labels=c("10","20","30","40","50")) +
tm_shape(union_buffers_2019, projection = CRS.new) +
#tm_fill(col = "blue", alpha = 0.2) +
tm_borders(col = "blue", alpha = 1, lty = 1, lwd = 1.8) +
layout_no_title +
tm_add_legend(title="") +
tm_credits(text="2019", position=c("center", "bottom"), size=1.5) +
tm_legend(show = FALSE)


In each of the visuals below, I have created 50 mi radius buffers around all UN bases that were withdrawn in South Darfur in 2017, 2018, and 2019. After 2019, no more bases were withdrawn from South Darfur. As with West Darfur in Visual 7, I have also calculated the total attacks on civilians and fatalities from such attacks inside these buffers and plotted them in a heat map table further below for analysis.


Whereas in 2016 (the year before the Edd El Fursan and Tulus UN bases were closed) there were 47 attacks on civilians and 58 fatalities within 50 mi of these bases, in 2017 there were 25 attacks on civilians (a 47% decrease) and 31 fatalities (an 87% decrease).

tmap_arrange(t_2016.1, t_2017.1, nrow = 1, ncol = 2, asp=1, outer.margins = 0)

Looking at the 2017 vs. 2018 maps below, whereas in 2017 (the year before the Buram and Graida UN compounds were closed) there were 32 attacks on civilians and 28 fatalities within a 50 mi radius of these bases, in 2018 there were 27 attacks and 28 fatalities - roughly the same level of violence. A caveat should be noted: it appears from the map that although troops were withdrawn from the Buram and Graida UN bases in 2018, there were actually more UN troops within the buffer in 2018 than in 2017. This is likely a result of transfering troops from the closed bases to other nearby bases that were still active.


tmap_arrange(t_2017.2, t_2018.2, nrow = 1, ncol = 2, asp=1, outer.margins = 0)

Further, while there were 21 attacks on civilians and 22 fatalities in 2018 within 50 mi radius of the Nyala UN compound before it was shut down, in 2019 there were only 10 attacks on civilians (a 110% decrease) and 6 fatalities (a 267% decrease).


tmap_arrange(t_2018.3, t_2019.3, nrow = 1, ncol = 2, asp=1, outer.margins = 0)

# CALCULATNG ATTACKS ON CIVILIANS & FATALITIES PER YEAR INSIDE 50 MILE RADIUS BUFFERS AROUND BASES CLOSED IN 2017, 2018, and 2019

# length() calculates the number of civilian attacks
# sum() calculates the number of fatalities from these attacks


# Calculating attacks (length) and fatalities (sum) in 2016 - Buffer 2017
fatality_sum_2016.buffer_2017 = sf::st_filter(acled_2016_joined_sf, union_buffers_2017, .predicate = st_within)
fatality_sum_2016.buffer_2017$fatalities = as.numeric(fatality_sum_2016.buffer_2017$fatalities)
attacks_2016_buffer_2017 <- length(fatality_sum_2016.buffer_2017$fatalities)
deaths_2016_buffer_2017 <- sum(fatality_sum_2016.buffer_2017$fatalities)

# Calculating attacks (length) and fatalities (sum) in 2017 - buffer 2017
fatality_sum_2017.buffer_2017 = sf::st_filter(acled_2017_joined_sf, union_buffers_2017, .predicate = st_within)
attacks_2017_buffer_2017 <- length(fatality_sum_2017.buffer_2017$fatalities)
deaths_2017_buffer_2017 <- sum(fatality_sum_2017.buffer_2017$fatalities)

# Calculating attacks (length) and fatalities (sum) in 2018 - buffer 2017
fatality_sum_2018.buffer_2017 = sf::st_filter(acled_2018_joined_sf, union_buffers_2017, .predicate = st_within)
attacks_2018_buffer_2017 <- length(fatality_sum_2018.buffer_2017$fatalities)
deaths_2018_buffer_2017 <- sum(fatality_sum_2018.buffer_2017$fatalities)

# Calculating attacks (length) and fatalities (sum) in 2019 - buffer 2017
fatality_sum_2019.buffer_2017 = sf::st_filter(acled_2019_joined_sf, union_buffers_2017, .predicate = st_within)
attacks_2019_buffer_2017 <- length(fatality_sum_2019.buffer_2017$fatalities)
deaths_2019_buffer_2017 <- sum(fatality_sum_2019.buffer_2017$fatalities)

# Calculating attacks (length) and fatalities (sum) in 2020 - buffer 2017
fatality_sum_2020.buffer_2017 = sf::st_filter(acled_2020_joined_sf, union_buffers_2017, .predicate = st_within)
attacks_2020_buffer_2017 <- length(fatality_sum_2020.buffer_2017$fatalities)
deaths_2020_buffer_2017 <- sum(fatality_sum_2020.buffer_2017$fatalities)

# Calculating attacks (length) and fatalities (sum) in 2021 - buffer 2017
fatality_sum_2021.buffer_2017 = sf::st_filter(acled_2021_joined_sf, union_buffers_2017, .predicate = st_within)
fatality_sum_2021.buffer_2017$fatalities = as.numeric(fatality_sum_2021.buffer_2017$fatalities)
attacks_2021_buffer_2017 <- length(fatality_sum_2021.buffer_2017$fatalities)
deaths_2021_buffer_2017 <- sum(fatality_sum_2021.buffer_2017$fatalities)


# Calculating attacks (length) and fatalities (sum) in 2016 - Buffer in 2018
fatality_sum_2016.buffer_2018 = sf::st_filter(acled_2016_joined_sf, union_buffers_2018, .predicate = st_within)
attacks_2016_buffer_2018 <- length(fatality_sum_2016.buffer_2018$fatalities)
deaths_2016_buffer_2018 <- sum(fatality_sum_2016.buffer_2018$fatalities)

# Calculating attacks (length) and fatalities (sum) in 2017 - Buffer in 2018
fatality_sum_2017.buffer_2018 = sf::st_filter(acled_2017_joined_sf, union_buffers_2018, .predicate = st_within)
attacks_2017_buffer_2018 <- length(fatality_sum_2017.buffer_2018$fatalities)
deaths_2017_buffer_2018 <- sum(fatality_sum_2017.buffer_2018$fatalities)

# Calculating attacks (length) and fatalities (sum) in 2018 - Buffer in 2018
fatality_sum_2018.buffer_2018 = sf::st_filter(acled_2018_joined_sf, union_buffers_2018, .predicate = st_within)
attacks_2018_buffer_2018 <- length(fatality_sum_2018.buffer_2018$fatalities)
deaths_2018_buffer_2018 <- sum(fatality_sum_2018.buffer_2018$fatalities)

# Calculating attacks (length) and fatalities (sum) in 2019 - Buffer in 2018
fatality_sum_2019.buffer_2018 = sf::st_filter(acled_2019_joined_sf, union_buffers_2018, .predicate = st_within)
attacks_2019_buffer_2018 <- length(fatality_sum_2019.buffer_2018$fatalities)
deaths_2019_buffer_2018 <- sum(fatality_sum_2019.buffer_2018$fatalities)

# Calculating attacks (length) and fatalities (sum) in 2020 - Buffer in 2018
fatality_sum_2020.buffer_2018 = sf::st_filter(acled_2020_joined_sf, union_buffers_2018, .predicate = st_within)
attacks_2020_buffer_2018 <- length(fatality_sum_2020.buffer_2018$fatalities)
deaths_2020_buffer_2018 <- sum(fatality_sum_2020.buffer_2018$fatalities)

# Calculating attacks (length) and fatalities (sum) in 2021 - Buffer in 2018
fatality_sum_2021.buffer_2018 = sf::st_filter(acled_2021_joined_sf, union_buffers_2018, .predicate = st_within)
fatality_sum_2021.buffer_2018$fatalities=as.numeric(fatality_sum_2021.buffer_2018$fatalities)
attacks_2021_buffer_2018 <- length(fatality_sum_2021.buffer_2018$fatalities)
deaths_2021_buffer_2018 <- sum(fatality_sum_2021.buffer_2018$fatalities)


# Calculating attacks (length) and fatalities (sum) in 2016 - Buffer 2019
fatality_sum_2016.buffer_2019 = sf::st_filter(acled_2016_joined_sf, union_buffers_2019, .predicate = st_within)
attacks_2016_buffer_2019 <- length(fatality_sum_2016.buffer_2019$fatalities)
deaths_2016_buffer_2019 <- sum(fatality_sum_2016.buffer_2019$fatalities)

# Calculating attacks (length) and fatalities (sum) in 2017 - Buffer 2019
fatality_sum_2017.buffer_2019 = sf::st_filter(acled_2017_joined_sf, union_buffers_2019, .predicate = st_within)
attacks_2017_buffer_2019 <- length(fatality_sum_2017.buffer_2019$fatalities)
deaths_2017_buffer_2019 <- sum(fatality_sum_2017.buffer_2019$fatalities)

# Calculating attacks (length) and fatalities (sum) in 2018 - Buffer 2019
fatality_sum_2018.buffer_2019 = sf::st_filter(acled_2018_joined_sf, union_buffers_2019, .predicate = st_within)
attacks_2018_buffer_2019 <- length(fatality_sum_2018.buffer_2019$fatalities)
deaths_2018_buffer_2019 <- sum(fatality_sum_2018.buffer_2019$fatalities)

# Calculating attacks (length) and fatalities (sum) in 2019 - Buffer 2019
fatality_sum_2019.buffer_2019 = sf::st_filter(acled_2019_joined_sf, union_buffers_2019, .predicate = st_within)
attacks_2019_buffer_2019 <- length(fatality_sum_2019.buffer_2019$fatalities)
deaths_2019_buffer_2019 <- sum(fatality_sum_2019.buffer_2019$fatalities)

# Calculating attacks (length) and fatalities (sum) in 2020 - Buffer 2019
fatality_sum_2020.buffer_2019 = sf::st_filter(acled_2020_joined_sf, union_buffers_2019, .predicate = st_within)
attacks_2020_buffer_2019 <- length(fatality_sum_2020.buffer_2019$fatalities)
deaths_2020_buffer_2019 <- sum(fatality_sum_2020.buffer_2019$fatalities)

# Calculating attacks (length) and fatalities (sum) in 2021 - Buffer 2019
fatality_sum_2021.buffer_2019 = sf::st_filter(acled_2021_joined_sf, union_buffers_2019, .predicate = st_within)
attacks_2021_buffer_2019 <- length(fatality_sum_2021.buffer_2019$fatalities)
deaths_2021_buffer_2019 <- sum(fatality_sum_2021.buffer_2019$fatalities)


We can see that in South Darfur these figures do not lend support to the theory that the withdrawal of UN compounds led to increasing patterns of violence against civilians in the immediate term. However, a similar trend as was seen in Visual 7 in West Darfur repeats itself in South Darfur with escalating violence years later and higher than usual fatality per attack rates. Looking at the table below, when attacks on civilians and fatalities are calculated in future years one can see that although violence was already prevalent in the area inside these buffers between the years 2016 and 2018, in the years 2020 and 2021 there was a noticeable increase in civilian killings beyond usual levels. Further, the relatively high levels of violence within all 3 buffers between the years 2016 and 2018 arguably should have raised reservations among UN officials about withdrawing at the pace they did.


In 2020, within 50 mi of the bases that were closed in 2017, there were just 12 attacks on civilians but there were 57 fatalities (an 84% increase compared to 2017), which equates to 4.75 fatalities per attack. Also in 2020, within 50 mi of the bases closed in 2018, there were only 16 attacks but there were 78 fatalities (a 179% increase compared to 2018), which equates to 4.88 fatalities per attack. Additionally, within 50 mi of the UN compound that closed in 2019, while there were just 10 attacks in 2020, but there were 32 fatalities (a 417% increase compared to 2019), which equates to 3.10 fatalities per attack. In 2021 in South Darfur, locations surrounding the UN base withdrawals in 2017 and 2018 had a slight reduction in fatalities per attack, with 1.83 fatalities per attack and 1.73 fatalities per attack, respectively - still higher than usual levels. This patterns of higher than normal fatality per attack ratios is similar to the experience of West Darfur in 2019 and 2021.


However, just as with the spikes in West Darfur violence in 2019 and 2021 following UN base withdrawals, the power vacuum from years of UN base withdrawals is likely only one of several factors which led to the increase in violence in South Darfur. Such factors which have been cited include intercommunal violence in retaliation for tribal cattle raiding and land disputes between displaced land owners and those occupying their land.


library(kableExtra)

# Create a dataframe
df <- data.frame(
  Attacks1 = c(attacks_2016_buffer_2017, attacks_2017_buffer_2017, attacks_2018_buffer_2017, attacks_2019_buffer_2017, attacks_2020_buffer_2017, attacks_2021_buffer_2017),  
  Fatalities1 = c(deaths_2016_buffer_2017, deaths_2017_buffer_2017, deaths_2018_buffer_2017, deaths_2019_buffer_2017, deaths_2020_buffer_2017, deaths_2021_buffer_2017),
  Attacks2 = c(attacks_2016_buffer_2018, attacks_2017_buffer_2018, attacks_2018_buffer_2018, attacks_2019_buffer_2018, attacks_2020_buffer_2018, attacks_2021_buffer_2018),
  Fatalities2 = c(deaths_2016_buffer_2018, deaths_2017_buffer_2018, deaths_2018_buffer_2018, deaths_2019_buffer_2018, deaths_2020_buffer_2018, deaths_2021_buffer_2018),
  Attacks3 = c(attacks_2016_buffer_2019, attacks_2017_buffer_2019, attacks_2018_buffer_2019, attacks_2019_buffer_2019, attacks_2020_buffer_2019, attacks_2021_buffer_2019),
  Fatalities3 = c(deaths_2016_buffer_2019, deaths_2017_buffer_2019, deaths_2018_buffer_2019, deaths_2019_buffer_2019, deaths_2020_buffer_2019, deaths_2021_buffer_2019),
  row.names = c("2016", "2017", "2018", "2019", "2020", "2021")
)

# Since dataframes in R cannot have duplicate column names, this code below changes the column names just for the display to allow for duplicate column names under "2017 Closed Base Areas", "2018 Closed Base Areas", and "2019 Closed Base Areas"

df_display <- df
names(df_display) <- rep(c("Attacks", "Fatalities"), 3)


# Below is the code for a heat map table in which the degree of redness in each cell depends on how large or small the integer values are in comparison to all of the other integers in the table. The value of having this type of table is that it quickly helps the reader understand trends by seeing which cells, columns, and/or rows are most important to look at. The font color is white for cell values > the 3rd quantile, and black otherwise. Through trial and error, I found that using such a high threshold helped the integer values inside the cells with both the darker as well as lighter shades of red to be clearly readable. A lower threshold will result in difficulty reading some integer values. Visual 7 used a threshold of the 8th decile for its heat map table, suggesting that tuning between the 7th and 8th deciles may be necessary to obtain the clearest table in most cases.

library(scales)

# Get min, max, and 3rd quantile values for the entire dataframe for scaling
min_val <- min(df)
max_val <- max(df)
quantile_val <- quantile(unlist(df), 0.75) # calculate the 3rd quantile from the original numeric dataframe

df_display[] <- lapply(df_display, function(x) {
  if(is.numeric(x)) {
    cell_spec(x, "html",
              color = ifelse(x > quantile_val, "white", "black"),
              background = scales::gradient_n_pal(c("white", "red"))(scales::rescale(x, from = c(min_val, max_val))))
  } else {
    x
  }
})


# Code for the table
kable_output <- kbl(df_display, format = "html", escape = FALSE, align = 'c') %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F) %>%
  add_header_above(c(" " = 1, "Area around 2017 Closed Bases" = 2,
                     "Area around 2018 Closed Bases" = 2,
                     "Area around 2019 Closed Bases" = 2))

# Center each column name, column contents, and each header
kable_output <- kable_output %>%
  column_spec(1, bold = TRUE) %>%
  column_spec(2, extra_css = "text-align: center;") %>%
  column_spec(3, extra_css = "text-align: center;", border_right = TRUE) %>%
  column_spec(4, extra_css = "text-align: center;") %>%
  column_spec(5, extra_css = "text-align: center;", border_right = TRUE) %>%
  column_spec(6, extra_css = "text-align: center;") %>%
  column_spec(7, extra_css = "text-align: center;", border_right = TRUE) %>%
  row_spec(0, extra_css = "width:150px;") # Ensures the headers are centered by having sufficient space to separate them
kable_output
Area around 2017 Closed Bases
Area around 2018 Closed Bases
Area around 2019 Closed Bases
Attacks Fatalities Attacks Fatalities Attacks Fatalities
2016 47 58 47 55 48 27
2017 25 31 32 28 35 22
2018 32 37 27 28 21 22
2019 11 16 13 16 10 6
2020 12 57 16 78 10 31
2021 18 33 26 45 18 17