# Iterate through shapefile and output an image for each feature using sf in R

I would like to create some statistical profiles for a range of features; specifically, 2016 Census profiles for the Census Divisions (CD) within British Columbia (BC), Canada. This would require a map showing each CD highlighted with respect to the other CDs (28 in total).

To do this using the sf package in R and the Statistics Canada Boundary files for 2016 Census Divisions (ZIP format) [~32MB] I would follow these pseudo-code steps:

2. Filter on BC (PRUID = 59)
3. Select the first feature/CD
4. Plot the entire set of features (black line with white fill)
5. Plot (with add = TRUE) the selected CD with some colour (black line with coloured fill)
6. Move to/select next feature/CD and go to Step 4

How would I code this in R?

The idea would be to generalize this so one can cycle through any set of features within a shapefile and output individual features highlighted within the entire geography.

The steps to perform this task using sf would be:

# Filter on British Columbia/BC (PRUID == 59) (Step 2)
bc_cd_2016 <- can_cd_2016 %>% filter(PRUID == '59')

# Set general plot parameters (if needed):
# MARgins to 0 lines c(bottom, left, top, right)
# Line WiDth to 0.1
par(mar = c(0,0,0,0), lwd = 0.1)

# Loop through each CD (Steps 3 and 6)
for (CD in 1:(nrow(bc_cd_2016))) {
# Open PNG stream to 'bc_cd_2016_CD#.png'
png(paste('bc_cd_2016_', bc_cd_2016\$CDUID[CD], '.png', sep = ''))
plot(bc_cd_2016)# Plot BC (Step 4)
plot(bc_cd_2016[CD,], col = 'red', add = TRUE)# Plot CD/feature (Step 5)
dev.off()# Close PNG stream
}

The output is a set of files with title bc_cd_2016_#.png where # denotes the CDUID. Example of bc_cd_2016_5947.png:

...or when cycling through bc_cd_2016_59??.png (using ImageMagick's convert to GIF):