Generate an adjacency matrix from Delaunay Triangulation in R

Question

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:

1. If the two nodes are NOT linked by an edge in the Delaunay triangulation: Value of the cell = 0
2. If the two nodes are linked by an edge in the Delaunay triangulation: Value of the cell = Geographic distance between the two nodes (using distm() from 'geosphere' package with DistVincenty option)

Answer 1

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.

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

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

   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)
   

   

3. Next, we extract the vertices from the Delaunay triangulation

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

4. We calculate Euclidean distances between all connected nodes, and reshape results in a 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.

5. The adjacency matrix is then simply

   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

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.

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)