# DRAFT :: Calculating Neighbouring Polygons in R

· by Mikey Harper · Read in about 3 min · (500 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

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())

Data scientist with interests in data visualisation, machine learning, geospatial modelling, energy and international development.

Latest Posts

Tag
Topic