Grid Sample

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’ fumnction 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\RtmpQJwpjm/rwanda_psu.shp' using driver `ESRI Shapefile'
## Writing 168 features with 3 fields and geometry type Point.

Additional draft coding needs to be brushed up

psu_polygons2<-psu_polygons

plot(psu_polygons)

psu_polygons$urban<-NULL

psu_polygons$RWA__2010_<-NULL

plot(psu_polygons)

psu_polygons$geometry<-NULL

psu_polygons$ADM1_PCODE<-paste0("RW", psu_polygons$strata)

read the shapefile of Rwanda

I have downloaded this shape file from website

mapshape <- st_read("shapefile/rwa_adm3_2006_NISR_WGS1984_20181002.shp")

## Reading layer `rwa_adm3_2006_NISR_WGS1984_20181002' from data source 
##   `C:\Nasif\UVA\Research\other\test\shapefile\rwa_adm3_2006_NISR_WGS1984_20181002.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 416 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 28.86173 ymin: -2.84023 xmax: 30.89975 ymax: -1.047167
## Geodetic CRS:  WGS 84

mapdata <- left_join(mapshape, psu_polygons, by="ADM1_PCODE")

Final plot

ggplot(mapdata) +
  geom_sf(fill = "gray", color = "white") +
  geom_point(data = psu_polygons2, aes(x = st_coordinates(geometry)[, 1], 
                                       y = st_coordinates(geometry)[, 2]), 
             color = "blue")+
  xlab("Longitude")+ylab("Latitude")+
  ggtitle("DHS survey household")+
  ggspatial::annotation_north_arrow(location = "tr", 
                                    pad_x = unit(5, "in"), 
                                  pad_y = unit(0.4, "in"),
            style = ggspatial::north_arrow_nautical(fill = c("grey40", "white"),line_col = "grey20",text_family = "ArcherPro Book"))+
  annotation_scale(location = "br", width_hint = 0.5) +
  coord_sf(crs = 4326)+
  theme(plot.title = element_text(hjust=0.5))+
  theme(legend.position="right")

It looks similar to the article’s plot.