Merge pull request #59 from CityRiverSpaces/52-corridor-valley-fn

Implement valley as initial method for corridor delineation

DESCRIPTION

@@ -37,7 +37,8 @@ Suggests:
     knitr,
     purrr,
     rmarkdown,
-    testthat (>= 3.0.0)
+    testthat (>= 3.0.0),
+    withr
 VignetteBuilder:
     knitr
 Config/testthat/edition: 3

NAMESPACE

@@ -23,7 +23,6 @@ export(get_slope)
 export(get_stac_asset_urls)
 export(get_utm_zone)
 export(get_valley)
-export(get_valley_mask)
 export(get_valley_polygon)
 export(get_valley_polygon_no_hole)
 export(get_valley_polygon_raw)

R/corridor.R

@@ -7,22 +7,25 @@
 #' @param bbox Bounding box defining the extent of the area of interest
 #' @param initial_method The method employed to define the initial river
 #'   corridor geometry. See [initial_corridor()] for the available methods
+#' @param buffer Buffer region to add to the river geometry to setup the initial
+#'   corridor (only used if `initial_method` is `"buffer"`)
+#' @param dem Digital elevation model (DEM) of the region (only used if
+#'   `initial_method` is `"valley"`)
 #' @param capping_method The method employed to connect the corridor edge end
 #'   points (i.e. to "cap" the corridor). See [cap_corridor()] for
 #'   the available methods
