use st_intersects rather than st_intersection for subsetting in point…

… indexing. fixes #11

NAMESPACE

@@ -105,7 +105,7 @@ importFrom(sf,st_crs)
 importFrom(sf,st_drop_geometry)
 importFrom(sf,st_geometry)
 importFrom(sf,st_geometry_type)
-importFrom(sf,st_intersection)
+importFrom(sf,st_intersects)
 importFrom(sf,st_is_longlat)
 importFrom(sf,st_join)
 importFrom(sf,st_linestring)

R/00_hydroloom.R

@@ -196,7 +196,7 @@ get_outlet_value <- function(x) {
 #' @importFrom dplyr ungroup n rename row_number arrange desc distinct mutate summarise
 #' @importFrom dplyr everything as_tibble pull group_split tibble bind_cols lag case_when
 #' @importFrom rlang :=
-#' @importFrom sf "st_geometry<-" st_drop_geometry st_geometry st_geometry_type st_intersection
+#' @importFrom sf "st_geometry<-" st_drop_geometry st_geometry st_geometry_type st_intersects
 #' @importFrom sf st_cast st_linestring st_is_longlat st_transform st_segmentize st_buffer
 #' @importFrom sf st_as_sf st_sf st_zm st_coordinates st_crs st_join st_reverse
 #' @importFrom pbapply pblapply pbsapply pbapply pboptions

R/index_points_to_lines.R

@@ -212,8 +212,9 @@ index_points_to_lines.hy <- function(x, points,
   search_radius <- as.numeric(search_radius) # everything in same units now
 
   if(!is.na(precision)) {
-    suppressWarnings(x <- st_intersection(x, point_buffer))
-    suppressWarnings(x <- st_cast(x, "LINESTRING", warn = FALSE))
+
+    x <- x[lengths(st_intersects(x, point_buffer, sparse = TRUE)) > 0, ]
+
   }
 
   x <- select(x, any_of(c(id, aggregate_id,
@@ -306,7 +307,7 @@ index_points_to_lines.hy <- function(x, points,
   x <- x |>
     add_index() |>
     filter(.data$L1 %in% matched$L1) |>
-    left_join(select(matched, all_of(c("L1", id))), by = "L1") |>
+    left_join(select(matched, all_of(c("L1", id))), by = "L1", relationship = "many-to-many") |>
     left_join(select(fline_atts, -"index"), by = id, relationship = "many-to-many")
 
   matched <- select(matched, point_id, node = "nn.idx", offset = "nn.dists", id)

tests/testthat/test_index.R

@@ -182,6 +182,35 @@ test_that("multipart indexing", {
 
 })
 
+test_that("no duplicates when using precision", {
+  # https://github.com/DOI-USGS/hydroloom/issues/11
+
+  # Define points
+  locs <- dplyr::tibble(id = c(1, 2),
+                 lon = c(-83.87865, -83.87975),
+                 lat = c(35.60989, 35.60963))
+  points <- sf::st_as_sf(locs, coords = c("lon", "lat"), crs = 4326) |>
+    sf::st_transform(5070)
+  flines_aoi <- sf::st_buffer(points, dist = units::as_units(5000, "m")) |>
+    sf::st_bbox() |>
+    sf::st_as_sfc()
+
+  # Define objects that we'd pass to get_flowline_index
+  search_rad <- units::as_units(500, "m")
+  max_match <- 3
+  precision <- 10
+  flines <- sf::read_sf(list.files(pattern = "index_precision.gpkg",
+                                   full.names = TRUE, recursive = TRUE))
+
+  check <-index_points_to_lines(x = flines, points = points,
+                        search_radius = search_rad,
+                        max_matches = max_match,
+                        precision = precision) |>
+    dplyr::filter(point_id == 2)
+
+  expect_equal(length(unique(check$comid)), 3)
+})
+
 test_that("disambiguate", {
 
   source(system.file("extdata", "sample_flines.R", package = "nhdplusTools"))