Merge pull request #14 from FishEcologyScience/CST-564-client-support

Handling points falling outside the polygon

Author: pbzonek

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: SAVM
 Title: Submerged aquatic vegetation model
-Version: 0.0.1.9004
-Date: 2025-07-15
+Version: 0.0.1.9005
+Date: 2025-08-12
 Authors@R: c(
     person(given = "Kevin",
            family = "Cazelles",

NAMESPACE

@@ -2,6 +2,7 @@
 
 S3method(preview_grid,sav_data)
 export(compute_fetch)
+export(identify_outsiders)
 export(invert_polygon)
 export(plot_sav_density)
 export(plot_sav_distribution)

NEWS.md

@@ -1,5 +1,6 @@
 # SAVM (devel)
 
+* `compute_fetch()` has a new argument `remove_outsider` to remove points falling outside the polygon. Also, `identify_outsiders()` helps visualize outsiders (see #13).
 * Model and plot functions have been adjusted to handle `sf` objects (see #11).
 * The element `mean_fetch` returned by `compute_fetch()` is now a `sf` object (see #9).
 * `preview_grid()` allows to preview grid (see #8).

R/compute_fetch.R

@@ -16,6 +16,10 @@
 #' @param crs {`object`}\cr{}
 #' Coordinate reference system (CRS) passed to [sf::st_crs()], used to
 #' transform `points` and `polygon`.
+#' @param remove_outsiders {`logical`}\cr{}
+#' Remove points outside the polygion. An error will be thrown if all points
+#' are rejected. Default is `FALSE`.
+#'
 #'
 #' @details
 #' Wind fetch is the unobstructed distance over which wind travels
@@ -88,7 +92,7 @@
 #' plot(res$transect_lines |> sf::st_geometry(), add = TRUE, col = 2, lwd = 0.5)
 #' }
 compute_fetch <- function(
-    points, polygon, max_dist = 15, n_bearings = 16, wind_weights = NULL, crs = NULL) {
+    points, polygon, max_dist = 15, n_bearings = 16, wind_weights = NULL, crs = NULL, remove_outsiders = FALSE) {
     valid_points(points)
     points$id_point <- seq_len(nrow(points))
     valid_polygon(polygon)
@@ -130,7 +134,7 @@ compute_fetch <- function(
         }
     }
 
-    valid_polygon_contains_points(points, polygon)
+    valid_polygon_contains_points(points, polygon, remove_outsiders)
 
     if (is.null(wind_weights)) {
         d_direction <- data.frame(
@@ -187,14 +191,32 @@ compute_fetch <- function(
                         )
                     ),
                 by = "id_point"
-            )  |>
+            ) |>
             dplyr::select(
                 c("id_point", "fetch_km", "weighted_fetch_km")
             ),
         transect_lines = transect_lines
     )
 }
 
+
+#' @describeIn compute_fetch Identify outsiders using a plot where points
+#' located outside the polygon are highlighted in red.
+#' @export
+identify_outsiders <- function(points, polygon) {
+    plot(polygon |> sf::st_geometry(), border = 1)
+    plot(
+        points |> sf::st_geometry(),
+        # there godd be more than one polygon
+        col = suppressMessages(
+            2 - apply(sf::st_within(points, polygon, sparse = FALSE), 1, any)
+        ),
+        pch = 19,
+        cex = 2,
+        add = TRUE
+    )
+}
+
 # HELPERS
 
 is_proj_unit_meter <- function(x) {
@@ -227,19 +249,33 @@ valid_polygon <- function(x) {
     }
 }
 
-valid_polygon_contains_points <- function(points, polygon) {
+valid_polygon_contains_points <- function(points, polygon, remove_outsiders = FALSE) {
+    if (remove_outsiders) {
+        points <- remove_points_outside_polygon(points, polygon)
+        if (!(points |> nrow())) {
+            rlang::abort("All points were outside the polygon considered.")
+        }
+    }
     chk <- suppressMessages({
         sf::st_contains(polygon, points, sparse = FALSE) |>
-            apply(2, any) |>
-            all()
+            apply(2, any)
     })
-    if (chk) {
+    if (all(chk)) {
         TRUE
     } else {
-        rlang::abort("`polygon` must include `points`.")
+        rlang::abort("`polygon` must include all points in `points`.
+        Use `remove_outsiders = TRUE` to remove points outside the polygon.
+        Alternatively use `identify_outsiders()` to vizualize outsiders.
+        ")
     }
 }
 
+remove_points_outside_polygon  <- function(points, polygon) {
+    suppressMessages(
+        points[apply(sf::st_within(points, polygon, sparse = FALSE), 1, any), ]
+    )
+}
+
 valid_direction <- function(direction) {
     if (!all(direction >= 0 & direction <= 360)) {
         rlang::abort("All directions must be within the range [0, 360].")

inst/example/le_in_out.geojson

@@ -0,0 +1,49 @@
+{
+    "type": "FeatureCollection",
+    "features": [
+        {
+            "type": "Feature",
+            "properties": {},
+            "geometry": {
+                "coordinates": [
+                    -81.36588837103685,
+                    42.242309732588126
+                ],
+                "type": "Point"
+            }
+        },
+        {
+            "type": "Feature",
+            "properties": {},
+            "geometry": {
+                "coordinates": [
+                    -80.77248459680594,
+                    42.1838784204125
+                ],
+                "type": "Point"
+            }
+        },
+        {
+            "type": "Feature",
+            "properties": {},
+            "geometry": {
+                "coordinates": [
+                    -81.4927603895497,
+                    42.96924260540226
+                ],
+                "type": "Point"
+            }
+        },
+        {
+            "type": "Feature",
+            "properties": {},
+            "geometry": {
+                "coordinates": [
+                    -80.45704257888741,
+                    43.00937083157342
+                ],
+                "type": "Point"
+            }
+        }
+    ]
+}

man/compute_fetch.Rd

@@ -12,6 +12,9 @@ le_pt_mid <- system.file("example", "le_middle.geojson", package = "SAVM") |>
 le_pt_out <- system.file("example", "le_land.geojson", package = "SAVM") |>
     sf::st_read(quiet = TRUE)
 
+le_pt_in_out <- system.file("example", "le_in_out.geojson", package = "SAVM") |>
+    sf::st_read(quiet = TRUE)
+
 
 test_that("helpers work", {
     expect_false(is_proj_unit_meter(4326))
@@ -39,13 +42,22 @@ test_that("helpers work", {
         "All directions must be within the range [0, 360].",
         fixed = TRUE
     )
+    #
+    expect_identical(
+        remove_points_outside_polygon(le_pt_in_out, le_bound),
+        le_pt_in_out[1:2, ]
+    )
+    expect_identical(
+        remove_points_outside_polygon(le_pt_out, le_bound) |> nrow(),
+        0L
+    )
 })
 
 
 test_that("helpers work", {
     withr::with_options(
         list(savm.verbose = "quiet"),
-        {   
+        {
             expect_error(
                 compute_fetch(le_pt, le_bound, n_bearings = 3),
                 "`n_bearings` should be equal or greater than 4."
@@ -59,8 +71,12 @@ test_that("helpers work", {
                 "Projection units must be meters."
             )
             expect_error(
-                compute_fetch(le_pt_out, le_bound_merc),
-                "`polygon` must include `points`"
+                compute_fetch(le_pt_in_out, le_bound_merc),
+                "`polygon` must include all points in `points`"
+            )
+            expect_error(
+                 compute_fetch(le_pt_out, le_bound_merc, remove_outsiders = TRUE),
+                 "All points were outside the polygon considered"
             )
         }
     )
@@ -82,14 +98,16 @@ test_that("compute_fetch() work", {
             expect_s3_class(res$transect_lines, "sf")
             expect_true(
                 all(
-                    c("direction", "weight", "transect_length", "rank") %in% 
-                    names(res$transect_lines)
+                    c("direction", "weight", "transect_length", "rank") %in%
+                        names(res$transect_lines)
                 )
             )
             # test with a point in the middle of the lake
             res2 <- compute_fetch(le_pt_mid, le_bound_merc, max_dist = 15)
             expect_equal(res2$mean_fetch$fetch_km, 15)
             expect_equal(res2$mean_fetch$weighted_fetch_km, 15)
+            res3 <- compute_fetch(le_pt_mid, le_bound_merc, max_dist = 15, remove_outsiders = TRUE)
+            expect_identical(res2, res3)
         }
     )
 })

tests/testthat/test-compute-fetch.R

@@ -127,6 +127,7 @@ study_depth <- stars::read_stars(system.file("example", "le_bathy.tiff", package
 > ***Note on spatial projections***: *The `SAVM` package works with spatial coordinates in meters and the WGS 84 / UTM zone 17N ([epsg:32617](https://epsg.io/32617)) projection as a default option since it was originally developped for use in the Great Lakes (WGS 84 / UTM zone 18N ([epsg:32618](https://epsg.io/32618)) is also a an option). Other viable options in a North American context are NAD83 / UTM zone 17N ([epsg:26917](https://epsg.io/26917)) and NAD83(CSRS) / UTM zone 17N ([egsg::2958](https://epsg.io/2958)).*
 
 
+
 ## Compute fetch and extract depth
 
 ### Compute fetch 
@@ -148,6 +149,9 @@ fetch$mean_fetch  |> sf::st_write("fetch.shp")
 ```
 
 
+> ***Note for points outside the polygon***: If any points fall outside the polygon, an error is raised. To remove outsiders, set `remove_outsider = TRUE`. Alternatively, you can identify outsiders using `identify_outsider()`.
+
+
 ### Extract depth 
 
 To extract depth data from a raster file for our study points, we use `st_extract()` from the package `stars`. Note that the points and the raster must use the same projection.