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)


In creating a before vs. after UN base withdrawal comparison of violence against civilians, since Geo-PKO data on UNAMID troop numbers only extends to 2020, I will compare Darfur’s 2016 and 2020 maps. It is important though to keep in mind from Visual 2 that civilian killings were higher in Darfur as a whole in 2021 than 2020, and from Visual 3 that a significant rise in civilian killings occurred in West Darfur in 2021 compared to 2020.


In the maps below, between 2016 and 2020 the pattern of civilian attacks in Darfur shifted notably. In 2016, the deadliest attacks were highly concentrated at the intersection of Central, South, and North Darfur. However, by 2020 there were more high-fatality attacks in South Darfur, while West Darfur saw a significant rise in generally lower-fatality attacks coinciding with the withdrawal of UN bases and peacekeepers.


Both maps reveal a consistent pattern: administrative districts with the most civilian attacks had UN bases averaging a yearly troop strength of 200 or fewer.


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

# In this plot, I create faceted choropleth maps using the tmap library for the years 2016 and 2020 in which Darfur administrative units are colored by the average number of UN troops from nearby bases in a given year. The geographic locations of active bases housing UN troops are plotted as nodes. At the same time, red bubbles are plotted indicating the locations of violent attacks on civilians, with larger size bubbles indicating larger numbers of civilian fatalities from a given attack.

library(tidyverse)
library(dplyr)
library(tmap)
library(sp)
options(sf_use_s2 = FALSE)
library(sf)
library(readxl)

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


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_2020 = subset(acled_short, year == 2020)

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

# Sources for the dates of UN base closures:
# https://oios.un.org/file/8866/download?token=F6tc9ZI5
# https://www.ecoi.net/en/file/local/2052225/S_2021_470_E.pdf


#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 


#2020

geo_pko_for_maps_20 = filter(geo_pko, year == 2020)
geo_pko_df_20 = geo_pko_for_maps_20[c("year","location","latitude", "longitude")]
geo_pko_20_avg.troops_per_base = geo_pko_for_maps_20 %>% group_by(location) %>% summarise(avg.troops = mean(no.troops, na.rm=TRUE))
geo_pko_df_20 = inner_join(geo_pko_20_avg.troops_per_base, geo_pko_df_20, by="location")
geo_pko_df_20 = distinct(geo_pko_df_20)
active_bases_2020_points = geo_pko_df_20 


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


# Converting dataframes to spatial objects, namely "simple feature objects":


# Convert data frames to simple feature objects
active_bases_2016_points <- sf::st_as_sf(active_bases_2016_points, coords = c("longitude", "latitude"), crs = 4326)

# When converting to simple feature objects data from the `longitude` and `latitude` columns are extracted and combined into a new `geometry` column. But for plotting later, it will be best to keep the `longitude` and `latitude` columns. I can do this by creating new columns (`longitude` and `latitude`) whose values are extracted from `geometry`. I do this below for the `active_bases_2016` simple feature object:


active_bases_2020_points <- sf::st_as_sf(active_bases_2020_points, coords = c("longitude", "latitude"), crs = 4326)

# Ensuring darfur_internal shape file uses the same CRS as the  active_bases point objects:

# First, convert 'SpatialPolygonsDataFrame' to an 'sf' object
darfur_internal <- st_as_sf(darfur_internal)

# Now you can transform the CRS to be the same as what active_bases_2016_points uses
darfur_internal <- st_transform(darfur_internal, st_crs(active_bases_2016_points))


# Ensure there are no NA values for the variable avg.troops for administrative region-based observations which have no UN bases with average troops numbers recorded. Any NAs should be converted to 0s.


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


