Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
R script using Simple Features library to identify 1st and 2nd degree neighbors; plot results using Leaflet
# 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))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.