5

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:

  1. Load the source .shp
  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.

Related: How do I generate an image of each polygon in a shapefile?

1 Answer 1

5

The steps to perform this task using sf would be:

library('sf') # Load sf package
library('dplyr') # Load dplyr package
# Load 2016 Census shape file for Canada (Step 1)
can_cd_2016 <- st_read('lcd_000b16a_e.shp')
# 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:

enter image description here

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

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.