# DRAFT :: Calculating Neighbouring Polygons in R

· by Mikey Harper · Read in about 3 min · (521 words) ·

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

# Load packages
library(igraph) # build network
library(spdep) # builds network
library(RColorBrewer)  # for plot colour palette
library(ggplot2) # plots results

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
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() ## Regions defined for each Polygons ggplot_powiaty <- merge(x = ggplot_powiaty, y = filtered, by = "id") names(ggplot_powiaty) ## [1] "id" "long" "lat" "order" "hole" "piece" "group" "dist" 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())