-#' TODO: add input arguments for the initialization method
 #'
 #' @return A simple feature geometry representing the river corridor
 #' @export
 get_corridor <- function(
-  network, river_centerline, river_surface, bbox, initial_method = "buffer",
-  capping_method = "direct"
+  network, river_centerline, river_surface, bbox, initial_method = "valley",
+  buffer = NULL, dem = NULL, capping_method = "direct"
 ) {
 
   # Draw the initial corridor geometry within the area of interest
   river <- c(river_centerline, river_surface)
-  # TODO: remove hardcoded buffer value here
-  corridor_init <- initial_corridor(river, initial_method, bbox, buffer = 1000)
+  corridor_init <- initial_corridor(river, method = initial_method,
+                                    buffer = buffer, dem = dem, bbox = bbox)
 
   # Find the corridor end points
   end_points <- corridor_end_points(river_centerline, bbox)
@@ -52,20 +55,33 @@ get_corridor <- function(
 #'
 #' @param river A simple feature geometry representing the river
 #' @param method The method employed to draw the initial river corridor:
-#'   - "buffer" (default): add a fixed buffer region to the river geometry
-#'     (see [river_buffer()])
-#' @param ... Additional arguments required by the function that implements the
-#'   chosen method (see `method`)
+#'   - "buffer": add a fixed buffer region to the river geometry (see
+#'     [river_buffer()])
+#'   - "valley" (default): use the river valley boundary, as estimated from the
+#'     provided digital elevation model (DEM, see [river_valley()])
+#' @param buffer Buffer region to add to the river geometry (only used if
+#'   `initial_method` is `"buffer"`)
+#' @param dem Digital elevation model (DEM) of the region (only used if
+#'   `initial_method` is `"valley"`)
 #' @param bbox Bounding box defining the extent of the area of interest
 #'
 #' @return A simple feature geometry
-initial_corridor <- function(river, method = "buffer", ..., bbox = NULL) {
-  args <- list(...)
+initial_corridor <- function(
+  river, method = "valley", buffer = NULL, dem = NULL, bbox = NULL
+) {
   if (method == "buffer") {
-    return(river_buffer(river, buffer = args$buffer, bbox = bbox))
+    if (is.null(buffer)) {
+      stop("Buffer should be provided if `method` is `'buffer'`")
+    }
+    return(river_buffer(river, buffer, bbox = bbox))
+  } else if (method == "valley") {
+    if (is.null(dem)) {
+      stop("DEM should be provided if `method` is `'valley'`")
+    }
+    return(river_valley(river, dem, bbox = bbox))
   } else {
     stop(
-      sprintf("Unknown method to initialize river corridor: {method}", method)
+      sprintf("Unknown method to initialize river corridor: %s", method)
     )
   }
 }
@@ -87,6 +103,20 @@ river_buffer <- function(river, buffer, bbox = NULL) {
   }
 }
 
+#' Draw a corridor as an estimate of the river valley.
+#'
+#' The valley is drawn on a digital elevation model (DEM) of the area,
+#' see [get_valley()] for details on the implementation.
+#'
+#' @param river A simple feature geometry representing the river
+#' @param dem The DEM of the area as a terra SpatRaster object
+#' @param bbox Bounding box defining the extent of the area of interest
+#'
+#' @return A simple feature geometry
+river_valley <- function(river, dem, bbox = NULL) {
+  get_valley(dem, river, bbox = bbox)
+}
+
 #' Find the corridor end points.
 #'
 #' Using the river center line and the bounding box of the area of interest,

R/delineate.R

@@ -8,6 +8,10 @@
 #'   effects close to its limits
 #' @param initial_method The method employed to define the initial river
 #'   corridor geometry. See [initial_corridor()] for the available methods
+#' @param initial_buffer Buffer region to add to the river geometry to setup the
+#'   initial corridor (only used if `initial_method` is `"buffer"`)
+#' @param dem Digital elevation model (DEM) of the region (only used if
+#'   `initial_method` is `"valley"`)
 #' @param capping_method The method employed to connect the corridor edge end
 #'   points (i.e. to "cap" the corridor). See [cap_corridor()] for
 #'   the available methods
@@ -16,28 +20,41 @@
 #'  [get_segments()]  and [rcoins::stroke()]. Only used if `segments` is TRUE.
 #' @param segments Whether to carry out the corridor segmentation
 #' @param riverspace Whether to carry out the riverspace delineation
+#' @param ... Additional (optional) input arguments for retrieving the DEM
+#'   dataset (see [get_dem()]). Only relevant if `initial_method` is `"valley"`
+#'   and `dem` is NULL
 #'
 #' @return A simple feature geometry
 #' @export
 delineate_corridor <- function(
   city_name, river_name, crs = NULL, bbox_buffer = NULL,
-  initial_method = "buffer", capping_method = "direct", angle_threshold = 90,
-  segments = FALSE, riverspace = FALSE
+  initial_method = "valley", initial_buffer = NULL, dem = NULL,
+  capping_method = "direct", angle_threshold = 90, segments = FALSE,
+  riverspace = FALSE, ...
 ) {
-  # Retrieve all relevant OSM datasets using the extended bounding box
-  osm_data <- get_osmdata(city_name, river_name, crs, bbox_buffer)
+  # Define the area of interest and (if not provided) the CRS
+  bbox <- get_osm_bb(city_name)
+  if (!is.null(bbox_buffer)) bbox <- buffer_bbox(bbox, buffer = bbox_buffer)
+  if (is.null(crs)) crs <- get_utm_zone(bbox)
 
-  # Define the area of interest
-  bbox <- sf::st_transform(osm_data$bb, sf::st_crs(osm_data$boundary))
+  # Retrieve all relevant OSM datasets within the area of interest
+  osm_data <- get_osmdata(bbox, city_name, river_name, crs = crs)
+
+  # If using the valley method, and the DEM is not provided, retrieve dataset
+  if (initial_method == "valley" && is.null(dem)) {
+    dem <- get_dem(bbox, crs = crs, ...)
+  }
 
   # Set up the combined street and rail network for the delineation
   network_edges <- dplyr::bind_rows(osm_data$streets, osm_data$railways)
   network <- as_network(network_edges)
 
   # Run the corridor delineation on the spatial network
+  bbox_repr <- reproject(bbox, crs)
   corridor <- get_corridor(
-    network, osm_data$river_centerline, osm_data$river_surface, bbox,
-    initial_method, capping_method
+    network, osm_data$river_centerline, osm_data$river_surface, bbox_repr,
+    initial_method = initial_method, buffer = initial_buffer, dem = dem,
+    capping_method = capping_method
   )
 
   if (segments) {

R/network.R

@@ -298,14 +298,20 @@ nearest_node <- function(network, target) {
 
 #' Subset a network keeping the only nodes that intersect a target geometry.
 #'
+#' If subsetting results in multiple disconnected components, we keep the main
+#' one.
+#'
 #' @param network A network object
 #' @param target The target geometry
 #'
 #' @return A spatial network object
 filter_network <- function(network, target) {
   network |>
     tidygraph::activate("nodes") |>
-    tidygraph::filter(sfnetworks::node_intersects(target))
+    tidygraph::filter(sfnetworks::node_intersects(target)) |>
+    # keep only the main connected component of the network
+    tidygraph::activate("nodes") |>
+    dplyr::filter(tidygraph::group_components() == 1)
 }
 
 #' Identify network edges that are intersecting a geometry

R/osmdata.R

@@ -37,37 +37,25 @@ get_osm_bb <- function(city_name) {
 #' the city boundary, the river centreline and surface, the streets, and the
 #' railways.
 #'
+#' @param bb Bounding box defining the area of interest
 #' @param city_name A character string with the name of the city.
 #' @param river_name A character string with the name of the river.
 #' @param crs An integer with the EPSG code for the projection. If no CRS is
 #'            specified, the default is the UTM zone for the city.
-#' @param buffer A numeric with the buffer distance in meters. By default,
-#'               no buffer is applied
 #'
 #' @return An list with the retrieved OpenStreetMap data sets for the
 #'         given location
 #' @export
 #'
 #' @examples
-#' get_osmdata("Bucharest", "Dambovita", buffer = 2000)
-get_osmdata <- function(city_name, river_name, crs = NULL, buffer = NULL) {
-  bb <- get_osm_bb(city_name)
-
-  if (is.null(crs)) crs <- get_utm_zone(bb)
-
-  if (!is.null(buffer)) {
-    bb <- bb |>
-      sf::st_as_sfc() |>
-      sf::st_transform(crs = crs) |>
-      sf::st_buffer(buffer) |>
-      sf::st_transform(crs = 4326) |>
-      sf::st_bbox()
-  }
-
-  boundary <- get_osm_city_boundary(city_name, bb, crs)
-  river <- get_osm_river(river_name, bb, crs)
-  srteets <- get_osm_streets(bb, crs)
-  railways <- get_osm_railways(bb, crs)
+#' bb <- get_osm_bb("Bucharest")
+#' crs <- get_utm_zone(bb)
+#' get_osmdata(bb, "Bucharest", "Dambovita", crs)
+get_osmdata <- function(bb, city_name, river_name, crs = NULL) {
+  boundary <- get_osm_city_boundary(bb, city_name, crs = crs)
+  river <- get_osm_river(bb, river_name, crs = crs)
+  srteets <- get_osm_streets(bb, crs = crs)
+  railways <- get_osm_railways(bb, crs = crs)
 
   osm_data <- list(
     bb = bb,
@@ -87,8 +75,8 @@ get_osmdata <- function(city_name, river_name, crs = NULL, buffer = NULL) {
 #' bounding box with the OSM tags "place:city" and "boundary:administrative".
 #' The result is filtered by the city name.
 #'
-#' @param city_name A character string with the name of the city
 #' @param bb Bounding box of class `bbox`
+#' @param city_name A character string with the name of the city
 #' @param crs Coordinate reference system as EPSG code
 #' @param multiple A logical indicating if multiple city boundaries should be
 #'                 returned. By default, only the first one is returned.
@@ -100,16 +88,15 @@ get_osmdata <- function(city_name, river_name, crs = NULL, buffer = NULL) {
 #' @examples
 #' bb <- get_osm_bb("Bucharest")
 #' crs <- get_utm_zone(bb)
-#' get_osm_city_boundary("Bucharest", bb, crs)
-get_osm_city_boundary <- function(city_name, bb, crs, multiple = FALSE) {
+#' get_osm_city_boundary(bb, "Bucharest", crs)
+get_osm_city_boundary <- function(bb, city_name, crs = NULL, multiple = FALSE) {
   # Define a helper function to fetch the city boundary
   fetch_boundary <- function(key, value) {
     CRiSp::osmdata_as_sf(key, value, bb)$osm_multipolygons |>
       dplyr::filter(
         .data$`name:en` == stringr::str_extract(city_name, "^[^,]+") |
           .data$name == stringr::str_extract(city_name, "^[^,]+")
       ) |>
-      sf::st_transform(crs) |>
       sf::st_geometry()
   }
 
@@ -128,6 +115,8 @@ get_osm_city_boundary <- function(city_name, bb, crs, multiple = FALSE) {
     stop("No city boundary found. The city name may be incorrect.")
   }
 
+  if (!is.null(crs)) city_boundary <- sf::st_transform(city_boundary, crs)
+
   if (length(city_boundary) > 1) {
     if (!multiple) {
       message("Multiple boundaries were found. Using the first one.")
@@ -142,8 +131,8 @@ get_osm_city_boundary <- function(city_name, bb, crs, multiple = FALSE) {
 
 #' Get the river centreline and surface from OpenStreetMap
 #'
-#' @param river_name The name of the river
 #' @param bb Bounding box of class `bbox`
+#' @param river_name The name of the river
 #' @param crs Coordinate reference system as EPSG code
 #'
 #' @return A list with the river centreline and surface
@@ -152,13 +141,14 @@ get_osm_city_boundary <- function(city_name, bb, crs, multiple = FALSE) {
 #' @examples
 #' bb <- get_osm_bb("Bucharest")
 #' crs <- get_utm_zone(bb)
-#' get_osm_river("Dâmbovița", bb, crs)
-get_osm_river <- function(river_name, bb, crs) {
+#' get_osm_river(bb, "Dâmbovița", crs)
+get_osm_river <- function(bb, river_name, crs = NULL) {
   # Get the river centreline
   river_centerline <- CRiSp::osmdata_as_sf("waterway", "river", bb)
   river_centerline <- river_centerline$osm_multilines |>
     dplyr::filter(.data$name == river_name) |>
-    sf::st_transform(crs) |>
+    # the query can return more features than actually intersecting the bb
+    sf::st_filter(sf::st_as_sfc(bb), .predicate = sf::st_intersects) |>
     sf::st_geometry()
 
   # Get the river surface
@@ -168,10 +158,14 @@ get_osm_river <- function(river_name, bb, crs) {
     sf::st_geometry() |>
     sf::st_as_sf() |>
     sf::st_crop(bb) |>
-    sf::st_transform(crs) |>
     sf::st_filter(river_centerline, .predicate = sf::st_intersects) |>
     sf::st_union()
 
+  if (!is.null(crs)) {
+    river_centerline <- sf::st_transform(river_centerline, crs)
+    river_surface <- sf::st_transform(river_surface, crs)
+  }
+
   return(list(centerline = river_centerline,
               surface = river_surface))
 }
@@ -191,7 +185,7 @@ get_osm_river <- function(river_name, bb, crs) {
 #' bb <- get_osm_bb("Bucharest")
 #' crs <- get_utm_zone(bb)
 #' get_osm_streets(bb, crs)
-get_osm_streets <- function(bb, crs, highway_values = NULL) {
+get_osm_streets <- function(bb, crs = NULL, highway_values = NULL) {
   if (is.null(highway_values)) {
     highway_values <- c("motorway", "trunk", "primary", "secondary", "tertiary")
   }
@@ -211,8 +205,9 @@ get_osm_streets <- function(bb, crs, highway_values = NULL) {
   streets_lines <- streets$osm_lines |>
     dplyr::bind_rows(poly_to_lines) |>
     dplyr::select("highway") |>
-    dplyr::rename(!!sym("type") := !!sym("highway")) |>
-    sf::st_transform(crs)
+    dplyr::rename(!!sym("type") := !!sym("highway"))
+
+  if (!is.null(crs)) streets_lines <- sf::st_transform(streets_lines, crs)
 
   return(streets_lines)
 }
@@ -230,12 +225,13 @@ get_osm_streets <- function(bb, crs, highway_values = NULL) {
 #' bb <- get_osm_bb("Bucharest")
 #' crs <- get_utm_zone(bb)
 #' get_osm_railways(bb, crs)
-get_osm_railways <- function(bb, crs) {
+get_osm_railways <- function(bb, crs = NULL) {
   railways <- osmdata_as_sf("railway", "rail", bb)
   railways_lines <- railways$osm_lines |>
     dplyr::select("railway") |>
-    dplyr::rename(!!sym("type") := !!sym("railway")) |>
-    sf::st_transform(crs)
+    dplyr::rename(!!sym("type") := !!sym("railway"))
+
+  if (!is.null(crs)) railways_lines <- sf::st_transform(railways_lines, crs)
 
   return(railways_lines)
 }