Generate an adjacency matrix from Delaunay Triangulation in R

I have a dataframe with a list of coordinates (lat, long) as below:

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


I already managed to generate the Delaunay triangulation using the code below:

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


What I would like to do now is to generate an adjacency matrix (10 by 10 matrix) having the following cell values:

How did you generate the triangulation? It would be better if we could start from where you are rather than have to figure that out too.
– G5W
2 days ago

2 Answers
2

To OP: Please note the comments; for future posts it is important to make your post and code self-contained. There is little point in asking a question based on a transformation (Delaunay triangulation) of your sample data, if you don't share that transformation code.

That aside, here is a step-by-step example how to construct an adjacency matrix according to your specifications. For simplicity, I here assume that by "distance between the two nodes" you mean Euclidean distance.

Let's load the sample data

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)


We first perform Delaunay triangulation using deldir from the deldir package.

deldir

deldir

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


Let's plot the result

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)


Next, we extract the vertices from the Delauney triangulation

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


We calculate Euclidean distances between all connected nodes, and reshape results in a data.frame

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.


Note that you can replace dist.eucl. with any other distance metric, e.g. Haversine, cosine, etc. I chose dist.eucl only because of convenience.

dist.eucl.

dist.eucl

The adjacency matrix is then simply

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


The adjacency matrix is essentially available in the the output of the Delaunay triangulation, it just needs a little reformatting. We avoid the distm function because we don't want to calculate the distance between all pairs of points, just adjacent pairs. It is more efficient to just call the distance function directly.

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


Using this data:

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)


By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.