R: Counting how many polygons between two

Question

I was trying to recreate a map showing how many municipals are you away from Cracow:

and to change the city from Cracow to Wrocław. The map was done in GIMP.

I got a shapefile (available here: http://www.gis-support.pl/downloads/powiaty.zip). I read the shapefile documentation packages like maptools, rgdal or sf, but I couldn't find an automatic function to count it, because I wouldn't like to do that manually.

Is there a function to do that?

Credits: The map was done by Hubert Szotek on https://www.facebook.com/groups/mapawka/permalink/1850973851886654/

Answer 1

I am not that experienced at network analysis, so I must confess not to understand every single line of code as follows. But it works! A lot of the material was adapted from here: https://cran.r-project.org/web/packages/spdep/vignettes/nb_igraph.html

This is the final results:

Code

# 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("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$dist == 0

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

To plot for Wrocław as shown at the top of the post, just change displaylayer <- "Ref0264" and update the title.