使用 tmap 可视化多边形邻居的网格

Question

我正在尝试使用 {{tmap}} 包制作如下所示的地图。

该图片来自 Michael Harper 的优秀博客文章。它是使用

plot()

函数制作的，如下所示。

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

我不知道如何复制 tmap 中的第二个

plot()

命令。有没有办法让 tmap 使用由

spdep::poly2nb()

或

tmaptools::get_neighbours()

创建的邻接列表在一组点之间绘制线条？

Answer 1

可以使用

spdep::nb2lines()

从邻居列表中构建线条，然后可以将其用作

tmap

的额外层。以下示例对所有几何图形使用

sf

，因为这是

tmap

更喜欢的，但它也应该与

sp

一起使用。

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...

^{创建于 2023-10-24，使用 reprex v2.0.2}