从 R 中的 Delaunay 三角剖分生成邻接矩阵

Question

我有一个数据框，其中包含如下坐标（纬度、经度）列表：

point lat long
1  51 31
2  52 31
3  52 30
4  56 28
5  57 29
6  53 32
7  54 35
8  52 32
9  48 30
10 49 27

我已经使用下面的代码成功生成了 Delaunay 三角剖分：

library(deldir)
vtess <- deldir(df$lat, df$long)
plot(vtess, wlines="triang", wpoints="none", number=FALSE, add=TRUE, lty=1)

我现在想做的是生成具有以下单元格值的邻接矩阵（10 x 10 矩阵）：

如果两个节点在 Delaunay 三角剖分中没有通过边连接：单元格的值 = 0
如果两个节点由 Delaunay 三角剖分中的一条边链接：单元格的值 = 两个节点之间的地理距离（使用 'geosphere' 包中的 distm() 和 DistVincenty 选项）

Answer 1

致OP：请注意评论；对于未来的 post，重要的是使您的 post 和代码独立。如果您不共享该转换代码，那么根据样本数据的转换（Delaunay 三角剖分）提出问题毫无意义。

除此之外，这是一个如何根据您的规范构建邻接矩阵的分步示例。为简单起见，我在这里假设 "distance between the two nodes" 你的意思是 Euclidean 距离。

让我们加载示例数据

df <- read.table(text = 
    "point lat long
1  51 31
2  52 31
3  52 30
4  56 28
5  57 29
6  53 32
7  54 35
8  52 32
9  48 30
10 49 27", header = T)

我们首先使用 deldir 包中的 deldir 执行 Delaunay 三角剖分。

library(deldir)
dxy <- deldir(df$lat, df$long)

让我们绘制结果

plot(df$lat, df$long, col = "red")
text(df$lat, df$long, df$point, cex = 0.5, col = "red", pos = 2)
plot(dxy, wlines = "triang", wpoints = "none", add = T)

接下来，我们从 Delaunay 三角剖分中提取顶点

# Extract the Delaunay vertices
vert <- data.frame(
    id1 = dxy$delsgs$ind1,
    id2 = dxy$delsgs$ind2)

我们计算所有连接节点之间的欧几里德距离，并将结果重塑为 data.frame

# Construct adjacency matrix
library(tidyverse)
dist.eucl <- function(x, y) sqrt(sum((x - y)^2))
df.adj <- vert %>%
    mutate_all(funs(factor(., levels = df$point))) %>%
    rowwise() %>%
    mutate(val = dist.eucl(df[id1, 2:3], df[id2, 2:3])) %>%
    ungroup() %>%
    complete(id1, id2, fill = list(val = 0)) %>%
    spread(id1, val)
## A tibble: 10 x 11
#   id2     `1`   `2`   `3`   `4`   `5`   `6`   `7`   `8`   `9`  `10`
#   <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 1      0.    1.00    0.  0.    0.    0.    0.    1.41    0.  0.
# 2 2      0.    0.      0.  0.    0.    0.    0.    1.00    0.  0.
# 3 3      1.41  1.00    0.  4.47  0.    2.24  0.    0.      4.  4.24
# 4 4      0.    0.      0.  0.    1.41  5.00  0.    0.      0.  0.
# 5 5      0.    0.      0.  0.    0.    5.00  6.71  0.      0.  0.
# 6 6      0.    1.41    0.  0.    0.    0.    3.16  1.00    0.  0.
# 7 7      0.    0.      0.  0.    0.    0.    0.    3.61    0.  0.
# 8 8      0.    0.      0.  0.    0.    0.    0.    0.      0.  0.
# 9 9      3.16  0.      0.  0.    0.    0.    7.81  4.47    0.  3.16
#10 10     0.    0.      0.  7.07  0.    0.    0.    0.      0.  0.

请注意，您可以将 dist.eucl. 替换为任何其他距离度量，例如haversine, cosine等，我选择dist.eucl只是为了方便

邻接关系matrix就是

df.adj %>% select(-id2) %>% as.matrix()
#            1        2 3        4        5        6        7        8 9
#[1,] 0.000000 1.000000 0 0.000000 0.000000 0.000000 0.000000 1.414214 0
#[2,] 0.000000 0.000000 0 0.000000 0.000000 0.000000 0.000000 1.000000 0
#[3,] 1.414214 1.000000 0 4.472136 0.000000 2.236068 0.000000 0.000000 4
#[4,] 0.000000 0.000000 0 0.000000 1.414214 5.000000 0.000000 0.000000 0
#[5,] 0.000000 0.000000 0 0.000000 0.000000 5.000000 6.708204 0.000000 0
#[6,] 0.000000 1.414214 0 0.000000 0.000000 0.000000 3.162278 1.000000 0
#[7,] 0.000000 0.000000 0 0.000000 0.000000 0.000000 0.000000 3.605551 0
#[8,] 0.000000 0.000000 0 0.000000 0.000000 0.000000 0.000000 0.000000 0
#[9,] 3.162278 0.000000 0 0.000000 0.000000 0.000000 7.810250 4.472136 0
#[10,] 0.000000 0.000000 0 7.071068 0.000000 0.000000 0.000000 0.000000 0
#           10
#[1,] 0.000000
#[2,] 0.000000
#[3,] 4.242641
#[4,] 0.000000
#[5,] 0.000000
#[6,] 0.000000
#[7,] 0.000000
#[8,] 0.000000
#[9,] 3.162278
#[10,] 0.000000

Answer 2

邻接矩阵在 Delaunay 三角剖分的输出中基本可用，只是需要稍微重新格式化。我们避免使用 distm 函数，因为我们不想计算所有点对之间的距离，只计算相邻点对之间的距离。直接调用距离函数效率更高

library(deldir)
library(geosphere)

del = deldir(dd$lat, dd$long)
del$delsgs$dist = with(del$delsgs, 
    distVincentySphere(p1 = cbind(y1, x1), p2 = cbind(y2, x2))
)
# we use y,x because the triangulation was lat,long but 
# distVincentySphere expects long,lat

# create empty adjacency matrix, fill in distances
adj = matrix(0, nrow = nrow(dd), ncol = nrow(dd))
adj[as.matrix(del$delsgs[c("ind1", "ind2")])] = del$delsgs$dist
round(adj)
#        [,1]  [,2]   [,3]   [,4]   [,5]   [,6]   [,7] [,8]   [,9]  [,10]
#  [1,]      0     0 131124      0      0      0      0    0 341685      0
#  [2,] 111319     0  68535      0      0 130321      0    0      0      0
#  [3,]      0     0      0      0      0      0      0    0      0      0
#  [4,]      0     0 464058      0      0      0      0    0      0 782155
#  [5,]      0     0      0 127147      0      0      0    0      0      0
#  [6,]      0     0 175378 422215 484616      0      0    0      0      0
#  [7,]      0     0      0      0 504301 227684      0    0 753748      0
#  [8,] 131124 68535      0      0      0 111319 299883    0 467662      0
#  [9,]      0     0 445278      0      0      0      0    0      0      0
# [10,]      0     0 395715      0      0      0      0    0 247685      0

使用此数据：

dd = read.table(text = "point lat long
1  51 31
2  52 31
3  52 30
4  56 28
5  57 29
6  53 32
7  54 35
8  52 32
9  48 30
10 49 27", header = T)

从 R 中的 Delaunay 三角剖分生成邻接矩阵

Generate an adjacency matrix from Delaunay Triangulation in R

r

graph

geo