my_box <- rgeos::bbox2SP(n = 16.12054, s = 8.64132, w = 24.61368, e = 27.58313, proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

# Preparing the ACLED dataset point objects to have data from the darfur_internal map joined to them:

# Convert ACLED data frames to simple feature objects (spatial point objects), and ensure crs = 4326 (a very common crs value) is applied to the ACLED point data. We know that this CRS has already been applied to the darfur_internal polygon map:

acled_2016 <- sf::st_as_sf(acled_2016, coords = c("longitude", "latitude"), crs = 4326)
acled_2020 <- sf::st_as_sf(acled_2020, coords = c("longitude", "latitude"), crs = 4326)


# Spatial joining data from darfur_internal to the ACLED point objects:

# .predicate is the type of spatial join. E.g., in the case below, acled_2016 and acled_2020 are point objects, and since I want to keep point objects as the end result these are placed first. darfur_internal is a polygon shape file I want to add data/info from to the point objects, so it is the second object listed. For .predicate, I choose "st_within" because the point objects are within the polygon object. left=FALSE is chosen because I do not want to do a "left join" because I only want to keep the ACLED point objects that fall within the darfur_internal map - not any falling outside the boundaries of the map:

acled_2016 <- st_join(acled_2016, .predicate="st_within", left=FALSE, darfur_internal)
acled_2020 <- st_join(acled_2020, .predicate="st_within", left=FALSE, darfur_internal)

#ACLED, 2016
acled_2016 <- st_as_sf(acled_2016, coords = c("longitude", "latitude"), crs = 4326) 

#ACLED, 2020
acled_2020 <- st_as_sf(acled_2020, coords = c("longitude", "latitude"), crs = 4326)


# Spatial joining data from active_bases_2016_points and active_bases_2020_points to the darfur_internal shapefile:

# The darfur_internal map shapefile contains second level administrative unit polygons within it. After conducting the spatial joins, the second level administrative units contain as attributes the then currently existing active UN bases. The second level administrative units will be colored by the average number of UN troops present from all bases within them during the given year.

# Recall that the .predicate is the type of spatial join. E.g., in the case below, darfur_internal is a polygon shape file, and since I want to keep a polygon-based map of darfur as the end result, this is placed first. active_bases_2016_points and active_bases_2020_points are point objects that I want to add data/info from to the darfur_internal map, so they are the second objects listed. For .predicate, I choose "st_contains" because the darfur_internal polygon object "contains" the point objects within it. left=TRUE is chosen because I want to keep all administrative boundaries of the darfur_internal map - even those administrative boundaries which contain no point objects from active_bases_2016_points and active_bases_2020_points within them. If I were to instead set left=FALSE, then the Darfur map's administrative unit polygons which do not contain any active bases within them would be discarded.


darfur_map_contains_2016_bases <- st_join(darfur_internal, .predicate="st_contains", left=TRUE, active_bases_2016_points)

darfur_map_contains_2020_bases <- st_join(darfur_internal, .predicate="st_contains", left=TRUE, active_bases_2020_points)


# 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_2020_bases$avg.troops[is.na(darfur_map_contains_2020_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_2020_bases$year[is.na(darfur_map_contains_2020_bases$year)] <- 2020


CRS.new <- st_crs("EPSG:4326")
darfur_map_contains_2016_bases <- st_as_sf(darfur_map_contains_2016_bases)
darfur_map_contains_2016_bases <- st_transform(darfur_map_contains_2016_bases, CRS.new)
darfur_map_contains_2020_bases <- st_as_sf(darfur_map_contains_2020_bases)
darfur_map_contains_2020_bases <- st_transform(darfur_map_contains_2020_bases, CRS.new)


acled_2016$fatalities <- as.numeric(acled_2016$fatalities)
acled_2020$fatalities <- as.numeric(acled_2020$fatalities)

darfur_map_contains_2016_bases$avg.troops <- as.numeric(darfur_map_contains_2016_bases$avg.troops)
darfur_map_contains_2020_bases$avg.troops <- as.numeric(darfur_map_contains_2020_bases$avg.troops)
active_bases_2016_points$avg.troops <- as.numeric(active_bases_2016_points$avg.troops)
active_bases_2020_points$avg.troops <- as.numeric(active_bases_2020_points$avg.troops)

# Now, for active_bases_2016 and active_bases_2020, filter these objects such that only the UN bases are listed 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 below:

active_bases_2016_points <- active_bases_2016_points[active_bases_2016_points$avg.troops != 0, ]
active_bases_2020_points <- active_bases_2020_points[active_bases_2020_points$avg.troops != 0, ]


tmap_mode("plot")


# Set the breaks for "Average Number of UN Troops"
troops_breaks <- c(0, 200, 400, 600, 800, 1000)
# Set the breaks and labels for "Fatalities"
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()

# Define the layout object without title
layout <- tm_layout(
          title = "Fatalities from Attacks on Darfur Civilians: 2016 vs. 2020\nBlue Dots are Active UN Bases",
          title.position = c("left","top"),
          title.size=20, title.color="steelblue", fontface = 7,
          legend.title.size = 1.45,
          #legend.title.size = 1,
          legend.title.color = "steelblue",
          legend.text.size = 1.156,
          #legend.text.size = 0.8,
          legend.text.color = "steelblue",
          legend.position = c(0.75,0.25),#legend.position = c("right","bottom"),
          legend.bg.color = "white",
          legend.bg.alpha = 1,
          bg.color="white",
          frame=FALSE,
          inner.margins = c(0.05, 0.05, 0.05, 0.00),
          outer.margins = c(2, 0.1, 1, 0.9),
          between.margin = 0,
          asp=1.5)

# Define the layout object without title
layout2 <- tm_layout(
          title = "",  
          title.size=20, title.color="steelblue", fontface = 7,
          legend.title.size = 1.45,
          #legend.title.size = 1,
          legend.title.color = "steelblue",
          legend.text.size = 1.156,
          #legend.text.size = 0.8,
          legend.text.color = "steelblue",
          legend.position = c(0.75,0.25),#legend.position = c("right","bottom"),
          legend.bg.color = "white",
          legend.bg.alpha = 1,
          bg.color="white",
          frame=FALSE,
          inner.margins = c(0.05, 0.05, 0.05, 0.00),
          outer.margins = c(2, 0.1, 1, 0.9),
          between.margin = 1,
          asp=1.5)


# Modify the tm_bubbles functions for both maps
t=tm_shape(darfur_map_contains_2016_bases, bbox =my_box) + tm_fill("avg.troops", title="Average Number\nof UN Troops", palette="Greens", breaks=troops_breaks) + 
layout +
tm_borders() +
tm_shape(darfur_map_1st_order_adm, bbox = my_box) +
tm_borders(col=adm1_borders_col, lwd=adm1_borders_lwd) +
tm_shape(active_bases_2016_points) + 
tm_symbols(size=0.5, col="steelblue", shapes.legend = "Active UN Bases", shapes=22) +
#tm_dots(col = "steelblue", size = 0.5, shape = NA) +
tm_shape(acled_2016)+
tm_bubbles(size="fatalities", col="red", alpha=0.4, jitter=0.05, title.size="Fatalities", breaks=c(0,20,40,60,80,100), scale=2.7, sizes.legend=c(0.95*6.25,0.95*12.5,0.95*18.75,0.95*25,0.95*31.25), sizes.legend.labels=c("10","20","30","40","50"))+
tm_credits(text="2016", position=c("center", "bottom"), size=1.5) + tm_layout(
  #legend.position = c("right", "center")
  legend.position = c(0.75,0.25)
  )


t2=tm_shape(darfur_map_contains_2020_bases, bbox =my_box) + tm_fill("avg.troops", title="Average Number\nof UN Troops", palette="Greens", breaks=troops_breaks) + 
layout2 +
tm_borders() +
tm_shape(darfur_map_1st_order_adm, bbox = my_box) +
tm_borders(col=adm1_borders_col, lwd=adm1_borders_lwd) +
tm_shape(active_bases_2020_points) + 
tm_symbols(size=0.5, col="steelblue", shapes.legend = "Active UN Bases", shapes=22) +
#tm_dots(col = "steelblue", size = 0.5, shape = NA) +
tm_shape(acled_2020)+
tm_bubbles(size="fatalities", col="red", alpha=0.4, jitter=0.05, title.size="Fatalities", breaks=c(0,20,40,60,80,100), scale=2, sizes.legend=c(20,40,60,80,100), sizes.legend.labels=c("10","20","30","40","50"))+
tm_credits(text="2020", position=c("center", "bottom"), size=1.5) + tm_legend(show = FALSE)

# Define the width of the maps. Adjust the width as needed
map_width <- unit(10, "cm")


t <- t + tm_layout(
  inner.margins = c(0.055, 0.385, .045, 0.260),   # Increase this value to zoom out more
  outer.margins = c(1.7, 1, 1, 0),  # Increase the first value to add space on the left
  between.margin = 0,
  asp = 1 
)


t2 <- t2 + tm_layout(
  inner.margins = c(0.055, 0.38, .045, 0.260),   # Increase this value to zoom out more
  outer.margins = c(1.7, 1, 1, 0), # Increase the first value to add space on the left
  between.margin = 0,
  asp = 1
)

# Combine the two maps using tmap_arrange()
grid_map <- tmap_arrange(t, t2, ncol = 2, nrow = 1, widths = c(1.1*map_width, 1.1*map_width))

# Display the grid map
grid_map

From Visual 3, recall that in 2016 the 3 states with the most civilian killings were Central Darfur (253), South Darfur (139), and North Darfur (114), totaling 506, with West Darfur only having 36. Yet, by 2020, the 3 states with the most killings were South Darfur (95), Central Darfur (35), and West Darfur (31), amounting to 161 fatalities between them - a significant overall decrease from 2016. However, in 2021 West Darfur’s civilian killings rose to 155, the highest for any region and a 400% increase from the previous year.


We can see from the maps above that although fatalities in Central and South Darfur decreased from 2016 to 2020, West Darfur maintained similar fatalities in both years but experienced a sharp increase in 2021, coinciding with UNAMID’s complete closure.


Also, recall that Visual 5’s 2016 heat maps show a high concentration of civilian attacks at the junction of Central, South, and North Darfur, regardless of fatality count. By 2020, the highest density of attacks shifted to north-eastern Central Darfur, with a moderate increase in West Darfur.


Visual 5’s 2021 Darfur heat map indicated a slight increase in the concentration of civilian attacks in West Darfur compared to 2020. This initially seems unexpected given Visual 3’s data showing West Darfur had the highest number of civilian killings of any state in 2021. However, this likely reflects a high fatality rate per attack, despite a less dense concentration of incidents. Similarly, the 2020 heatmap shows a mild concentration of attacks in South Darfur, contrasting with this visual’s indication that the most fatal attacks occurred there.


The analysis so far lends itself to the view that it is plausible the withdrawal of UN bases in these 2 regions may have been premature and allowed violence to spread.


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