Applying a function by looping over two tables

We can take advantage of vectorization in R.

dists <- outer(seq_len(nrow(d1)), seq_len(nrow(d2)),
               function(i1, i2) sqrt( (d1$x1[i1] - d2$x2[i2])^2 + (d1$y1[i1] - d2$y2[i2])^2 ))
dists
#           [,1]     [,2]      [,3]     [,4]     [,5]      [,6]     [,7]      [,8]      [,9]    [,10]
#  [1,] 17.00680 16.62061 15.460223 19.05221 16.62596 13.787420 21.89659 15.518574 14.009078 20.04872
#  [2,] 10.91539 11.16448  9.307341 12.94013 10.49420  7.907477 15.80520  9.769963  7.975472 13.88798
#  [3,] 12.69908 12.35378 10.997895 14.71615 12.26294  9.341896 17.41379 11.123316  9.529537 15.61447
#  [4,] 14.02981 12.82053 12.119169 16.01171 13.54346 10.136304 18.36197 11.792099 10.467381 16.77871
#  [5,] 12.06452 14.20809 11.261432 14.12690 11.86161 10.671116 17.89579 12.611623 10.471053 15.47810
#  [6,] 13.69528 15.16659 12.696312 15.76172 13.44900 11.784031 19.32244 13.693289 11.706553 17.03575
#  [7,] 12.39585 14.42022 11.559708 14.45912 12.18613 10.909225 18.19487 12.845719 10.732153 15.79858
#  [8,] 15.48903 17.25052 14.668507 17.54801 15.29772 13.858714 21.30320 15.772602 13.757822 18.91564
#  [9,] 14.24652 13.16141 12.382742 16.23886 13.77327 10.446511 18.65747 12.116814 10.759542 17.03651
# [10,] 13.57476 14.44892 12.344098 15.63823 13.26645 11.181340 18.93203 13.057110 11.191215 16.79776
d1$dist <- apply(dists, 1, min)
d1$point2 <- d2$point[ apply(dists, 1, which.min) ]
d1
#           x1        y1 point      dist point2
# 1  -36.70045 -35.57741     a 13.787420      p
# 2  -31.71920 -31.95160     b  7.907477      p
# 3  -33.42082 -32.52115     c  9.341896      p
# 4  -35.01071 -32.14848     d 10.136304      p
# 5  -30.82340 -35.35468     e 10.471053      s
# 6  -32.48995 -35.91305     f 11.706553      s
# 7  -31.13327 -35.50498     g 10.732153      s
# 8  -33.14196 -37.89258     h 13.757822      s
# 9  -35.11322 -32.48860     i 10.446511      p
# 10 -33.04845 -34.96242     j 11.181340      p

The outer function creates an cartesian product of the two arguments, which we'll use as row indices on the two frames. The anonymous function is called once with (in this case) 100 integers in i1 and 100 integers in i2; if we debug the function, we'd see:

outer(seq_len(nrow(d1)), seq_len(nrow(d2)), function(i1, i2) { browser(); sqrt( (d1$x1[i1] - d2$x2[i2])^2 + (d1$y1[i1] - d2$y2[i2])^2 ); })
# Called from: FUN(X, Y, ...)
# Browse[1]> 
debug at #1: sqrt((d1$x1[i1] - d2$x2[i2])^2 + (d1$y1[i1] - d2$y2[i2])^2)
# Browse[2]> 
i1
#   [1]  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3
#  [64]  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10  1  2  3  4  5  6  7  8  9 10
# Browse[2]> 
i2
#   [1]  1  1  1  1  1  1  1  1  1  1  2  2  2  2  2  2  2  2  2  2  3  3  3  3  3  3  3  3  3  3  4  4  4  4  4  4  4  4  4  4  5  5  5  5  5  5  5  5  5  5  6  6  6  6  6  6  6  6  6  6  7  7  7
#  [64]  7  7  7  7  7  7  7  8  8  8  8  8  8  8  8  8  8  9  9  9  9  9  9  9  9  9  9 10 10 10 10 10 10 10 10 10 10

(Note the full-combination of i1-v-i2.)

This will calculate the distance from all points in d1 to all points in d2. Since you're talking euclidean, I thought it more efficient to do it manually instead of using dist, but you can use whatever function you like, with one caveat: it must accept all values at once.

For instance, if the distance function you want to use requires that its four arguments are length-1, then you will need to loop over all of the values (i1[1]-v-i2[1], then i1[2]-v-i2[1], etc). I don't think this is necessary or even advised if you can stick with vectorized calcs like this.

The output from outer is a matrix, with the rows here coming from the first argument (row indices for d1). (Pro tip: knowing which is which is easier to know and confirm when the two arguments are different lengths.)

We can then apply over the rows to find the shortest distance and the d2$point associated with it.

Using tidyverse strategy, it may be easier to understand (perhaps)

• I used purrr::map2_* functions to iterate d1$x1 and d1$y1 over all values of d2
• syntax inside map function is also written in tidyverse style
  • First, a distance variable say dist is mutated inside d2 for the given two arguments
  • Second, d2 is then filtered for min of dist
  • lastly, point is pulled for point and dist is pulled for distance.

