Request for finding closest feature

[MHaringa]

Firstly, great work with the package! For example I have the following cities:

library(s2)
city_names <- c("Vatican City", "San Marino", "Luxembourg")
cities <- s2_data_cities(city_names)

And I would like to find the closest city to each city:

city_names[s2_closest_feature(cities, cities)]
#> [1] "Vatican City" "San Marino"   "Luxembourg"

This gives the city name itself. I am not sure how I can use s2_closest_feature() to find the closest city, and not the city itself.
Ideally I would like to use for example (with no specification for the vector y):

city_names[s2_closest_feature(cities)]

Can this be added?

^{Created on 2021-05-08 by the reprex package (v1.0.0)}

[edzer]

Maybe this?

> sapply(seq_along(cities), function(i) setdiff(seq_along(cities), i)[s2_closest_feature(cities[i], cities[-i])])
[1] 2 1 2

[paleolimbot]

You're right that this is difficult to reproduce in s2! Is your use-case dealing with so many points that an vapply() or sapply()-based solution is going to be unreasonably slow? The following might be slightly faster:

library(s2)

city_names <- c("Vatican City", "San Marino", "Luxembourg")
cities <- s2_data_cities(city_names)
city_names[s2_closest_feature(cities, cities)]
#> [1] "Vatican City" "San Marino"   "Luxembourg"

ind <- vapply(
  seq_along(cities),
  function(i) s2_closest_feature(cities[i], cities[-i]),
  integer(1)
)

ifelse(ind >= seq_along(cities), ind + 1, ind)
#> [1] 2 1 2

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

[MHaringa]

Thanks for your solutions! My use-case indeed deals with so many points that looping takes a little while. An option to solve this kind of problems using s2 easily would be very much appreciated.

[edzer]

We can do this now in sf, which calls s2 for unprojected data - seems a common enough use case:

> library(sf)
Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
> library(s2)
> city_names <- c("Vatican City", "San Marino", "Luxembourg")
> cities <- st_as_sfc(s2_data_cities(city_names))
> st_nearest_feature(cities)
[1] 2 1 2

[paleolimbot]

This is related to k-nearest neighbours, which is also a common request and wasn't possible quite yet. To implement this I added s2_closest_edges() in #112...I think the current function has some rough edges but it can do what you'd like it to do about 20x faster than the vapply() version:

library(s2)
n <- 1e4
geog <- as_s2_geography(s2_point(runif(n, -1, 1), runif(n, -1, 1), runif(n, -1, 1)))

next_closest_native <- function(geog) unlist(s2_closest_edges(geog, geog, k = 2, min_distance = 0))

bench::mark(next_closest_native(geog))
#> # A tibble: 1 x 6
#>   expression                     min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr>                <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 next_closest_native(geog)    349ms    357ms      2.80     147KB        0

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

[edzer]

I believe you hit linear scaling here, where my Q&D does quadratic or worse...

[paleolimbot]

I think there's an even faster way to do it too ( https://github.com/r-spatial/s2/blob/master/inst/include/s2/s2closest_point_query.h ), but is more conducive to a dedicated S2Point data structure which we don't really have a system for yet.