In dark-peak-analytics/who_decide_AP: R Package to facilitate easy access to geospatial model data

knitr::opts_chunk$set(
  eval = FALSE,
  collapse = TRUE,
  comment = "#>"
)

The aim of this Vignette is to demonstrate how the data retrieved using the bolsteR package can be integrated with other R packages into a more advanced analysis pipeline.

Note: It is assumed that you are already familiar with the bolsteR package and its main functions. See the 'Simple Application Programming Interface (API) querys' vignette for a brief intorduction to the bolsteR package

Sample analysis

The objective of this sample analysis will be to determine what proportion of the population in Kenya lives within a 5km radius to a hospital, which is taken as a proxy for access to health care.

First, we load the bolsteR package, and two other packages from the WorldPop Open Population Repository (WOPR), which provides additional functions to query and process geo-spatial information.

knitr::opts_chunk$set(warning = FALSE, message = FALSE)
# devtools::install_github("wpgp/wopr")
# devtools::install_github("wpgp/wpCPR")
library(bolsteR)
library(wopr)
library(wpCPR)
library(GADMTools) 

Retrieve regional boundaries for Kenya:

ken_gadm <- gadm_sp_loadCountries(
  fileNames = 'KEN', 
  level = 2, 
  basefile = './shapefiles', 
  simplify=0.025
  )$spdf
# adjust projection
proj4string(ken_gadm) <- CRS("+proj=longlat +datum=WGS84")
# map regional boundaries
plot(ken_gadm)

Subset the data:

For simplicity, the analysis will focus on just one particular region: Baringo Central

Baringo = subset(ken_gadm, ken_gadm$NAME_2 == "Baringo Central")
Baringo_map = st_as_sf(Baringo) # transform into sf
plot(ken_gadm)
plot(Baringo_map$geometry, add = T, col = "red")

As a first step, we retrieve a population estimate for the Baringo region.

file_path = paste0(tempdir(), "/", "test2.shp")
writeOGR(
  obj=Baringo, 
  dsn=file_path, 
  layer="", 
  driver="ESRI Shapefile"
  )
Baringo_pop = wpCPRPopulation(
  year = 2020,
  shapeFilePath = file_path,
  verbose = T
  )
Baringo_pop[,c("NAME_0","NAME_1","NAME_2", "pop")]

Next, we use functions from the bolsteR package to retrieve the locations of health facilities in Kenya.

# set api key for healthsites.io
api_key <- set_hs_API_key("123xyz")
Kenya <- Country_R6$new("Kenya", api_key)
# api query to fetch the locations of health facilities in Kenya
Kenya$initiate_HS_API()

We now want to filter the facility data to obtain only those facilities that are hospitals located in Baringo. For this, we first select all hospitals, then determine the location of each hospital, and finally subset the data set based on the location.

# select hospitals
hospitals_kenya = Kenya$health_facilities_data[
  Kenya$health_facilities_data$attributes.amenity == "hospital" &
    !is.na(Kenya$health_facilities_data$attributes.amenity),
  ]

# create sf object
pnts_sf <- st_as_sf(hospitals_kenya,coords = c('lng','lat'))
st_crs(pnts_sf) <- CRS("+proj=longlat +datum=WGS84") # set projection

# determine the region for each health facility in the data set
facilities_by_region = st_intersects(
  st_geometry(pnts_sf), 
  st_as_sf(ken_gadm), 
  sparse = F
  )
colnames(facilities_by_region) <- ken_gadm$NAME_2

# select facilities in Baringo
facilities_in_Baringo <- subset(
  pnts_sf, 
  facilities_by_region[,colnames(facilities_by_region) == "Baringo Central"]
  )

# create a 5km radius around each facility
facilities_in_Baringo_circles <- st_buffer(facilities_in_Baringo, dist = 5000)

# visualise results: shown are the locations of the hospitals in Baringo (red) and the 5km buffer around them
plot(Baringo)
plot(st_geometry(facilities_in_Baringo), col="red", add = T)
plot(st_geometry(facilities_in_Baringo_circles),  add = T)

Now we can query the API again to retrieve population estimates for each hospital catchment area (i.e. 5km buffer).

# 
facilities_in_Baringo_circles = facilities_in_Baringo_circles[,c("geometry","attributes.name","attributes.amenity")]
names(facilities_in_Baringo_circles) = c("geometry","name","type")
write_sf(facilities_in_Baringo_circles, file_path, append = F)
facility_Baringo_pop = wpCPRPopulation(
  year = 2020,
  shapeFilePath = file_path,
  verbose = T
)
facility_Baringo_pop

To derive the total population within all catchment areas combined, we compute the union of all buffers and query the API again.

facility_Baringo_union = sf::st_union(facilities_in_Baringo_circles)
facility_Baringo_union <- st_cast(facility_Baringo_union, "POLYGON")
write_sf(facility_Baringo_union, file_path, append = F)
facility_union_Baringo_pop = wpCPRPopulation(
  year = 2020,
  shapeFilePath = file_path,
  verbose = T
)
facility_union_Baringo_pop

Finally, we can compute the percentage of the population that lives within 5km of any of the hospitals as with

round( sum(facility_union_Baringo_pop$pop) / Baringo_pop$pop, 4)*100

dark-peak-analytics/who_decide_AP documentation built on May 25, 2022, 8:31 p.m.