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?

3 Answers 3


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
  • Thank you, but I get a warning when running the cbind: number of rows of result is not a multiple of vector length (arg 2). And my pljoin is then a list. It works if I take away the st_drop_geometry though, but with an extra shape field. How to avoid this?
    – Erik
    Feb 7, 2020 at 19:00
  • nearest should be the same length as rows in p, so l[nearest,] should be the same number of rows as p so the cbind should work. Hard to figure out without your data or more details.
    – Spacedman
    Feb 7, 2020 at 21:26

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.

  • This package is really helpful. I was getting very cryptic error messages when I try to use a combination of sf::st_distance() and sf::st_nearest_feature() to achieve the same thing. Mar 27, 2023 at 18:35

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')

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.