Merge pull request #112 from r-spatial/dewey-dev

provide a path for k-nearest neighbours for points

Author: Dewey Dunnington

NAMESPACE

@@ -58,6 +58,7 @@ export(s2_bounds_rect)
 export(s2_buffer_cells)
 export(s2_centroid)
 export(s2_centroid_agg)
+export(s2_closest_edges)
 export(s2_closest_feature)
 export(s2_closest_point)
 export(s2_contains)

NEWS.md

@@ -1,5 +1,7 @@
 # s2 (development version)
 
+* Added `s2_closest_edges()` to make k-nearest neighbours calculation
+  possible on the sphere (#111, #112).
 * Added `s2_interpolate()`, `s2_interpolate_normalized()`, 
   `s2_project()`, and `s2_project_normalized()` to provide linear
   referencing support on the sphere (#96, #110).

R/RcppExports.R

@@ -125,6 +125,10 @@ cpp_s2_farthest_feature <- function(geog1, geog2) {
     .Call(`_s2_cpp_s2_farthest_feature`, geog1, geog2)
 }
 
+cpp_s2_closest_edges <- function(geog1, geog2, n, min_distance) {
+    .Call(`_s2_cpp_s2_closest_edges`, geog1, geog2, n, min_distance)
+}
+
 cpp_s2_may_intersect_matrix <- function(geog1, geog2, maxEdgesPerCell, maxFeatureCells, s2options) {
     .Call(`_s2_cpp_s2_may_intersect_matrix`, geog1, geog2, maxEdgesPerCell, maxFeatureCells, s2options)
 }

R/s2-matrix.R

@@ -11,6 +11,11 @@
 #' @inheritParams s2_contains
 #' @param x,y Geography vectors, coerced using [as_s2_geography()].
 #'   `x` is considered the source, where as `y` is considered the target.
+#' @param k The number of closest edges to consider when searching. Note
+#'   that in S2 a point is also considered an edge.
+#' @param min_distance The minimum distance to consider when searching for
+#'   edges. This filter is applied after the search is complete (i.e.,
+#'   may cause fewer than `k` values to be returned).
 #' @param max_edges_per_cell For [s2_may_intersect_matrix()],
 #'   this values controls the nature of the index on `y`, with higher values
 #'   leading to coarser index. Values should be between 10 and 50; the default
@@ -42,6 +47,11 @@
 #' # for each feature in x
 #' country_names[s2_farthest_feature(cities, countries)]
 #'
+#' # use s2_closest_edges() to find the k-nearest neighbours
+#' nearest <- s2_closest_edges(cities, cities, k = 2, min_distance = 0)
+#' city_names
+#' city_names[unlist(nearest)]
+#'
 #' # predicate matrices
 #' country_names[s2_intersects_matrix(cities, countries)[[1]]]
 #'
@@ -53,6 +63,13 @@ s2_closest_feature <- function(x, y) {
   cpp_s2_closest_feature(as_s2_geography(x), as_s2_geography(y))
 }
 
+#' @rdname s2_closest_feature
+#' @export
+s2_closest_edges <- function(x, y, k, min_distance = -1, radius = s2_earth_radius_meters()) {
+  stopifnot(k >= 1)
+  cpp_s2_closest_edges(as_s2_geography(x), as_s2_geography(y), k, min_distance / radius)
+}
+
 #' @rdname s2_closest_feature
 #' @export
 s2_farthest_feature <- function(x, y) {

man/s2_closest_feature.Rd

@@ -368,6 +368,20 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// cpp_s2_closest_edges
+List cpp_s2_closest_edges(List geog1, List geog2, int n, double min_distance);
+RcppExport SEXP _s2_cpp_s2_closest_edges(SEXP geog1SEXP, SEXP geog2SEXP, SEXP nSEXP, SEXP min_distanceSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< List >::type geog1(geog1SEXP);
+    Rcpp::traits::input_parameter< List >::type geog2(geog2SEXP);
+    Rcpp::traits::input_parameter< int >::type n(nSEXP);
+    Rcpp::traits::input_parameter< double >::type min_distance(min_distanceSEXP);
+    rcpp_result_gen = Rcpp::wrap(cpp_s2_closest_edges(geog1, geog2, n, min_distance));
+    return rcpp_result_gen;
+END_RCPP
+}
 // cpp_s2_may_intersect_matrix
 List cpp_s2_may_intersect_matrix(List geog1, List geog2, int maxEdgesPerCell, int maxFeatureCells, List s2options);
 RcppExport SEXP _s2_cpp_s2_may_intersect_matrix(SEXP geog1SEXP, SEXP geog2SEXP, SEXP maxEdgesPerCellSEXP, SEXP maxFeatureCellsSEXP, SEXP s2optionsSEXP) {
@@ -920,6 +934,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_s2_data_frame_from_s2_lnglat", (DL_FUNC) &_s2_data_frame_from_s2_lnglat, 1},
     {"_s2_cpp_s2_closest_feature", (DL_FUNC) &_s2_cpp_s2_closest_feature, 2},
     {"_s2_cpp_s2_farthest_feature", (DL_FUNC) &_s2_cpp_s2_farthest_feature, 2},
+    {"_s2_cpp_s2_closest_edges", (DL_FUNC) &_s2_cpp_s2_closest_edges, 4},
     {"_s2_cpp_s2_may_intersect_matrix", (DL_FUNC) &_s2_cpp_s2_may_intersect_matrix, 5},
     {"_s2_cpp_s2_contains_matrix", (DL_FUNC) &_s2_cpp_s2_contains_matrix, 3},
     {"_s2_cpp_s2_within_matrix", (DL_FUNC) &_s2_cpp_s2_within_matrix, 3},

src/RcppExports.cpp

@@ -165,6 +165,39 @@ IntegerVector cpp_s2_farthest_feature(List geog1, List geog2) {
   return op.processVector(geog1);
 }
 
+// [[Rcpp::export]]
+List cpp_s2_closest_edges(List geog1, List geog2, int n, double min_distance) {
+
+  class Op: public IndexedBinaryGeographyOperator<List, IntegerVector> {
+  public:
+    IntegerVector processFeature(Rcpp::XPtr<Geography> feature, R_xlen_t i) {
+      S2ClosestEdgeQuery query(this->geog2Index.get());
+      query.mutable_options()->set_max_results(n);
+      S2ClosestEdgeQuery::ShapeIndexTarget target(feature->ShapeIndex());
+      const auto& result = query.FindClosestEdges(&target);
+
+      // this code searches edges, which may come from the same feature
+      std::unordered_set<int> features;
+      for (S2ClosestEdgeQuery::Result res : result) {
+        if (res.distance().radians() > this->min_distance) {
+          features.insert(this->geog2IndexSource[res.shape_id()] + 1);
+        }
+      }
+
+      return IntegerVector(features.begin(), features.end());
+    }
+
+    int n;
+    double min_distance;
+  };
+
+  Op op;
+  op.n = n;
+  op.min_distance = min_distance;
+  op.buildIndex(geog2);
+  return op.processVector(geog1);
+}
+
 // ----------- indexed binary predicate operators -----------
 
 class IndexedMatrixPredicateOperator: public IndexedBinaryGeographyOperator<List, IntegerVector> {

src/s2-matrix.cpp

@@ -14,6 +14,36 @@ test_that("s2_closest|farthest_feature() works", {
   expect_identical(s2_data_tbl_countries$name[country_match_farthest], "New Zealand")
 })
 
+test_that("s2_closest_edges() works", {
+  expect_identical(
+    s2_closest_edges(
+      "POINT (0 0)",
+      c("POINT (0 0)", "POINT (0 1)", "POINT (0 2)", "POINT (0 3)"),
+      k = 1
+    ),
+    list(1L)
+  )
+
+  expect_identical(
+    s2_closest_edges(
+      "POINT (0 0)",
+      c("POINT (0 0)", "POINT (0 1)", "POINT (0 2)", "POINT (0 3)"),
+      k = 2,
+      min_distance = 0
+    ),
+    list(2L)
+  )
+
+  expect_identical(
+    s2_closest_edges(
+      "POINT (0 0)",
+      c("POINT (0 0)", "POINT (0 1)", "POINT (0 2)", "POINT (0 3)"),
+      k = 5
+    ) %>% lapply(sort),
+    list(1:4)
+  )
+})
+
 test_that("matrix predicates work", {
   expect_identical(
     s2_contains_matrix(