expose ways of iterating over vectors of s2 geometries?

[rsbivand]

I've started looking at supplementing code for finding point neighbours within distance thresholds on the sphere in spdep through sf::st_is_within_distance(), which is OK. However, to return distances for the chosen neighbours on the sphere, in line 64 in https://github.com/r-spatial/spdep/blob/master/R/nbdists.R I iterate over subsets of points through lapply() - see also the examples at the foot of https://r-spatial.github.io/spdep/reference/dnearneigh.html - it would be better to iterate in compiled code rather than use lapply(). Is LinkingTo feasible?

Finding k-nearest neighbours on the sphere looks even more complex for k > 1. Is there any such query in the s2 code itself? I looked at the documentation but couldn't see obvious opportunities, though I think they may be there somewhere, given that indexing is being used.

Current implementations use legacy C code for spdep::dnearneigh(), spdep::nbdists() and spdep::knearneigh() with Great Circle WGS84 ellipsoid distances and no indexing.

Any suggestions would be very welcome!

[edzer]

http://s2geometry.io/devguide/cpp/quickstart suggests that S2ClosestPointQuery can do both knn and distance-based queries, indexed.

[paleolimbot]

First, I think might find s2::s2_closest_edges() useful, which, despite the name, can do k-nearest neighbours (see #111 (comment) ). That could be modified to be more useful if it doesn't do exactly what you're trying to do.

Second, while it will probably be a long time before it's stable enough to be public, I've started a C API for internal use ( https://github.com/r-spatial/s2/blob/master/src/s2-c-api.cpp ). If S2 is going to gain traction it needs drop-in GEOS functions. Right now I'm just poking away at that as I run into things I need to implement in C...if you have ideas for things you'd like to see there I'd love ideas.

[edzer]

Looks like it! Here is an example use:

library(sf)
# Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(s2)
up = function(n) {
  x = runif(n)
  y = runif(n)
  st_as_sf(data.frame(x = x, y = y), coords = c("x", "y")) 
}
u1 = up(1000)
s2 = st_as_s2(u1)
s2_closest_edges(s2[1:10], s2, k = 3)
# [[1]]
# [1] 28 85  1

# [[2]]
# [1] 688 406   2

# [[3]]
# [1]  57 958   3

# [[4]]
# [1] 760 975   4

# [[5]]
# [1] 586 224   5

# [[6]]
# [1] 429 661   6

# [[7]]
# [1] 795 664   7

# [[8]]
# [1] 205  75   8

# [[9]]
# [1] 274 502   9

# [[10]]
# [1] 370 487  10

[rsbivand]

Thanks, implemented and pushed to github. With prototypes in place, I'll try to check timings (which look OK for distance and k-nearest neighbours, not OK for finding the distances for known neighbours - will look at multicore as a way round) and differences from the legacy code.

[paleolimbot]

If you give me a reprex of what is slow I can look into how to make it faster. Seems like a common/reasonable use case to make sure s2 can handle!

[rsbivand]

Thanks for being willing to take a look. This is an implicit reprex via remotes::install_github("r-spatial/spdep"):

library(spdep)
data(state)
us48.fipsno <- read.geoda(system.file("etc/weights/us48.txt",
 package="spdep")[1])
if (as.numeric(paste(version$major, version$minor, sep="")) < 19) {
 m50.48 <- match(us48.fipsno$"State.name", state.name)
} else {
 m50.48 <- match(us48.fipsno$"State_name", state.name)
}
xy <- as.matrix(as.data.frame(state.center))[m50.48,]
xy1 <- st_as_sf((as.data.frame(state.center))[m50.48,], coords=1:2,
  crs=st_crs("+proj=longlat +ellps=GRS80"))
old_use_s2 <- sf_use_s2()
sf_use_s2(TRUE)
gck1b <- knn2nb(knearneigh(xy1, k=1))
system.time(o <- nbdists(gck1b, xy1))
sf_use_s2(FALSE)
gck1c <- knn2nb(knearneigh(xy1, k=1))
system.time(o <- nbdists(gck1c, xy1))
sf_use_s2(old_use_s2)

and maybe debug(nbdists) to see:

https://github.com/r-spatial/spdep/blob/531b0d5e830653bb04197780f1cf4ffd1937aec7/R/nbdists.R#L61-L66

with a list of indices passed to nbdists() as "nb" (gck1* above) used to iterate over the distances to report (for point i, return distances (finally in "km") to neighbours id i in nb[[i]]. It did get a good deal faster after doing the coercion to "s2" from "sfc" once, not n times, but the n subsets and runs of s2_distance() in R seem to be costly. Please ask if this isn't clear. I think the same look-ups are what gstat would need too, for non-default nmax=, min= etc. for local kriging, after finding the k-nearest neighbours and where the metric distance is needed. In nbdists(), the "km" metric can be assigned in a separate lapply(), I think with little cost.

[paleolimbot]

I think you can get away with a single s2_distance() and a single units::set_units() call by flattening the neighbours object temporarily:

library(s2)
n <- 1e3
k <- 4
s2x <- as_s2_geography(s2_point(runif(n, -1, 1), runif(n, -1, 1), runif(n, -1, 1)))
nb <- s2_closest_edges(s2x, s2x, k + 1, min_distance = 0)

d1 <- function() {
    lapply(seq_along(nb), function(i) {
    s2::s2_distance(s2x[i], s2x[nb[[i]]])
  })
}

d2 <- function() {
  # you can abbreviate this if you know all the lengths are the same
  nb_flat_len <- vapply(nb, length, integer(1))
  dist_unnest <- s2::s2_distance(
    vctrs::vec_rep_each(s2x, nb_flat_len),
    s2x[unlist(nb)]
  )
  
  # maybe a faster way to do this nesting
  di_end <- cumsum(nb_flat_len)
  di_start <- c(1L, di_end[-length(di_end)] + 1L)
  d2 <- Map("[", list(dist_unnest), Map(":", di_start, di_end))
}

identical(d1(), d2())
#> [1] TRUE

bench::mark(d1(), d2())
#> Warning: Some expressions had a GC in every iteration; so filtering is disabled.
#> # A tibble: 2 x 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 d1()        93.24ms 101.96ms      9.94    6.38MB     17.9
#> 2 d2()         5.96ms   6.47ms    143.    529.28KB     11.9

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

[rsbivand]

Thanks! This is my 1970s Algol/Fortran playbook! A longer vector with start/end points - it should work fine. I'll need to work around length 1 neighbour vectors with 0L signalling no-neighbours on the "nb" object, but that is commonplace. Teaching today, will report back later on progress.

[rsbivand]

I have something of a puzzle:

> library(sf)
Linking to GEOS 3.10.0dev, GDAL 3.3.0, PROJ 8.0.1
> usc71k <- st_read("df_tracts.gpkg")
Reading layer `df_tracts' from data source 
  `/home/rsb/und/uam21/units_2/df_tracts.gpkg' using driver `GPKG'
Simple feature collection with 71353 features and 28 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -124.7332 ymin: 24.5447 xmax: -66.9499 ymax: 49.38436
Geodetic CRS:  NAD83
> usc71k_pts <- st_centroid(st_geometry(usc71k), of_largest_polygon=TRUE)
library(spdep)
system.time(k1 <- knn2nb(knearneigh(usc71k_pts, k=1)))
> library(spdep)
Loading required package: sp
Loading required package: spData
> system.time(k1 <- knn2nb(knearneigh(usc71k_pts, k=1)))
   user  system elapsed 
  6.006   0.026   6.405 
> system.time(k6 <- knn2nb(knearneigh(usc71k_pts, k=6)))
system.time(o <- nbdists(k1, usc71k_pts))
(all.linked <- max(unlist(o)))
   user  system elapsed 
  8.192   0.022   8.533 
> system.time(o <- nbdists(k1, usc71k_pts))
   user  system elapsed 
  0.583   0.001   0.594 
> (all.linked <- max(unlist(o)))
[1] 95.09906
> system.time(gc_nb <- dnearneigh(usc71k_pts, 0, all.linked+(.Machine$double.eps^(1/2))))
    user   system  elapsed 
3401.182    2.298 3431.844 
> st_write(usc71k_pts, "usc71k_pts.gpkg")
Writing layer `usc71k_pts' to data source `usc71k_pts.gpkg' using driver `GPKG'
Writing 71353 features with 0 fields and geometry type Point.

s2-based nbdists() and knearneigh() look pretty good (improved after dropping units, which turned out to be costly and not needed). But dnearneigh() is slower than the old GC by almost an order of magnitude, even after changing from st_is_within_distance() to s2_dwithin_matrix() (recent commit/push). Is s2_dwithin_matrix() using indexing? Points file (centroids of 71k US census tracts):
usc71k_pts.zip

[paleolimbot]

I've also found s2_dwithin() (regular or matrix) to be slow! It should be using indexing but might not be in the way I would expect.

As a workaround, I've used s2_buffer_cells() + s2_intersects() to preselect prior to running s2_dwithin(). (see https://argocanada.github.io/blog/posts/2021-06-04-finding-coastal-argo-trajectories/ , towards thee top).

[rsbivand]

Thanks. The case here is, however, finding neighbours within dist of each other. Maybe I'll rather try k-nearest neighbours then subset by distance, because s2::s2_closest_edges() obviously does use the indices. BruteForceMatrixPredicateOperator looks a bit brutal in cpp_s2_dwithin_matrix(), doesn't it?

[paleolimbot]

I think there's an easy option for that!

https://github.com/r-spatial/s2/blob/master/src/s2-matrix.cpp#L175

Can be updated to include

https://github.com/r-spatial/s2/blob/master/inst/include/s2/s2closest_edge_query.h#L135

[paleolimbot]

It might be too late to put this in version 1.0.6 (I'm running revdep checks now), but it will certainly make the next release.

[paleolimbot]

It's probably also worth checking if there's a faster way of doing the current dwithin, which should be much faster than it is.

s2/src/s2-predicates.cpp

Lines 134 to 136 in 1162a25

	S2ClosestEdgeQuery query(feature1->ShapeIndex());
	S2ClosestEdgeQuery::ShapeIndexTarget target(feature2->ShapeIndex());
	return query.IsDistanceLessOrEqual(&target, S1ChordAngle::Radians(this->distance[i]));

[rsbivand]

OK - I'll experiment with k-nearest neighbours, but that could be messy with this data set, there are very many distance neighbours in the urban areas compared to the longest distance between first nearest neighbours. I'll also look at buffering and intersecting, which seemed wasteful on a small data set but may use indexing.