R script using Simple Features library to identify 1st and 2nd degree neighbors; plot results using Leaflet

Find_nth_Degree_Neighbors.r

# Filename: Find_nth_Degree_Neighbors.r
# Author: Anthony Louis D'Agostino 
# Date Created: 01/18/2017
# Last Edited: 
# Purpose: Tutorial example of identifying n-th degree neighbors using North Carolina county boundaries 
# Notes: Borrows from Edzer Pebesma's sf Vignette #5 (https://cran.r-project.org/web/packages/sf/vignettes/sf5.html)

# -- preamble; installs necessary packages 

remove(list = ls(all.names = TRUE))


# -- nice preamble courtesy of Timo Grossenbacher (https://timogrossenbacher.ch/2016/12/beautiful-thematic-maps-with-ggplot2-only/#general-ggplot2-theme-for-map)
detachAllPackages <- function() {
  basic.packages.blank <-  c("stats", 
                             "graphics", 
                             "grDevices", 
                             "utils", 
                             "datasets", 
                             "methods", 
                             "base")
  basic.packages <- paste("package:", basic.packages.blank, sep = "")
  
  package.list <- search()[ifelse(unlist(gregexpr("package:", search())) == 1, 
                                  TRUE, 
                                  FALSE)]
  
  package.list <- setdiff(package.list, basic.packages)
  
  if (length(package.list) > 0)  for (package in package.list) {
    detach(package, character.only = TRUE)
    print(paste("package ", package, " detached", sep = ""))
  }
}

detachAllPackages()

if (!require(Matrix)) {
  install.packages("Matrix", repos = "http://cran.us.r-project.org")
  require(Matrix)
}
if (!require(sf)) {
  install.packages("sf", repos = "http://cran.us.r-project.org")
  require(sf)
}
if (!require(ggplot2)) {
  install.packages("ggplot2", repos = "http://cran.us.r-project.org")
  require(ggplot2)
}
if (!require(leaflet)) {
  install.packages("leaflet", repos = "http://cran.us.r-project.org")
  require(leaflet)
}



FIRSTdegreeNeighbors <- function(x) { 
  
  # -- use sf functionality and output to sparse matrix format for minimizing footprint 
  first.neighbor <- st_touches(x, x, sparse = TRUE) 
  
  # -- convert results to data.frame (via sparse matrix)
  n.ids <- sapply(first.neighbor, length)
  vals <- unlist(first.neighbor) 
  out <- sparseMatrix(vals, rep(seq_along(n.ids), n.ids))
  
  out.summ <- summary(out)  # -- this is currently only generating row values, need to map to actual obs [next line]
  data.frame(county = x[out.summ$j,]$NAME, 
                       countyid = x[out.summ$j,]$CNTY_ID, 
                       firstdegreeneighbors = x[out.summ$i,]$NAME, 
                       firstdegreeneighborid = x[out.summ$i,]$CNTY_ID, 
                       stringsAsFactors = FALSE)
}


# -- read in NC shapefile 
nc <- st_read(system.file("shape/nc.shp", package="sf"))

# -- data frame of all first degree neighbors 
nc.neighbors <- FIRSTdegreeNeighbors(nc)

# -- let's focus on a single county 
county.of.interest <- 'Chatham'

# -- subset to first degree neighbors of specified county   
county.1st <- subset(nc.neighbors, county == county.of.interest)

# -- rename column to avoid problems with 2nd iteration
colnames(county.1st)[3] <- "firstiterationfirstdegree"
county.2nd <- merge(county.1st, nc.neighbors, by.x = "firstiterationfirstdegree", by.y = "county") 

# -- find second degree neighbors that are not in first degree neighbor set 
county.2nd.only <- county.2nd[!county.2nd$firstdegreeneighbors %in% county.1st$firstiterationfirstdegree,]
county.2nd.only <- subset(county.2nd.only, !firstdegreeneighbors %in% county.of.interest) 

# -- plot verifying we've captured all the adjacent neighbors of N degree
par(mar=c(0,0,0,0))
plot(st_geometry(nc), col = c(grey(0.9)), axes = TRUE)

# -- add county of interest 
plot(st_geometry(subset(nc, NAME %in% county.of.interest)), col = 'purple', add = TRUE)

# -- add 1st degree neighbors 
plot(st_geometry(subset(nc, NAME %in% county.1st$firstiterationfirstdegree)), col = 'red', add = TRUE) 

# -- add 2nd degree neighbors 
plot(st_geometry(subset(nc, NAME %in% county.2nd.only$firstdegreeneighbors)), col = 'blue', add = TRUE) 

# -- alternatively, produce a nice leaflet plot 
               
chatham.county <- st_geometry(subset(nc, NAME %in% county.of.interest)) %>%
  st_transform(crs = 4326) %>%
  as("Spatial")

first.degree <- subset(nc, NAME %in% county.1st$firstiterationfirstdegree) %>% 
  st_transform(4236) %>% 
    as("Spatial")

second.degree <- subset(nc, NAME %in% county.2nd.only$firstdegreeneighbors) %>%
  st_transform(4236) %>% 
    as("Spatial")


chatham.base <- leaflet(chatham.county) %>% 
  addTiles() %>% 
  addPolygons(fillColor = "#e41a1c",
              opacity = 0.4, fillOpacity = 0.5,
              weight = 1, 
              smoothFactor = 0.2, 
              highlightOptions = highlightOptions(color = "white", weight = 0.2,
                                                  bringToFront = TRUE))

chatham.base %>% 
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(data = first.degree, stroke = TRUE, 
              opacity = 0.4, fillOpacity = 0.5,
              weight = 1, 
              smoothFactor = 0.2, 
              fillColor = "#fc8d62",
              highlightOptions = highlightOptions(color = "black", weight = 0.2,
                                                  bringToFront = TRUE)) %>%
  addPolygons(data = second.degree, stroke = TRUE,
              opacity = 0.3, fillOpacity = 0.5,
              weight = 1, 
              smoothFactor = 0.2, 
              fillColor = "#8da0cb",
              highlightOptions = highlightOptions(color = "white", weight = 1,
                                                  bringToFront = TRUE))