Malaria risk assessment for Burundi

Part 1: Data access

For our malaria risk assessment, we will use data from the Malaria Atlas Project and Worldpop.

Administrative boundaries of our area of interest (AOI)

# Downloads the geometric shape of our AOI on admin 2 level from the Malaria Atlas Project. 
AOI_shp <- malariaAtlas::getShp(ISO = "BDI", admin_level = "admin2")  
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\f.lamy\AppData\Local\Temp\RtmpE9lM8t\shp\shp2580a416db4\mapadmin_2_2018.shp", layer: "mapadmin_2_2018"
## with 47 features
## It has 16 fields
# Creates a simple plot of the data that we have just downloaded
plot(AOI_shp)

Malaria incidence

# Downloads malaria indicator data for our AOI from Malaria Atlas Project
AOI_Pf2020 <- getRaster(surface = "Plasmodium  falciparum Incidence version 2020", shp = AOI_shp)
# Makes a basic plot of the downloaded malaria indicator data.
plot(AOI_Pf2020)

Accessibility to healthcare services

# Downloads data for estimated walking time/minutes to the nearest healthcare facility for our AOI from Amalaria Atlas project
AOI_access <- getRaster(surface = "Walking-only travel time to healthcare map without access to motorized transport", shp = AOI_shp)
# Makes a basic plot of the downloaded accessibility indicator data
plot(AOI_access)

# Population statistics
This dataset was produced based on the 2020 population census/projection-based estimates for 2020.More Info.

# Downloads the population data (covariate) and saves it to disk
pop_path <- wpgpGetCountryDataset(ISO3 = "BDI",
                             covariate = "ppp_2020",
                             destDir = tempdir())
# Uses the "raster" function to load the downloaded data from disk into R
pop <- raster(pop_path)
# The downloaded data has a spatial resolution of 100m. For our analysis we do not need such fine data, it only slows our analysis down.
# Therefore, we aggregate the data by the factor 10 to get a spatial resolution of 1km. 
AOI_pop <- aggregate(pop, fact = 10, fun = sum)
# Clears the "pop" data from memory. We don't need it anymore as we continue working with the aggregated data.
remove(pop)
# Makes a basic plot of the population indicator data
plot(AOI_pop)

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot. # Part 2: Data manipulation and analysis # Turn the raw data into indicators

#----------------
# Analysis - Count Pop, Incidence, Avg. dist to healthcare per Admin area
# Convert administrative borders to a sf object
AOI_sf <- st_as_sf(AOI_shp)
# Sums the population per district
pop_sum <- terra::extract(AOI_pop, AOI_sf, fun = sum, ID = TRUE, na.rm = TRUE) %>% 
  round(digits =2)
pop_district <- cbind(pop_sum,AOI_sf)
# Mean plasmodium falciparum indicence per district
Pf20_mean <- terra::extract(AOI_Pf2020, AOI_sf, fun = mean, ID = TRUE, na.rm = TRUE) %>% 
  round(digits = 2)
Pf20_district <- cbind(Pf20_mean, AOI_sf)
# Mean distance to the nearest healthcare facility per district
access_mean <- terra::extract(AOI_access, AOI_sf, fun = mean, ID = TRUE, na.rm = TRUE) %>% 
  round(digits = 2)
access_district <- cbind(access_mean, AOI_sf)
# Bind all the indicators to the same geospatial outlines
AOI_district <- cbind(pop_district, Pf20_district, access_district)
# Removes the Lake Tanganyika district because it covers only lake areas.
AOI_district <- AOI_district[!(AOI_district$name_2=="Lake Tanganyika"),]
remove(pop_sum, Pf20_mean, access_mean)

Create composite indicator

# Perform min-max normalization on population indicator
AOI_district$pop_norm <- (AOI_district$pop_sum - min(AOI_district$pop_sum))/(max(AOI_district$pop_sum) - min(AOI_district$pop_sum))
# Perform min-max normalization on malaria incidence indicator
AOI_district$pf_norm <- (AOI_district$Pf20_mean - min(AOI_district$Pf20_mean))/(max(AOI_district$Pf20_mean) - min(AOI_district$Pf20_mean))
# Perform min-max normalization on accessibility indicator
AOI_district$access_norm <- (AOI_district$access_mean - min(AOI_district$access_mean))/(max(AOI_district$access_mean) - min(AOI_district$access_mean))
# Sum all three sub-indicators 
AOI_district$composite_indicator <- (AOI_district$pop_norm + AOI_district$pf_norm + AOI_district$access_norm)/3  
AOI_district$composite_indicator <- round(AOI_district$composite_indicator, digits = 2)

Part 3: Visualizing results as maps

tmap_mode("view")
pal <- c('#60b564','#fed400')
pop <- tm_shape(AOI_district) +
  tm_polygons(col = "pop_sum", palette = pal, popup.vars = c("Pop"="pop_sum")) +
  tm_borders(lwd = 0.2) + 
  tm_layout(title = "Total Population") + 
  tm_fill(col = '#BFE1F4') + 
  tm_borders(lwd = 1, col = "black") 
access <- tm_shape(AOI_district) +
  tm_polygons(col = "access_mean", palette = pal, id = "name_2", popup.vars = c("Access"="access_mean")) +
  tm_borders(lwd = 0.2) + 
  tm_layout(title = "Healthcare accessibility") + 
  tm_fill(col = '#BFE1F4') + 
  tm_borders(lwd = 1, col = "black") 
pf20 <- tm_shape(AOI_district) +
  tm_polygons(col = "Pf20_mean", palette = pal, id = "name_2", popup.vars = c("Malaria incidence"="Pf20_mean")) +
  tm_borders(lwd = 0.2) + 
  tm_layout(title = "Malaria incidence") +
  tm_fill(col = '#BFE1F4') + 
  tm_borders(lwd = 1, col = "black")
# Visualize composite indicator as a map
risk <- tm_shape(AOI_district) +
  tm_polygons(col = "composite_indicator", palette = pal, id = "composite_indicator", popup.vars = c("Name" = "name_2", "Risk index"="composite_indicator", "Population" = "pop_norm", "Accessibility" = "access_norm", "Malaria indcidence" = "pf_norm")) +
  tm_borders(lwd = 0.2) + 
  tm_layout(title = "Composite risk indicator") + 
  tm_fill(col = '#BFE1F4') + 
  tm_borders(lwd = 1, col = "black") 
x <- list(pop, access, pf20, risk)
tmap_arrange(x, sync = TRUE)