d1 <- read.table(text = '       x1        y1      point
1  -36.70045 -35.57741     a
2  -31.71920 -31.95160     b
3  -33.42082 -32.52115     c
4  -35.01071 -32.14848     d
5  -30.82340 -35.35468     e
6  -32.48995 -35.91305     f
7  -31.13327 -35.50498     g
8  -33.14196 -37.89258     h
9  -35.11322 -32.48860     i
10 -33.04845 -34.96242     j', header = T)

d2 <- read.table(text = '      x2        y2      point
1  -21.84797 -27.29278     k
2  -27.96603 -21.43688     l
3  -24.57574 -25.98523     m
4  -20.21674 -26.02366     n
5  -22.54094 -26.86359     o
6  -27.66791 -25.16078     p
7  -20.78323 -20.54069     q
8  -27.26773 -23.25467     r
9  -26.76333 -25.70280     s
10 -20.38453 -23.92657     t', header = T)

suppressMessages(library(tidyverse))
d1 %>% mutate(nearest_point = map2_chr(x1, y1, ~ d2 %>% mutate(dist = sqrt((x2 - .x)^2 + (y2 - .y)^2)
                                                               ) %>% 
                                         filter(dist == min(dist)) %>%
                                         pull(point)
                                       ),
              shortest_dist = map2_dbl(x1, y1, ~ d2 %>% mutate(dist = sqrt((x2 - .x)^2 + (y2 - .y)^2)
                                                               ) %>% 
                                         filter(dist == min(dist)) %>%
                                         pull(dist)
                                       )
              )
#>           x1        y1 point nearest_point shortest_dist
#> 1  -36.70045 -35.57741     a             p     13.787420
#> 2  -31.71920 -31.95160     b             p      7.907477
#> 3  -33.42082 -32.52115     c             p      9.341896
#> 4  -35.01071 -32.14848     d             p     10.136304
#> 5  -30.82340 -35.35468     e             s     10.471053
#> 6  -32.48995 -35.91305     f             s     11.706553
#> 7  -31.13327 -35.50498     g             s     10.732153
#> 8  -33.14196 -37.89258     h             s     13.757822
#> 9  -35.11322 -32.48860     i             p     10.446511
#> 10 -33.04845 -34.96242     j             p     11.181340

^{Created on 2021-04-29 by the reprex package (v2.0.0)}

If you are interested in using other packages: Using laGP package:

dist_matrix <- sqrt(laGP::distance(d1[,1:2], d2[,1:2]))
id <- max.col(-dist_matrix)
min_distance <- dist_matrix[cbind(seq(nrow(d1)), id)]
cbind(d1, min_distance, min_point = d2[id, 3])

          x1        y1 point min_distance min_point
1  -36.70045 -35.57741     a    13.787420         p
2  -31.71920 -31.95160     b     7.907477         p
3  -33.42082 -32.52115     c     9.341896         p
4  -35.01071 -32.14848     d    10.136304         p
5  -30.82340 -35.35468     e    10.471053         s
6  -32.48995 -35.91305     f    11.706553         s
7  -31.13327 -35.50498     g    10.732153         s
8  -33.14196 -37.89258     h    13.757822         s
9  -35.11322 -32.48860     i    10.446511         p
10 -33.04845 -34.96242     j    11.181340         p

Here is another approach using only dplyr. It is based on a cross join d1 %>% full_join(d2, by = character()) to get all combinations and then we can calculate the distance (vectorized) and finally we filter the row with the smallest distance per point.x.

library(dplyr)

d1 %>%
  full_join(d2, by = character()) %>% 
  mutate(dist = sqrt((x1 - x2)^2 + (y1 - y2)^2)) %>% 
  group_by(point.x) %>% 
  filter(dist == min(dist))

#> # A tibble: 10 x 7
#> # Groups:   point.x [10]
#>       x1    y1 point.x    x2    y2 point.y  dist
#>    <dbl> <dbl> <chr>   <dbl> <dbl> <chr>   <dbl>
#>  1 -36.7 -35.6 a       -27.7 -25.2 p       13.8 
#>  2 -31.7 -32.0 b       -27.7 -25.2 p        7.91
#>  3 -33.4 -32.5 c       -27.7 -25.2 p        9.34
#>  4 -35.0 -32.1 d       -27.7 -25.2 p       10.1 
#>  5 -30.8 -35.4 e       -26.8 -25.7 s       10.5 
#>  6 -32.5 -35.9 f       -26.8 -25.7 s       11.7 
#>  7 -31.1 -35.5 g       -26.8 -25.7 s       10.7 
#>  8 -33.1 -37.9 h       -26.8 -25.7 s       13.8 
#>  9 -35.1 -32.5 i       -27.7 -25.2 p       10.4 
#> 10 -33.0 -35.0 j       -27.7 -25.2 p       11.2


# the data
d1 <- read.table(text = '       x1        y1      point
1  -36.70045 -35.57741     a
2  -31.71920 -31.95160     b
3  -33.42082 -32.52115     c
4  -35.01071 -32.14848     d
5  -30.82340 -35.35468     e
6  -32.48995 -35.91305     f
7  -31.13327 -35.50498     g
8  -33.14196 -37.89258     h
9  -35.11322 -32.48860     i
10 -33.04845 -34.96242     j', header = T)

d2 <- read.table(text = '      x2        y2      point
1  -21.84797 -27.29278     k
2  -27.96603 -21.43688     l
3  -24.57574 -25.98523     m
4  -20.21674 -26.02366     n
5  -22.54094 -26.86359     o
6  -27.66791 -25.16078     p
7  -20.78323 -20.54069     q
8  -27.26773 -23.25467     r
9  -26.76333 -25.70280     s
10 -20.38453 -23.92657     t', header = T)

^{Created on 2021-04-29 by the reprex package (v0.3.0)}