Mountain View

DRAFT :: Calculating Neighbouring Polygons in R

· by Mikey Harper · Read in about 3 min · (500 words) ·
GIS R


A while back on StackOverflow, I answered a question on how: https://stackoverflow.com/questions/45682226/r-counting-how-many-polygons-between-two/47007246#47007246

Code

knitr::opts_chunk$set(eval=FALSE)
# Load packages
library(raster) # loads shapefile
library(igraph) # build network
library(spdep) # builds network
library(RColorBrewer)  # for plot colour palette
library(ggplot2) # plots results
    
# Load Data
powiaty <- shapefile("C:/Users/Mikey/Downloads/powiaty/powiaty")

Firstly the poly2nb function is used to calculate neighbouring regions:

# Find neighbouring areas
nb_q <- poly2nb(powiaty)

This creates our spatial mesh, which we can see here:

# Plot original results
coords <- coordinates(powiaty)
plot(powiaty)
plot(nb_q, coords, col="grey", add = TRUE)

This is the bit where I am not 100% sure what is happening. Basically, it is working out the shortest distance between all the shapefiles in the network, and returns a matrix of these pairs.

# Sparse matrix
nb_B <- nb2listw(nb_q, style="B", zero.policy=TRUE)
B <- as(nb_B, "symmetricMatrix")
    
# Calculate shortest distance
g1 <- graph.adjacency(B, mode="undirected")
dg1 <- diameter(g1)
sp_mat <- shortest.paths(g1)

Having made the calculations, the data can now be formatted to get into plotting format, so the shortest path matrix is merged with the spatial dataframe. I wasn’t sure what would be best to use as the ID for referring to datasets so I chose the jpt_kod_je variable.

# Name used to identify data
referenceCol <- powiaty$jpt_kod_je
  
# Rename spatial matrix
sp_mat2 <- as.data.frame(sp_mat)
sp_mat2$id <- rownames(powiaty@data)
names(sp_mat2) <- paste0("Ref", referenceCol)
          
# Add distance to shapefile data
powiaty@data <- cbind(powiaty@data, sp_mat2)
powiaty@data$id <- rownames(powiaty@data)

The data is now in a suitable format to display. Using the basic function spplot we can get a graph quite quickly:

displaylayer <- "Ref1261" # id for Krakow

# Plot the results as a basic spplot
spplot(powiaty, displaylayer)

I prefer ggplot for plotting more complex graphs as you can control the styling easier. However it is a bit more picky about how the data is fed into it, so we need to reformat the data for it before we build the graph:

# Or if you want to do it in ggplot
filtered <- data.frame(id = sp_mat2[,ncol(sp_mat2)], dist = sp_mat2[[displaylayer]]) 

ggplot_powiaty <- powiaty %>% fortify()
ggplot_powiaty <- merge(x = ggplot_powiaty, y = filtered, by = "id")
names(ggplot_powiaty)

And the plot. I have customised it a bit by removing elements which aren’t required and added a background. Also, to make the region at the centre of the search black, I subset the data using ggplot_powiaty[ggplot_powiaty$dist == 0, ], and then plot this as another polygon.

ggplot(ggplot_powiaty, aes(x = long, y = lat, group = group, fill = dist)) +
  geom_polygon(colour = "black") +
  geom_polygon(data =ggplot_powiaty[ggplot_powiaty$dist == 0, ],
               fill = "grey60") +
    labs(title = "Distance of Counties from Krakow", caption = "Mikey Harper") +
  scale_fill_gradient2(low = "#d73027", mid = "#fee08b", high = "#1a9850", midpoint = 10) +
  theme(
    axis.line = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank(),
    plot.background = element_rect(fill = "#f5f5f2", color = NA), 
    panel.background = element_rect(fill = "#f5f5f2", color = NA), 
    legend.background = element_rect(fill = "#f5f5f2", color = NA),
    panel.border = element_blank())