I am using nngeo to identify for each town the nearest city. Library nngeo is available here.

I want to also report the distance between each town and its nearest city.

I wanted to use gDistance but am not sure how to properly specify it. Does anyone have an idea?

towns$dist<-apply(t(gDistance(towns,cities,byid = TRUE)),MARGIN = 1, FUN = "min") 


matched = st_join(towns, cities, join = nngeo::st_nn, maxdist = 5000, k = 1, progress = TRUE)
  • (there's only 6 nearest-neighbours within the range of 5000, is that what you want?)
    – Spacedman
    Aug 8, 2022 at 15:53
  • 1
    Use st_nn to get a list of nearest neighbours indexes, then lapply and other loops to get the distance via st_distance to the cities by indexing using the match vector. But what do you want done when there's no nearest point within the threshold in st_nn?
    – Spacedman
    Aug 8, 2022 at 16:01

1 Answer 1


I suggest using the solution presented in this answer (using the sf package): Getting a new column with distance to the nearest feature in R

I complete as requested in the comment. It is also possible to use such a code. It is slower because it is executed row by row. I assume that source_geometry and target_geometry are data in the form of sf with geometry column named geometry.

#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

#creating data

target_geometry <- st_sfc(st_point(c(1,0)), st_point(c(2,0)), st_point(c(3,0)), st_point(c(4,0)),
                 st_point(c(5,0)), st_point(c(6,0))) %>% 
  st_as_sf(tibble(name = LETTERS[1:6]))
source_geometry <- st_sfc(st_point(c(2,1)), st_point(c(5,1))) %>% 
  st_as_sf(tibble(name = 1:2))

st_geometry(target_geometry) <- "geometry"
st_geometry(source_geometry) <- "geometry"

source_geometry %>% 
  rowwise() %>% 
  mutate(nearest_id = st_nearest_feature(geometry, target_geometry),
         nearest_name = target_geometry$name[nearest_id],
         nearest_geom = st_geometry(target_geometry$geometry[nearest_id]),
         distance = st_distance(geometry, nearest_geom)) %>% 
  ungroup() %>% 
#> Simple feature collection with 2 features and 4 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 2 ymin: 1 xmax: 5 ymax: 1
#> CRS:           NA
#> # A tibble: 2 × 5
#>    name geometry nearest_id nearest_name distance[,1]
#>   <int>  <POINT>      <int> <chr>               <dbl>
#> 1     1    (2 1)          2 B                       1
#> 2     2    (5 1)          5 E                       1

Created on 2022-11-18 with reprex v2.0.2

  • Posting a link to another answer with no code or explanation from you is not an appropriate way to answer as indicated at this link: gis.stackexchange.com/help/how-to-answer
    – John Polo
    Nov 15, 2022 at 22:48
  • I supplemented the answer with another code that I was using for this purpose. I am sorry that your previous answer did not comply with the rules.
    – wacekk
    Nov 18, 2022 at 15:42
  • Did you mean to type "... my previous answer"? I did not give an answer. I was commenting on yours. No apologies are necessary. My comment was intended to be constructive, not punitive.
    – John Polo
    Nov 18, 2022 at 15:49
  • Yes, of course it was about my earlier answer. Sorry. Unfortunately, I can't edit a comment anymore.
    – wacekk
    Nov 18, 2022 at 15:51
  • No problem at all. :)
    – John Polo
    Nov 18, 2022 at 15:54

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.