From 04597bf7bf0befcee663e36d9332023dc809674b Mon Sep 17 00:00:00 2001
From: Luuk van der Meer <luukvandermeer@live.nl>
Date: Sun, 17 Nov 2024 20:03:28 +0100
Subject: [PATCH] refactor: Make st_network_iso less dependent on GEOS version
 :construction:

---
 NAMESPACE                   |  1 +
 R/iso.R                     | 28 +++++++++++++++++++++-------
 man/st_network_iso.Rd       | 20 +++++++++++++++-----
 vignettes/sfn05_routing.qmd | 25 ++++---------------------
 4 files changed, 41 insertions(+), 33 deletions(-)

diff --git a/NAMESPACE b/NAMESPACE
index 6edab0d..2057b4b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -274,6 +274,7 @@ importFrom(sf,st_combine)
 importFrom(sf,st_concave_hull)
 importFrom(sf,st_contains)
 importFrom(sf,st_contains_properly)
+importFrom(sf,st_convex_hull)
 importFrom(sf,st_coordinates)
 importFrom(sf,st_covered_by)
 importFrom(sf,st_covers)
diff --git a/R/iso.R b/R/iso.R
index 97f8ff0..84e89a0 100644
--- a/R/iso.R
+++ b/R/iso.R
@@ -35,7 +35,8 @@
 #'
 #' @param ratio The ratio of the concave hull. Defaults to \code{1}, meaning
 #' that the convex hull is computed. See \code{\link[sf]{st_concave_hull}} for
-#' details. Ignored if \code{delineate = FALSE}.
+#' details. Ignored if \code{delineate = FALSE}. Setting this to a value
+#' smaller than 1 requires a GEOS version of at least 3.11.
 #'
 #' @param allow_holes May the concave hull have holes? Defaults to \code{FALSE}.
 #' Ignored if \code{delineate = FALSE}.
@@ -56,12 +57,21 @@
 #' oldpar = par(no.readonly = TRUE)
 #' par(mar = c(1,1,1,1))
 #'
-#' # Note that this function requires GEOS >= 3.11!
-#' if (compareVersion(sf_extSoftVersion()[["GEOS"]], "3.11.0") > -1) {
+#' center = st_centroid(st_combine(st_geometry(roxel)))
+#'
+#' net = as_sfnetwork(roxel, directed = FALSE)
+#'
+#' iso = net |>
+#'   st_network_iso(node_is_nearest(center), c(1000, 500, 250))
 #'
-#'   center = st_centroid(st_combine(st_geometry(roxel)))
+#' colors = c("#fee6ce90", "#fdae6b90", "#e6550d90")
 #'
-#'   net = as_sfnetwork(roxel, directed = FALSE)
+#' plot(net)
+#' plot(st_geometry(iso), col = colors, add = TRUE)
+#'
+#' # The level of detail can be increased with the ratio argument.
+#' # This requires GEOS >= 3.11.
+#' if (compareVersion(sf_extSoftVersion()[["GEOS"]], "3.11.0") > -1) {
 #'
 #'   iso = net |>
 #'     st_network_iso(node_is_nearest(center), c(1000, 500, 250), ratio = 0.3)
@@ -82,7 +92,7 @@ st_network_iso = function(x, node, cost, weights = edge_length(), ...,
 
 #' @importFrom methods hasArg
 #' @importFrom rlang enquo
-#' @importFrom sf st_combine st_concave_hull st_sf
+#' @importFrom sf st_combine st_concave_hull st_convex_hull st_sf
 #' @importFrom units as_units deparse_unit
 #' @export
 st_network_iso.sfnetwork = function(x, node, cost, weights = edge_length(),
@@ -115,7 +125,11 @@ st_network_iso.sfnetwork = function(x, node, cost, weights = edge_length(),
     in_iso = matrix[1, ] <= k
     iso = st_combine(node_geom[in_iso])
     if (delineate) {
-      iso = st_concave_hull(iso, ratio = ratio, allow_holes = allow_holes)
+      if (ratio == 1) {
+        iso = st_convex_hull(iso)
+      } else {
+        iso = st_concave_hull(iso, ratio = ratio, allow_holes = allow_holes)
+      }
     }
     iso
   }
diff --git a/man/st_network_iso.Rd b/man/st_network_iso.Rd
index 3cb7dc0..aca80f3 100644
--- a/man/st_network_iso.Rd
+++ b/man/st_network_iso.Rd
@@ -44,7 +44,8 @@ geometry is returned instead. Defaults to \code{TRUE}.}
 
 \item{ratio}{The ratio of the concave hull. Defaults to \code{1}, meaning
 that the convex hull is computed. See \code{\link[sf]{st_concave_hull}} for
-details. Ignored if \code{delineate = FALSE}.}
+details. Ignored if \code{delineate = FALSE}. Setting this to a value
+smaller than 1 requires a GEOS version of at least 3.11.}
 
 \item{allow_holes}{May the concave hull have holes? Defaults to \code{FALSE}.
 Ignored if \code{delineate = FALSE}.}
