Getting a new column with distance to the nearest feature in R

Question

Having two sf data frames (a is points and b is polygons), I can get columns from the nearest feature by for instance:

c <- st_join(a, b, st_nearest_feature), 

But how is the best way to get an additional column with the distance as well?

Answer 1

Don't do the join. st_nearest_feature(a,b) will get you the index (row number) of the nearest feature in b to each feature in a.

EG using data p and l from ?st_nearest_feature made into sf data frame:

> (nearest = st_nearest_feature(p,l))
[1] 1 2 2 3

Then use st_distance to get the element-wise distances between each element of p and the corresponding element of l:

> (dist = st_distance(p, l[nearest,], by_element=TRUE))
[1] 0.10 0.01 0.01 0.10

You could use the nearest index to do the join like this:

> (pljoin = cbind(p, st_drop_geometry(l)[nearest,]))
Simple feature collection with 4 features and 2 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 0.1 ymin: -0.1 xmax: 0.1 ymax: 0.9
epsg (SRID):    NA
proj4string:    NA
  id st_drop_geometry.l..nearest...         geometry
1  a                              A POINT (0.1 -0.1)
2  b                              B POINT (0.1 0.11)
3  c                              B POINT (0.1 0.09)
4  d                              C  POINT (0.1 0.9)

(Column names are a bit mashed up but maybe you can work with that)

Then add the distance:

> pljoin$dist = dist
> 

Answer 2

You can use nngeo package to do that.

In particular, st_nn(a, b, k = 1, returnDist = T) returns both the nearest neighbor and distance. More generally, you can find j nearest neighbors with the corresponding distance by setting k = j in the st_nn argument.

Answer 3

So I found the above solutions to be extremely helpful, but very slow on large datasets. Also I found that weird implementations like st_distance in sf put off newcomers that want to use R for spatial analysis. Anyway, here is my solution that has a lot more code but trust me it's a lot faster than the standard st_distance function.

#add the coordinates for the points in SF dataframe A and B --------------
a_coord <- st_coordinates(a)
a <- cbind(a, a_coord)

b_coord <- st_coordinates(b)
b <- cbind(b, b_coord)


#get closest feature in B to A -----------------------------------------
A_B <- a %>%
  st_join(b %>%
          select(B_ID, X, Y) %>%
            rename(B_X = X, B_Y = Y), join = st_nearest_feature)



#create a WKT from the coords of A and closest feature in B --------------
A_B$line_wkt <- paste('linestring(',A_B$X,A_B$Y,',',A_B$B_X,A_B$B_Y,')')


#Convert WKT into Geom--------------------------------------
A_B <- A_B %>% 
  st_drop_geometry() %>%
  st_as_sf( wkt = 'line_wkt ', crs= 4326)


#Get the length (distance) of each line ----------------------------------
A_B$length <- as.numeric(st_length(A_B) )


#Join results with original A --------------------------------------------
a <- a %>%
  left_join(A_B %>%
              st_drop_geometry() %>%
              select(A_ID, B_ID, length), by = 'A_ID')