visualize mesh of polygon neighbors with tmap

I'm trying to make a map like the one below, using the {{tmap}} package.

That image is from Michael Harper's excellent blog post. It was made using the plot() function, as shown below.

# Loading example data
library(raster) # loads shapefile

# Data Analysis
library(spdep) # builds network

# Load Data
boundaries <- raster::getData(name = "GADM", country = "HUN", level = 2)

# Find neighbouring areas
nb_q <- poly2nb(boundaries)

# Plot original results
coords <- coordinates(boundaries)

# Show the results
plot(boundaries)
plot(nb_q, coords, col="grey", add = TRUE)

I can't figure out how to replicate the second plot() command in tmap. Is there a way to make tmap draw lines between a set of points using an adjacency list as created by spdep::poly2nb() or tmaptools::get_neighbours()?

Solution

There's spdep::nb2lines() to construct lines from the neighbours list, which you can then use as an extra layer for tmap. Following example uses sf for all geometries as this is what tmap preferes, but it should work with sp as well.

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(tmap)
library(spdep)

# cache downloaded GADM data in ~/geodata_dl
boundaries <- 
  geodata::gadm(country = "HUN", level = 2, path = "~/geodata_dl") |> 
  st_as_sf()

centroids <- st_geometry(boundaries) |> st_centroid()
nb_q <- poly2nb(boundaries)
nb_lines <- nb2lines(nb = nb_q, coords = centroids)

tm_shape(boundaries) +
  tm_polygons(col = "grey60", border.col = "grey75") +
tm_shape(nb_lines) + 
  tm_lines(alpha = .5, col = "gold") +
tm_shape(centroids) +
  tm_dots(size = .1, alpha = .5)

# input boundaries:
print(boundaries[,1:3], n = 3)
#> Simple feature collection with 168 features and 3 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 16.11384 ymin: 45.74783 xmax: 22.90558 ymax: 48.58638
#> Geodetic CRS:  WGS 84
#> First 3 features:
#>       GID_2 GID_0 COUNTRY                       geometry
#> 1 HUN.1.1_1   HUN Hungary POLYGON ((19.11987 46.03626...
#> 2 HUN.1.2_1   HUN Hungary POLYGON ((18.81133 45.91128...
#> 3 HUN.1.3_1   HUN Hungary POLYGON ((19.46291 46.17018...

# polygon centroids:
print(centroids, n = 3)
#> Geometry set for 168 features 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 16.32867 ymin: 45.85162 xmax: 22.66388 ymax: 48.45529
#> Geodetic CRS:  WGS 84
#> First 3 geometries:
#> POINT (19.31522 46.10193)
#> POINT (18.99711 46.14293)
#> POINT (19.31279 46.27723)

# nb2lines result:
print(nb_lines, n = 3)
#> Simple feature collection with 906 features and 5 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 16.32867 ymin: 45.85162 xmax: 22.66388 ymax: 48.45529
#> Geodetic CRS:  WGS 84
#> First 3 features:
#>   i j i_ID j_ID wt                       geometry
#> 1 1 2    1    2  1 LINESTRING (19.31522 46.101...
#> 2 1 3    1    3  1 LINESTRING (19.31522 46.101...
#> 3 1 8    1    8  1 LINESTRING (19.31522 46.101...

^{Created on 2023-10-24 with reprex v2.0.2}