Load the relevant packages

library(raster)

Read the population raster (image) data

population_raster <- raster("RWA_ppp_2010_adj_v2.tif")

Plot of the population raster

plot(population_raster)

Read the shape file

A population raster and a shapefile RWAshp defining the various strata in the census. Here, the strata shapefile consists of the 30 districts.

load("RWAshp.rda")

Rasterize by strata and create rasterize plot for strata

strata_raster <- rasterize(RWAshp,population_raster,field = "STL.2")
plot(strata_raster)

Urban VS Rural

To classify urban and rural areas, a raster will be created based on population density. Urban areas will be defined by selecting the densest cells.

total_pop <- cellStats(population_raster, stat = "sum")
urban_pop_value <- total_pop * .16
pop_df <- data.frame(index = 1:length(population_raster[]),
                     pop = population_raster[])
pop_df <- pop_df[!is.na(pop_df$pop), ]
pop_df <- pop_df[order(pop_df$pop,decreasing = T), ]
pop_df$cumulative_pop <- cumsum(pop_df$pop)
pop_df$urban <- 0
pop_df$urban[which(pop_df$cumulative_pop <= urban_pop_value)] <- 1
urban_raster <- population_raster >= min(subset(pop_df,urban == 1)$pop)
plot(urban_raster)

Required additional packages

library(sf)
library(tidyverse)
library(ggplot2)
library(ggspatial)

load the ‘sample_psu’ function, similar function of gs_sample from gridsample package

load("sample_psu.RData")

sample calculation using gridsample algorithm with arbitrary value

psu_polygons <- sample_psu(
  population_raster = population_raster,
  strata_raster = strata_raster,
  urban_raster = urban_raster,
  cfg_min_pop_per_cell = 0.01,
  cfg_hh_per_stratum = 416,
  cfg_hh_per_urban = 26,
  cfg_hh_per_rural = 26,
  cfg_max_psu_size = 10,
  cfg_pop_per_psu = 610,
  cfg_sample_rururb = TRUE,
  cfg_sample_spatial = FALSE,
  cfg_sample_spatial_scale = 100,
  output_path = tempdir(),
  sample_name = "rwanda_psu"
)
## Writing layer `rwanda_psu' to data source 
##   `C:\Users\nasif\AppData\Local\Temp\RtmpySZ5S1/rwanda_psu.shp' using driver `ESRI Shapefile'
## Writing 168 features with 3 fields and geometry type Point.

Grid sample cells with PPES

Reassign population data with a new name

pop_data <- population_raster

Convert the raster to polygons

pop_polygons <- rasterToPolygons(pop_data, dissolve = TRUE)
pop_sf <- st_as_sf(pop_polygons)

Visualization of primary sampling units (PSUs)

ggplot() +
  geom_sf(data = pop_sf, aes(fill = RWA_ppp_2010_adj_v2), color = "white", size = 0.1) +
  geom_point(data = psu_polygons, aes(x = st_coordinates(geometry)[, 1], 
                                       y = st_coordinates(geometry)[, 2]), 
             color = "blue")+
  scale_fill_gradient(name = "Population Density", low = "mistyrose", high = "brown") +
  labs(title = "Primary sampling units of Rwanda",
       x = "Longitude", y = "Latitude") +
  theme(plot.title = element_text(hjust=0.5))+
  theme(legend.position="right")+
  annotation_scale(location = "bl") +
  annotation_north_arrow(location = "tl", which_north = "true", style = north_arrow_orienteering())

It looks similar to the article’s Figure 3 (top).

Downscaling

Define the extent to zoom into

y_min <- -2.00  # Example longitude minimum
y_max <- -1.90  # Example longitude maximum
x_min <- 30.01   # Example latitude minimum
x_max <- 30.15   # Example latitude maximum

Filter population polygons within this extent

zoomed_area <- st_crop(pop_sf, xmin = x_min, xmax = x_max, ymin = y_min, ymax = y_max)
zoomed_ss <- st_crop(psu_polygons, xmin = x_min, xmax = x_max, ymin = y_min, ymax = y_max)

Final downscaling plot for Kigali city, Rwanda

ggplot() +
  geom_sf(data = zoomed_area, aes(fill = RWA_ppp_2010_adj_v2), color = "white", size = 0.1) +
  geom_point(data = zoomed_ss, aes(x = st_coordinates(geometry)[, 1], 
                                   y = st_coordinates(geometry)[, 2]), 
             color = "blue", size = 3)+
  scale_fill_gradient(name = "Population Density", low = "mistyrose", high = "brown") +
  coord_sf(xlim = c(x_min, x_max), ylim = c(y_min, y_max)) +
  labs(title = "Zoomed-In Population Grid (Kigali Rwanda)",
       x = "Longitude", y = "Latitude") +
  theme(plot.title = element_text(hjust=0.5))+
  theme(legend.position="right")+
  annotation_scale(location = "bl") +
  annotation_north_arrow(location = "tl", which_north = "true", style = north_arrow_orienteering())

It looks similar to the article’s Figure 5.