This question is, in part, a follow up to: https://stackoverflow.com/questions/57269254/how-to-impute-missing-neighbours-of-a-spatial-weight-matrix-queen-contiguity

I have a .shp similar in structure to: Plot 1


columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
sp.sample <- columbus[c(1:5, 20, 21, 24, 27, 33, 40, 41, 42, 47),]

plot(st_geometry(sp.sample), border = "grey");
text(st_coordinates(st_centroid(sp.sample)), as.character(sp.sample$POLYID))

I use poly2nb to determine neighbours:

coords = st_coordinates(st_centroid(st_geometry(sp.sample)))
queen_nb = poly2nb(sp.sample, row.names=sp.sample$ID, queen=TRUE)

plot(st_geometry(sp.sample), border = "grey");
plot(queen_nb, coords, add = T)
text(st_coordinates(st_centroid(sp.sample)), as.character(sp.sample$POLYID))

Some polygons [e.g, 42] contain no neighbours. Other polygons [e.g, 21, 24] do share at least one neighbour, but are unconnected to the overall 'lattice'.

In the end, I want all polygons in the lattice to be connected in some way: plot of anticipated outcome

I have a function adapted from: https://stackoverflow.com/questions/57269254/how-to-impute-missing-neighbours-of-a-spatial-weight-matrix-queen-contiguity:

add_nb <- function(x){
  queen_nb <- poly2nb(x, queen = TRUE)
  count = card(queen_nb)
  ## get nearest neighbour index, use centroids:
  nnbs = knearneigh(st_coordinates(st_centroid(x)))$nn
  no_edges_from = which(count==0)
  for(i in no_edges_from){
    queen_nb[[i]] = nnbs[i]

new_queen_nb <- add_nb(sp.sample)

This allows me to join unlinked polygons to the lattice (i.e., those with no neighbours), but when polygons do have a neighbour (but are still unlinked to the main lattice) it does not join them to a neighbour.

The solution I am currently working with is to manually change the value in the matrix:

W <- nb2mat(new_queen_nb, style='B', zero.policy = T) # binary (0/1) weights 
colnames(W) <- sp.sample$POLYID
rownames(W) <- sp.sample$POLYID

W["24", "5"] <- 1
W["5", "24"] <- 1

This works fine, but is more impractical on a larger dataset. I have tried using snap = within poly2nb() and also st_buffer(), but in both cases, these are applied to all polygons, which I would prefer to avoid.

How do I adapt the function to join the remaining polygons to the main lattice?

I am new to working with spatial data in R so not that familiar with the different packages/functions.

Plot of real data for visualisation:

Plot of real data

  • 1
    In the case that one of your 'island' polygons is equidistant from multiple disconnected 'neighbors', how do you want the behavior of the algorithm to work?
    – Kartograaf
    Commented Oct 4, 2021 at 22:19
  • In that case, preferably, I would like the 'island' polygon to be joined to all 'neighbours' if they are of equal distance. Although, in my current function, the polygons with no neighbours just join to a single neighbour. Commented Oct 4, 2021 at 22:26
  • @Kartograaf Just to add, also, when I run: W <- nb2mat(new_queen_nb, style='B', zero.policy = T), the resulting matrix is non-symmetrical. i.e., in the case, above, row 42, col 21 <- 1, but row 21, col 42 <- 0. I am not too sure how to rectify this, but the matrix needs to be symmetrical in the analyses. Commented Oct 4, 2021 at 22:53

1 Answer 1


You can use the igraph package to construct a graph that connects all polygons in some way, then add this to the adjacency graph.

For example, use igraph to create a minimum spanning tree (MST) based on centroid distance, then add that to the adjacency graph. Duplicate edges (from the MST and the adjacency) will become single edges. All objects will be connected because the MST is connected and graph addition can't make a graph disconnected.

This does require a bit of data fiddling to get the distance matrix structure into and out of igraph correctly, but in the end you get a connected network of regions.

Here's the code - it requires the igraph package as well as spdep:

mstconnect <- function(polys, nb, distance="centroid"){
    if(distance == "centroid"){
        coords = sf::st_coordinates(sf::st_centroid(sf::st_geometry(polys)))
        dmat = as.matrix(dist(coords))
    }else if(distance == "polygon"){
        dmat = sf::st_distance(polys) + 1 # offset for adjacencies
        diag(dmat) = 0 # no self-intersections
        stop("Unknown distance method")
    gfull = igraph::graph.adjacency(dmat, weighted=TRUE, mode="undirected")
    gmst = igraph::mst(gfull)
    edgemat = as.matrix(igraph::as_adj(gmst))
    edgelistw = spdep::mat2listw(edgemat)
    edgenb = edgelistw$neighbour
    attr(edgenb,"region.id") = attr(nb, "region.id")
    allnb = spdep::union.nb(nb, edgenb)

It takes your polygons and an original neighbour object as input, and returns a new neighbour object with the additional links to make it fully connected. Make sure the polygons are the same ones used to make the second argument neighbour object, otherwise things will not line up properly.

mcadd = mstconnect(sp.sample,queen_nb)


plot(st_geometry(sp.sample), border = "grey");
plot(mcadd, coords, add = T,lwd=6, col="red")
plot(queen_nb, coords, add = T, col="blue",lwd=2)

enter image description here

Note the red lines are the fully connected network, the blue lines are the source "queen" neighbourhood. Looks icky purple where there's both but I think you get the idea.

(Note: This edited version can also compute the distance matrix based on the st_distance function, which returns the closest distance between two polygons. Sometimes this might be a better thing that using the centroid distance. For your test sample, this joins the three main islands in a slightly different way because the ordering of centroid distance between the three polygons joining the islands is different to the ordering of their nearest distance.

This is one of those "I'm sure I've done this before" codes but I couldn't find it, so I wrote it again.

  • Hi @Spacedman, this sounds like a promising solution. I will be out of the office for a couple of days, so won't be able to attempt this right away. If you do get round to generating and posting a working code later, I would be grateful, and I can then test it when I am next at my computer. Cheers. Commented Oct 5, 2021 at 10:34
  • @Spacedman thanks for the push here! I was working on a solution using igraph last night but got stuck because I didn't know how to describe the MST part. I will try to finish this tonight if nobody else has.
    – Kartograaf
    Commented Oct 5, 2021 at 20:59
  • @Kartograaf done! let me know if it works on your actual data!
    – Spacedman
    Commented Oct 6, 2021 at 7:37
  • @Spacedman it works like a charm, yet another job well done! I copied the line defining the variable 'coords' outside the function to replicate your plots.
    – Kartograaf
    Commented Oct 6, 2021 at 16:02
  • Spacedman, works perfectly! Thank you so much to you and @Kartograaf for your time on answering this question. Commented Oct 6, 2021 at 21:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.