Merge pull request #61 from r-spatial/dewey-attempt-matrix-index

Attempt indexed matrix

Author: paleolimbot

.gitignore

@@ -3,3 +3,4 @@
 .DS_Store
 .vscode
 docs
+configure.log

NAMESPACE

@@ -93,6 +93,7 @@ export(s2_make_line)
 export(s2_make_polygon)
 export(s2_max_distance)
 export(s2_max_distance_matrix)
+export(s2_may_intersect_matrix)
 export(s2_minimum_clearance_line_between)
 export(s2_num_points)
 export(s2_options)

R/RcppExports.R

@@ -17,6 +17,10 @@ cpp_s2_num_points <- function(geog) {
     .Call(`_s2_cpp_s2_num_points`, geog)
 }
 
+cpp_s2_is_empty <- function(geog) {
+    .Call(`_s2_cpp_s2_is_empty`, geog)
+}
+
 cpp_s2_area <- function(geog) {
     .Call(`_s2_cpp_s2_area`, geog)
 }
@@ -101,6 +105,10 @@ cpp_s2_farthest_feature <- function(geog1, geog2) {
     .Call(`_s2_cpp_s2_farthest_feature`, geog1, geog2)
 }
 
+cpp_s2_may_intersect_matrix <- function(geog1, geog2, maxEdgesPerCell, maxFeatureCells, s2options) {
+    .Call(`_s2_cpp_s2_may_intersect_matrix`, geog1, geog2, maxEdgesPerCell, maxFeatureCells, s2options)
+}
+
 cpp_s2_contains_matrix <- function(geog1, geog2, s2options) {
     .Call(`_s2_cpp_s2_contains_matrix`, geog1, geog2, s2options)
 }
@@ -113,14 +121,14 @@ cpp_s2_intersects_matrix <- function(geog1, geog2, s2options) {
     .Call(`_s2_cpp_s2_intersects_matrix`, geog1, geog2, s2options)
 }
 
-cpp_s2_disjoint_matrix <- function(geog1, geog2, s2options) {
-    .Call(`_s2_cpp_s2_disjoint_matrix`, geog1, geog2, s2options)
-}
-
 cpp_s2_equals_matrix <- function(geog1, geog2, s2options) {
     .Call(`_s2_cpp_s2_equals_matrix`, geog1, geog2, s2options)
 }
 
+cpp_s2_touches_matrix <- function(geog1, geog2, s2options) {
+    .Call(`_s2_cpp_s2_touches_matrix`, geog1, geog2, s2options)
+}
+
 cpp_s2_dwithin_matrix <- function(geog1, geog2, distance) {
     .Call(`_s2_cpp_s2_dwithin_matrix`, geog1, geog2, distance)
 }
@@ -133,6 +141,26 @@ cpp_s2_max_distance_matrix <- function(geog1, geog2) {
     .Call(`_s2_cpp_s2_max_distance_matrix`, geog1, geog2)
 }
 
+cpp_s2_contains_matrix_brute_force <- function(geog1, geog2, s2options) {
+    .Call(`_s2_cpp_s2_contains_matrix_brute_force`, geog1, geog2, s2options)
+}
+
+cpp_s2_within_matrix_brute_force <- function(geog1, geog2, s2options) {
+    .Call(`_s2_cpp_s2_within_matrix_brute_force`, geog1, geog2, s2options)
+}
+
+cpp_s2_intersects_matrix_brute_force <- function(geog1, geog2, s2options) {
+    .Call(`_s2_cpp_s2_intersects_matrix_brute_force`, geog1, geog2, s2options)
+}
+
+cpp_s2_disjoint_matrix_brute_force <- function(geog1, geog2, s2options) {
+    .Call(`_s2_cpp_s2_disjoint_matrix_brute_force`, geog1, geog2, s2options)
+}
+
+cpp_s2_equals_matrix_brute_force <- function(geog1, geog2, s2options) {
+    .Call(`_s2_cpp_s2_equals_matrix_brute_force`, geog1, geog2, s2options)
+}
+
 s2_point_from_numeric <- function(x, y, z) {
     .Call(`_s2_s2_point_from_numeric`, x, y, z)
 }