@@ -74,12 +75,21 @@ library(sf, quietly = TRUE)
 oldpar = par(no.readonly = TRUE)
 par(mar = c(1,1,1,1))
 
-# Note that this function requires GEOS >= 3.11!
-if (compareVersion(sf_extSoftVersion()[["GEOS"]], "3.11.0") > -1) {
+center = st_centroid(st_combine(st_geometry(roxel)))
+
+net = as_sfnetwork(roxel, directed = FALSE)
+
+iso = net |>
+  st_network_iso(node_is_nearest(center), c(1000, 500, 250))
 
-  center = st_centroid(st_combine(st_geometry(roxel)))
+colors = c("#fee6ce90", "#fdae6b90", "#e6550d90")
 
-  net = as_sfnetwork(roxel, directed = FALSE)
+plot(net)
+plot(st_geometry(iso), col = colors, add = TRUE)
+
+# The level of detail can be increased with the ratio argument.
+# This requires GEOS >= 3.11.
+if (compareVersion(sf_extSoftVersion()[["GEOS"]], "3.11.0") > -1) {
 
   iso = net |>
     st_network_iso(node_is_nearest(center), c(1000, 500, 250), ratio = 0.3)
diff --git a/vignettes/sfn05_routing.qmd b/vignettes/sfn05_routing.qmd
index 0e3346b..7f25e80 100644
--- a/vignettes/sfn05_routing.qmd
+++ b/vignettes/sfn05_routing.qmd
@@ -402,12 +402,14 @@ ggraph(net, "sf") +
   theme_void()
 ```
 
-<!-- ### Isodistance polygons
+### Isodistance polygons
 
 With respect to a given point $p$ and a given distance $d$, an isodistance is the line for which it holds that the distance from any point on the line from $p$ is equal to $d$. In spatial network analysis, it is common to find all nodes that fall within an isodistance line computed for a given source node $p$. By enclosing this nodes with a concave hull, we obtain a *isodistance polygon*. 
 
 This workflow is implemented in the function `st_network_iso()`. It will first compute the distance matrix from a given source node $p$ (which can be specified as explained [here](#specifying-origins-and-destinations)) to all other nodes in the network using `st_network_cost()`, and then filter those nodes that are within a given distance $d$. Using `sf::st_concave_hull()`, it will compute the concave hull around those nodes to obtain an isodistance polygon. Multiple distance thresholds can be given, resulting in multiple isodistance polygons being drawn. The output of the function is an `sf` object with one row per requested polygon.
 
+By setting the `ratio` argument (a value between 0 and 1, with 1 being the convex hull), we can increase or decrease the level of detail of the polygons. By default, no holes are allowed in the polygons, but this can be changed setting `allow_holes = TRUE`. Do note that changing the ratio requires a GEOS version of at least 3.11. GEOS is a C/C++ library for computational geometry used by `{sf}`.
+
 ```{r}
 # Define the source.
 # This point will be snapped to its nearest node.
@@ -427,25 +429,6 @@ ggraph(net, "sf") +
   theme_void()
 ```
 
-By setting the `ratio` argument (a value between 0 and 1, with 1 being the convex hull), we can increase or decrease the level of detail of the polygons. By default, no holes are allowed in the polygons, but this can be changed setting `allow_holes = TRUE`.
-
-```{r}
-polys = st_network_iso(
-  net, pt, rev(c(100, 250, 500, 1000)),
-  ratio = 0.3, allow_holes = TRUE
-)
-polys
-```
-
-```{r}
-ggraph(net, "sf") +
-  geom_edge_sf(color = "darkgrey") +
-  geom_node_sf(color = "darkgrey") +
-  geom_sf(data = polys, aes(fill = drop_units(cost)), alpha = 0.5) +
-  scale_fill_continuous("distance [m]") +
-  theme_void()
-```
-
 A different approach can be taken with the morpher `to_spatial_neighborhood()`. Instead of drawing polygons, this function will subset the network to contain only those nodes inside the isodistance line, and the edges that connect them.
 
 ```{r}
@@ -457,7 +440,7 @@ plot(net, col = "darkgrey")
 plot(sub_net, col = "orange", lwd = 2, cex = 1, add = TRUE)
 ```
 
-Both `st_network_iso()` and `to_spatial_neighborhood()` allow to specify other edge weights than geographic distance. This can be done through the `weights` argument, as explained [here](#specifying-edge-weights). When using time as edge weight, we talk about *isochrones* instead of *isodistances*. -->
+Both `st_network_iso()` and `to_spatial_neighborhood()` allow to specify other edge weights than geographic distance. This can be done through the `weights` argument, as explained [here](#specifying-edge-weights). When using time as edge weight, we talk about *isochrones* instead of *isodistances*.
 
 ### Custom routing