@@ -157,6 +185,10 @@ cpp_s2_contains <- function(geog1, geog2, s2options) {
     .Call(`_s2_cpp_s2_contains`, geog1, geog2, s2options)
 }
 
+cpp_s2_touches <- function(geog1, geog2, s2options) {
+    .Call(`_s2_cpp_s2_touches`, geog1, geog2, s2options)
+}
+
 cpp_s2_dwithin <- function(geog1, geog2, distance) {
     .Call(`_s2_cpp_s2_dwithin`, geog1, geog2, distance)
 }

R/s2-accessors.R

@@ -85,7 +85,7 @@ s2_num_points <- function(x) {
 #' @rdname s2_is_collection
 #' @export
 s2_is_empty <- function(x) {
-  s2_num_points(x) == 0
+  cpp_s2_is_empty(as_s2_geography(x))
 }
 
 #' @rdname s2_is_collection

R/s2-matrix.R

@@ -11,6 +11,16 @@
 #' @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 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
+#'   of 50 is adequate for most use cases, but for specialized operations users
+#'   may wish to use a lower value to increase performance.
+#' @param max_feature_cells For [s2_may_intersect_matrix()], this value
+#'   controls the approximation of `x` used to identify potential intersections
+#'   on `y`. The default value of 4 gives the best performance for most operations,
+#'   but for specialized operations users may wish to use a higher value to increase
+#'   performance.
 #'
 #' @return A vector of length `x`.
 #' @export
@@ -19,7 +29,7 @@
 #' See pairwise predicate functions (e.g., [s2_intersects()]).
 #'
 #' @examples
-#' city_names <- c("Vatican City", "San Marino", "Luxembourg", "Palikir")
+#' city_names <- c("Vatican City", "San Marino", "Luxembourg")
 #' cities <- s2_data_cities(city_names)
 #' country_names <- s2_data_tbl_countries$name
 #' countries <- s2_data_countries()
@@ -32,6 +42,9 @@
 #' # for each feature in x
 #' country_names[s2_farthest_feature(cities, countries)]
 #'
+#' # predicate matrices
+#' country_names[s2_intersects_matrix(cities, countries)[[1]]]
+#'
 #' # distance matrices
 #' s2_distance_matrix(cities, cities)
 #' s2_max_distance_matrix(cities, countries[1:4])
@@ -60,13 +73,13 @@ s2_max_distance_matrix <- function(x, y, radius = s2_earth_radius_meters()) {
 
 #' @rdname s2_closest_feature
 #' @export
-s2_contains_matrix <- function(x, y, options = s2_options()) {
+s2_contains_matrix <- function(x, y, options = s2_options(model = 0)) {
   cpp_s2_contains_matrix(as_s2_geography(x), as_s2_geography(y), options)
 }
 
 #' @rdname s2_closest_feature
 #' @export
-s2_within_matrix <- function(x, y, options = s2_options()) {
+s2_within_matrix <- function(x, y, options = s2_options(model = 0)) {
   cpp_s2_within_matrix(as_s2_geography(x), as_s2_geography(y), options)
 }
 
@@ -91,7 +104,11 @@ s2_intersects_matrix <- function(x, y, options = s2_options()) {
 #' @rdname s2_closest_feature
 #' @export
 s2_disjoint_matrix <- function(x, y, options = s2_options()) {
-  cpp_s2_disjoint_matrix(as_s2_geography(x), as_s2_geography(y), options)
+  # disjoint is the odd one out, in that it requires a negation of intersects
+  # this is inconvenient to do on the C++ level, and is easier to maintain
+  # with setdiff() here (unless somebody complains that this is slow)
+  intersection <- cpp_s2_intersects_matrix(as_s2_geography(x), as_s2_geography(y), options)
+  Map(setdiff, list(seq_along(y)), intersection)
 }
 
 #' @rdname s2_closest_feature
@@ -103,24 +120,51 @@ s2_equals_matrix <- function(x, y, options = s2_options()) {
 #' @rdname s2_closest_feature
 #' @export
 s2_touches_matrix <- function(x, y, options = s2_options()) {
-  x <- as_s2_geography(x)
-  y <- as_s2_geography(y)
-
-  options_closed <- options
-  options_closed$polygon_model <- 2
-  options_closed$polyline_model <- 2
-
-  options_open <- options
-  options_open$polygon_model <- 0
-  options_open$polyline_model <- 0
-
-  intersects_closed <- cpp_s2_intersects_matrix(x, y, options_closed)
-  intersects_open <- cpp_s2_intersects_matrix(x, y, options_open)
-  Map(setdiff, intersects_closed, intersects_open)
+  cpp_s2_touches_matrix(as_s2_geography(x), as_s2_geography(y), options)
 }
 
 #' @rdname s2_closest_feature
 #' @export
 s2_dwithin_matrix <- function(x, y, distance, radius = s2_earth_radius_meters()) {
   cpp_s2_dwithin_matrix(as_s2_geography(x), as_s2_geography(y), distance / radius)
 }
+
+#' @rdname s2_closest_feature
+#' @export
+s2_may_intersect_matrix <- function(x, y, max_edges_per_cell = 50, max_feature_cells = 4) {
+  cpp_s2_may_intersect_matrix(
+    as_s2_geography(x), as_s2_geography(y),
+    max_edges_per_cell, max_feature_cells,
+    s2_options()
+  )
+}
+
+# ------- for testing, non-indexed versions of matrix operators -------
+
+s2_contains_matrix_brute_force <- function(x, y, options = s2_options()) {
+  cpp_s2_contains_matrix_brute_force(as_s2_geography(x), as_s2_geography(y), options)
+}
+
+s2_within_matrix_brute_force <- function(x, y, options = s2_options()) {
+  cpp_s2_within_matrix_brute_force(as_s2_geography(x), as_s2_geography(y), options)
+}
+
+s2_covers_matrix_brute_force <- function(x, y, options = s2_options(model = 2)) {
+  cpp_s2_contains_matrix_brute_force(as_s2_geography(x), as_s2_geography(y), options)
+}
+
+s2_covered_by_matrix_brute_force <- function(x, y, options = s2_options(model = 2)) {
+  cpp_s2_within_matrix_brute_force(as_s2_geography(x), as_s2_geography(y), options)
+}
+
+s2_intersects_matrix_brute_force <- function(x, y, options = s2_options()) {
+  cpp_s2_intersects_matrix_brute_force(as_s2_geography(x), as_s2_geography(y), options)
+}
+
+s2_disjoint_matrix_brute_force <- function(x, y, options = s2_options()) {
+  cpp_s2_disjoint_matrix_brute_force(as_s2_geography(x), as_s2_geography(y), options)
+}
+
+s2_equals_matrix_brute_force <- function(x, y, options = s2_options()) {
+  cpp_s2_equals_matrix_brute_force(as_s2_geography(x), as_s2_geography(y), options)
+}

R/s2-predicates.R

@@ -159,18 +159,8 @@ s2_intersects_box <- function(x, lng1, lat1, lng2, lat2, detail = 1000, options
 #' @rdname s2_contains
 #' @export
 s2_touches <- function(x, y, options = s2_options()) {
-  x <- as_s2_geography(x)
-  y <- as_s2_geography(y)
-
-  options_closed <- options
-  options_closed$polygon_model <- 2
-  options_closed$polyline_model <- 2
-
-  options_open <- options
-  options_open$polygon_model <- 0
-  options_open$polyline_model <- 0
-
-  s2_intersects(x, y, options_closed) & !s2_intersects(x, y, options_open)
+  recycled <- recycle_common(as_s2_geography(x), as_s2_geography(y))
+  cpp_s2_touches(recycled[[1]], recycled[[2]], options)
 }
 
 #' @rdname s2_contains