diff --git a/DESCRIPTION b/DESCRIPTION
index c766f61..86dc3b4 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -37,7 +37,8 @@ Suggests:
     knitr,
     purrr,
     rmarkdown,
-    testthat (>= 3.0.0)
+    testthat (>= 3.0.0),
+    withr
 VignetteBuilder:
     knitr
 Config/testthat/edition: 3
diff --git a/NAMESPACE b/NAMESPACE
index c022bd7..7ee30b3 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/corridor.R b/R/corridor.R
index b61068d..f30a717 100644
--- a/R/corridor.R
+++ b/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,
diff --git a/R/delineate.R b/R/delineate.R
index f50bf3d..29f3a7b 100644
--- a/R/delineate.R
+++ b/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) {
diff --git a/R/network.R b/R/network.R
index 938b4f2..22e4274 100644
--- a/R/network.R
+++ b/R/network.R
@@ -298,6 +298,9 @@ 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
 #'
@@ -305,7 +308,10 @@ nearest_node <- function(network, target) {
 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
diff --git a/R/osmdata.R b/R/osmdata.R
index a62559d..e775335 100644
--- a/R/osmdata.R
+++ b/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,8 +88,8 @@ 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 |>
@@ -109,7 +97,6 @@ get_osm_city_boundary <- function(city_name, bb, crs, multiple = FALSE) {
         .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)
 }
diff --git a/R/utils.R b/R/utils.R
index 74be207..4566da1 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -30,12 +30,13 @@ get_utm_zone <- function(x) {
   return(epsg_code)
 }
 
-#' Get the bounding box from the x object. If the x does not have a CRS, WGS84
-#' is assumed.
+#' Get the bounding box from the x object
 #'
-#' @param x Simple feature object or a bounding box, provided either as a
-#'   matrix (with x, y as rows and min, max as columns) or as a vector (xmin,
-#'   ymin, xmax, ymax)
+#' If the x does not have a CRS, WGS84 is assumed.
+#'
+#' @param x Simple feature object (or compatible) or a bounding box, provided
+#'   either as a matrix (with x, y as rows and min, max as columns) or as a
+#'   vector (xmin, ymin, xmax, ymax)
 #' @return A bounding box as returned by [`sf::st_bbox()`]
 #' @export
 as_bbox <- function(x) {
@@ -45,6 +46,52 @@ as_bbox <- function(x) {
   }
   bbox <- sf::st_bbox(x)
   crs <- sf::st_crs(bbox)
-  if (is.na(crs)) sf::st_crs(bbox) <- sf::st_crs(4326)
+  if (is.na(crs)) sf::st_crs(bbox) <- sf::st_crs("EPSG:4326")
   return(bbox)
 }
+
+#' Apply a buffer region to a bounding box
+#'
+#' If the input bbox is in lat/lon coordinates, the buffer is approximately
+#' applied by first transforming the bbox in a suitable projected coordinate
+#' reference system, expanding it with the given buffer, transforming it back to
+#' the lat/lon system, and finally taking the bounding box of the obtained area.
+#'
+#' @param bbox Bonding box as a simple feature object
+#' @param buffer Buffer region in meters
+#' @return Expanded bounding box as a simple feature object
+buffer_bbox <- function(bbox, buffer) {
+  is_bbox_longlat <- sf::st_is_longlat(bbox)
+  bbox_sfc <- sf::st_as_sfc(bbox)
+  if (is_bbox_longlat) {
+    crs_meters <- get_utm_zone(bbox)
+    bbox_sfc <- sf::st_transform(bbox_sfc, crs_meters)
+  }
+  bbox_sfc_buffer <- sf::st_buffer(bbox_sfc, buffer)
+  if (is_bbox_longlat) {
+    bbox_sfc_buffer <- sf::st_transform(bbox_sfc_buffer, sf::st_crs(bbox))
+  }
+  bbox_buffer <- sf::st_bbox(bbox_sfc_buffer)
+  return(bbox_buffer)
+}
+
+#' Reproject a raster or vector dataset to the specified
+#' coordinate reference system (CRS)
+#'
+#' @param x Raster or vector object
+#' @param crs CRS to be projected to
+#' @param ... Optional arguments for raster or vector reproject functions
+#'
+#' @return Object reprojected to specified CRS
+#' @export
+reproject <- function(x, crs, ...) {
+  if (inherits(x, "SpatRaster")) {
+    # terra::crs does not support a numeric value as CRS, convert to character
+    if (inherits(crs, "numeric")) crs <- sprintf("EPSG:%s", crs)
+    return(terra::project(x, crs, ...))
+  } else if (inherits(x, c("bbox", "sfc", "sf"))) {
+    return(sf::st_transform(x, crs, ...))
+  } else {
+    stop(sprintf("Cannot reproject object type: %s", class(x)))
+  }
+}
diff --git a/R/valley.R b/R/valley.R
index bf0fa4a..8e810be 100644
--- a/R/valley.R
+++ b/R/valley.R
@@ -1,77 +1,106 @@
-#' Load dem from a STAC endpoint
+#' Endpoint and collection ID of the default STAC collection where to access
+#' digital elevation model (DEM) data. This is the global Copernicus DEM 30
+#' dataset hosted on AWS, as listed in the EarthSearch STAC API endpoint.
+#' Note that AWS credentials need to be set up in order to access the data (not
+#' the catalog).
+# nolint start
+#' References:
+#'  - [EarthSearch STAC API](https://element84.com/earth-search/)
+#'  - [Copernicus DEM](https://spacedata.copernicus.eu/collections/copernicus-digital-elevation-model)
+#'  - [AWS Copernicus DEM datasets](https://copernicus-dem-30m.s3.amazonaws.com/readme.html)
+#'  - [Data license](https://docs.sentinel-hub.com/api/latest/static/files/data/dem/resources/license/License-COPDEM-30.pdf)
+# nolint end
+default_stac_dem <- list(
+  endpoint = "https://earth-search.aws.element84.com/v1",
+  collection = "cop-dem-glo-30"
+)
+
+#' Access digital elevation model (DEM) for a given region
 #'
 #' @param bb A bounding box, provided either as a matrix (rows for "x", "y",
-#'   columns for "min", "max") or as a vector ("xmin", "ymin", "xmax", "ymax")
-#' @param resource from which to source dem. Can be "STAC".
-#' if "STAC" the parameters the following parameters
-#' must be supplied as named parameters. If omitted defaults
-#' will be used.
-#' @param ... Additional parameters to be passed depending on the resource.
-#'            In case `resource = "STAC"`, the arguments `endpoint` and
-#'            `collection` need to be passed to `get_stac_asset_urls()`.
-#'
-#' @return dem
+#'   columns for "min", "max") or as a vector ("xmin", "ymin", "xmax", "ymax"),
+#'   in lat/lon coordinates (WGS84 coordinate referece system)
+#' @param dem_source Source of the DEM:
+#'   - If "STAC" (default), DEM tiles are searched on a SpatioTemporal Asset
+#'     Catalog (STAC) end point, then accessed and mosaicked to the area of
+#'     interest
+#' @param stac_endpoint URL of the STAC API endpoint (only used if `dem_source`
+#'   is `"STAC"`). To be provided together with `stac_collection`, or leave
+#'   blank to use defaults (see [`default_stac_dem`])
+#' @param stac_collection Identifier of the STAC collection to be queried (only
+#'   used if `dem_source` is `"STAC"`). To be provided together with
+#'   `stac_endpoint`, or leave blank to use defaults (see [`default_stac_dem`])
+#' @param crs Coordinate reference system (CRS) which to transform the DEM to
+#'
+#' @return DEM as a terra SpatRaster object
 #' @export
-get_dem <- function(bb, resource = "STAC", ...) {
+get_dem <- function(bb, dem_source = "STAC", stac_endpoint = NULL,
+                    stac_collection = NULL, crs = NULL) {
   bbox <- as_bbox(bb)
-  args <- list(...)
-  if (resource == "STAC") {
-    if (length(args) && !is.null(...)) {
-      endpoint <- args$endpoint
-      collection <- args$collection
-      asset_urls <- get_stac_asset_urls(bbox,
-                                        endpoint = endpoint,
-                                        collection = collection)
-    } else {
-      asset_urls <- get_stac_asset_urls(bbox)
+  if (dem_source == "STAC") {
+    if (is.null(stac_endpoint) && is.null(stac_collection)) {
+      stac_endpoint <- default_stac_dem$endpoint
+      stac_collection <- default_stac_dem$collection
+    } else if (is.null(stac_endpoint) || is.null(stac_collection)) {
+      stop("Provide both or neither of `stac_endpoint` and `stac_collection`")
     }
+    asset_urls <- get_stac_asset_urls(bbox, endpoint = stac_endpoint,
+                                      collection = stac_collection)
     dem <- load_raster(bbox, asset_urls)
-    return(dem)
   } else {
-    stop(sprintf("Resource %s unknown", resource))
+    stop(sprintf("DEM source %s unknown", dem_source))
   }
+  if (!is.null(crs)) dem <- reproject(dem, crs)
+  return(dem)
 }
 
-#' Create vector/polygon representation of valley from dem and river polygon
-#' for a crs
+#' Extract the river valley from the DEM
 #'
-#' @param dem of region
-#' @param river vector/polygon representation of river area
-#' @param crs coordiante reference system to use
+#' The slope of the digital elevation model (DEM) is used as friction (cost)
+#' surface to compute the cost distance from any grid cell of the raster to
+#' the river. A characteristic value (default: the mean) of the cost distance
+#' distribution in a region surrounding the river (default: a buffer region of
+#' 2 km) is then calculated, and used to threshold the cost-distance surface.
+#' The resulting area is then "polygonized" to obtain the valley boundary as a
+#' simple feature geometry.
 #'
-#' @return (multi)polygon representation of valley
-#' area as st_geometry without holes
+#' @param dem Digital elevation model of the region
+#' @param river A simple feature geometry representing the river
+#' @param bbox Bounding box defining the extent of the area of interest
+#'
+#' @return River valley as a simple feature geometry
 #' @export
-get_valley <- function(dem, river, crs) {
-  dem_repr <- reproject(dem, crs)
-  river_repr <- reproject(river, crs)
-  cd_masked <- smooth_dem(dem_repr) |>
+get_valley <- function(dem, river, bbox = NULL) {
+  if (!terra::same.crs(dem, sf::st_crs(river)$wkt)) {
+    stop("DEM and river geometry should be in the same CRS")
+  }
+  if (!is.null(bbox)) {
+    bbox <- as_bbox(bbox)
+    dem <- terra::crop(dem, terra::ext(bbox))
+  }
+  cd_masked <- smooth_dem(dem) |>
     get_slope() |>
-    get_cost_distance(river_repr) |>
-    mask_cost_distance(river_repr)
+    get_cost_distance(river) |>
+    mask_cost_distance(river)
 
   cd_thresh <- get_cd_char(cd_masked)
 
-  valley <- get_valley_mask(cd_masked, cd_thresh) |>
-    get_valley_polygon()
+  valley <- get_valley_polygon(cd_masked < cd_thresh)
   return(valley)
 }
 
-#' Retrieve asset urls for the intersection of a bounding box with a
-#' remote STAC endpoint
+#' Retrieve the URLs of all the assets intersecting a bbox from a STAC API
 #'
 #' @param bb A bounding box, provided either as a matrix (rows for "x", "y",
-#'   columns for "min", "max") or as a vector ("xmin", "ymin", "xmax", "ymax")
-#' @param endpoint url of (remote) STAC endpoint
+#'   columns for "min", "max") or as a vector ("xmin", "ymin", "xmax", "ymax"),
+#'   in lat/lon coordinates (WGS84 coordinate referece system)
+#' @param endpoint URL of the STAC API endpoint
 #' @param collection STAC collection to be queried
 #'
-#' @return A list of urls for the assets in the collection
-#' overlapping with the specified bounding box
+#' @return A list of URLs for the assets in the collection overlapping with
+#'   the specified bounding box
 #' @export
-get_stac_asset_urls <- function(
-    bb,
-    endpoint = "https://earth-search.aws.element84.com/v1",
-    collection = "cop-dem-glo-30") {
+get_stac_asset_urls <- function(bb, endpoint, collection) {
   bbox <- as_bbox(bb)
   it_obj <- rstac::stac(endpoint) |>
     rstac::stac_search(collections = collection, bbox = bbox) |>
@@ -92,13 +121,18 @@ get_stac_asset_urls <- function(
 #' @export
 load_raster <- function(bb, raster_urlpaths) {
   bbox <- as_bbox(bb)
-  raster_urlpaths |>
+  cropped <- raster_urlpaths |>
     lapply(terra::rast) |>
-    lapply(terra::crop, terra::ext(bbox)) |>
-    do.call(terra::merge, args = _)
+    # snap spatial extent outward to include pixels crossed by the boundary
+    lapply(terra::crop, terra::ext(bbox), snap = "out")
+  if (length(cropped) > 1) {
+    return(do.call(terra::merge, args = cropped))
+  } else {
+    return(cropped[[1]])
+  }
 }
 
-#' Write dem to cloud optimized GeoTiff file as specified location
+#' Write DEM to cloud optimized GeoTiff file as specified location
 #'
 #' @param dem to write to file
 #' @param fpath filepath for output. If no output directory is specified
@@ -107,6 +141,8 @@ load_raster <- function(bb, raster_urlpaths) {
 #' @param output_directory where file should be written.
 #' If specified fpath is treated as filename only.
 #'
+#' @return The input DEM. This function is used for the side-effect of writing
+#'   values to a file.
 #' @export
 dem_to_cog <- function(dem, fpath, output_directory = NULL) {
   if (is.null(output_directory)) {
@@ -124,31 +160,11 @@ dem_to_cog <- function(dem, fpath, output_directory = NULL) {
                      overwrite = TRUE)
 }
 
-#' Reproject a raster or vector dataset to the specified
-#' coordinate reference system (CRS)
-#'
-#' @param x Raster or vector object
-#' @param crs CRS to be projected to
-#' @param ... Optional arguments for raster or vector reproject functions
-#'
-#' @return Object reprojected to specified CRS
-#' @export
-reproject <- function(x, crs, ...) {
-  if (inherits(x, "SpatRaster")) {
-    wkt <- sf::st_crs(crs)$wkt
-    return(terra::project(x, wkt, ...))
-  } else if (inherits(x, c("bbox", "sfc", "sf"))) {
-    return(sf::st_transform(x, crs, ...))
-  } else {
-    stop(sprintf("Cannot reproject object type: %s", class(x)))
-  }
-}
-
 #' Spatially smooth dem by (window) filtering
 #'
 #' @param dem raster data of dem
-#' @param method smoothing function to be used, e.g. "median".
-#' As accepted by focal
+#' @param method smoothing function to be used, e.g. "median",
+#'   as accepted by [terra::focal()]
 #' @param window size of smoothing kernel
 #'
 #' @return filtered dem
@@ -159,7 +175,8 @@ smooth_dem <- function(dem, method = "median", window = 5) {
   return(dem_smoothed)
 }
 
-#' Derive slope as percentage from dem
+#' Derive slope as percentage from DEM
+#'
 #' This makes use of the terrain function of the terra package
 #'
 #' @param dem raster data of dem
@@ -184,17 +201,17 @@ get_slope <- function(dem) {
 #'
 #' @export
 mask_slope <- function(slope, river, lthresh = 1.e-3, target = 0) {
-  slope_masked <- terra::mask(
-                              slope,
+  slope_masked <- terra::mask(slope,
                               terra::ifel(slope <= lthresh, NA, 1),
                               updatevalue = lthresh)
-
-  slope_masked <- terra::mask(
-                              slope_masked,
-                              terra::vect(river),
-                              inverse = TRUE,
-                              updatevalue = target,
-                              touches = TRUE)
+  for (ngeom in seq_len(length(sf::st_geometry(river)))) {
+    slope_masked <- terra::mask(slope_masked,
+                                terra::vect(river[ngeom]),
+                                inverse = TRUE,
+                                updatevalue = target,
+                                touches = TRUE)
+  }
+  return(slope_masked)
 }
 
 #' Derive cost distance function from masked slope
@@ -247,17 +264,6 @@ get_cd_char <- function(cd, method = "mean") {
   }
 }
 
-#' Select valley pixels from cost distance based on threshold
-#'
-#' @param cd cost distance raster
-#' @param thresh threshold cost distance value below which pixels are assuemd
-#' to belong to the valley
-#' @export
-get_valley_mask <- function(cd, thresh) {
-  valley_mask <- (cd < thresh)
-  return(valley_mask)
-}
-
 #' Create vector/polygon representation of valley raster mask
 #'
 #' @param valley_mask raster mask of valley pixels
@@ -292,8 +298,8 @@ get_valley_polygon_no_hole <- function(valley_polygon) {
 #'
 #' @param valley_mask raster mask of valley pixels
 #'
-#' @return (multi)polygon representation of valley area
-#' as st_geometry without holes
+#' @return (multi)polygon representation of valley area as a simple feature
+#'   geometry without holes
 #' @export
 get_valley_polygon <- function(valley_mask) {
   val_poly <- CRiSp::get_valley_polygon_raw(valley_mask) |>
diff --git a/README.Rmd b/README.Rmd
index 1690bce..8054410 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -43,11 +43,11 @@ river_name <- "Dâmbovița"
 epsg_code <- 32635
 
 # Get base layer for plotting
-bucharest_bb <- get_osm_bb(city_name) 
+bucharest_bb <- get_osm_bb(city_name)
 bucharest_streets <- get_osm_streets(bucharest_bb, epsg_code)[, "geometry"]
 
 # Delineate river corridor
-bucharest_river <- delineate_corridor("Bucharest", "Dâmbovița", epsg_code)
+bucharest_river <- delineate_corridor("Bucharest", "Dâmbovița", crs = epsg_code)
 
 # Plot results
 plot(bucharest_river, border = "orange", lwd = 3)
diff --git a/README.md b/README.md
index b2a1799..c619da4 100644
--- a/README.md
+++ b/README.md
@@ -34,11 +34,11 @@ river_name <- "Dâmbovița"
 epsg_code <- 32635
 
 # Get base layer for plotting
-bucharest_bb <- get_osm_bb(city_name) 
+bucharest_bb <- get_osm_bb(city_name)
 bucharest_streets <- get_osm_streets(bucharest_bb, epsg_code)[, "geometry"]
 
 # Delineate river corridor
-bucharest_river <- delineate_corridor("Bucharest", "Dâmbovița", epsg_code)
+bucharest_river <- delineate_corridor("Bucharest", "Dâmbovița", crs = epsg_code)
 #> Warning: to_spatial_subdivision assumes attributes are constant over geometries
 #> Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
 
diff --git a/data-raw/bucharest.R b/data-raw/bucharest.R
index d53d788..4f6b5b0 100644
--- a/data-raw/bucharest.R
+++ b/data-raw/bucharest.R
@@ -5,8 +5,10 @@ epsg_code <- 32635
 bbox_buffer <- 2000
 
 # Fetch the OSM data
-bucharest_osm <- CRiSp::get_osmdata(city_name, river_name,
-                                    crs = epsg_code, buffer = bbox_buffer)
+bbox <- get_osm_bb(city_name)
+bbox_expanded <- buffer_bbox(bbox, bbox_buffer)
+bucharest_osm <- CRiSp::get_osmdata(bbox_expanded, city_name, river_name,
+                                    crs = epsg_code)
 
 # Fix encoding issue in the WKT strings
 fix_wkt_encoding <- function(x) {
diff --git a/data/bucharest_dem.rda b/data/bucharest_dem.rda
index 1b7f2d8..75e70b9 100644
Binary files a/data/bucharest_dem.rda and b/data/bucharest_dem.rda differ
diff --git a/data/bucharest_osm.rda b/data/bucharest_osm.rda
index fd186d4..214d805 100644
Binary files a/data/bucharest_osm.rda and b/data/bucharest_osm.rda differ
diff --git a/man/as_bbox.Rd b/man/as_bbox.Rd
index f7f4a7f..fcff840 100644
--- a/man/as_bbox.Rd
+++ b/man/as_bbox.Rd
@@ -2,20 +2,18 @@
 % Please edit documentation in R/utils.R
 \name{as_bbox}
 \alias{as_bbox}
-\title{Get the bounding box from the x object. If the x does not have a CRS, WGS84
-is assumed.}
+\title{Get the bounding box from the x object}
 \usage{
 as_bbox(x)
 }
 \arguments{
-\item{x}{Simple feature object or a bounding box, provided either as a
-matrix (with x, y as rows and min, max as columns) or as a vector (xmin,
-ymin, xmax, ymax)}
+\item{x}{Simple feature object (or compatible) or a bounding box, provided
+either as a matrix (with x, y as rows and min, max as columns) or as a
+vector (xmin, ymin, xmax, ymax)}
 }
 \value{
 A bounding box as returned by \code{\link[sf:st_bbox]{sf::st_bbox()}}
 }
 \description{
-Get the bounding box from the x object. If the x does not have a CRS, WGS84
-is assumed.
+If the x does not have a CRS, WGS84 is assumed.
 }
diff --git a/man/buffer_bbox.Rd b/man/buffer_bbox.Rd
new file mode 100644
index 0000000..7c99730
--- /dev/null
+++ b/man/buffer_bbox.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{buffer_bbox}
+\alias{buffer_bbox}
+\title{Apply a buffer region to a bounding box}
+\usage{
+buffer_bbox(bbox, buffer)
+}
+\arguments{
+\item{bbox}{Bonding box as a simple feature object}
+
+\item{buffer}{Buffer region in meters}
+}
+\value{
+Expanded bounding box as a simple feature object
+}
+\description{
+If the input bbox is in lat/lon coordinates, the buffer is approximately
+applied by first transforming the bbox in a suitable projected coordinate
+reference system, expanding it with the given buffer, transforming it back to
+the lat/lon system, and finally taking the bounding box of the obtained area.
+}
diff --git a/man/default_stac_dem.Rd b/man/default_stac_dem.Rd
new file mode 100644
index 0000000..1dfcc71
--- /dev/null
+++ b/man/default_stac_dem.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/valley.R
+\docType{data}
+\name{default_stac_dem}
+\alias{default_stac_dem}
+\title{Endpoint and collection ID of the default STAC colleciton where to access
+digital elevation model (DEM) data. This is the global Copernicus DEM 30
+dataset hosted on AWS, as listed in the EarthSearch STAC API endpoint.
+Note that AWS credentials need to be setup in order to access the data (not
+the catalog).
+References:
+\itemize{
+\item \href{https://element84.com/earth-search/}{EarthSearch STAC API}
+\item \href{https://spacedata.copernicus.eu/collections/copernicus-digital-elevation-model}{Copernicus DEM}
+\item \href{https://copernicus-dem-30m.s3.amazonaws.com/readme.html}{AWS Copernicus DEM datasets}
+\item \href{https://docs.sentinel-hub.com/api/latest/static/files/data/dem/resources/license/License-COPDEM-30.pdf}{Data license}
+}}
+\format{
+An object of class \code{list} of length 2.
+}
+\usage{
+default_stac_dem
+}
+\description{
+Endpoint and collection ID of the default STAC colleciton where to access
+digital elevation model (DEM) data. This is the global Copernicus DEM 30
+dataset hosted on AWS, as listed in the EarthSearch STAC API endpoint.
+Note that AWS credentials need to be setup in order to access the data (not
+the catalog).
+References:
+\itemize{
+\item \href{https://element84.com/earth-search/}{EarthSearch STAC API}
+\item \href{https://spacedata.copernicus.eu/collections/copernicus-digital-elevation-model}{Copernicus DEM}
+\item \href{https://copernicus-dem-30m.s3.amazonaws.com/readme.html}{AWS Copernicus DEM datasets}
+\item \href{https://docs.sentinel-hub.com/api/latest/static/files/data/dem/resources/license/License-COPDEM-30.pdf}{Data license}
+}
+}
+\keyword{datasets}
diff --git a/man/delineate_corridor.Rd b/man/delineate_corridor.Rd
index 3d979af..d2a7c2d 100644
--- a/man/delineate_corridor.Rd
+++ b/man/delineate_corridor.Rd
@@ -9,11 +9,14 @@ delineate_corridor(
   river_name,
   crs = NULL,
   bbox_buffer = NULL,
-  initial_method = "buffer",
+  initial_method = "valley",
+  initial_buffer = NULL,
+  dem = NULL,
   capping_method = "direct",
   angle_threshold = 90,
   segments = FALSE,
-  riverspace = FALSE
+  riverspace = FALSE,
+  ...
 )
 }
 \arguments{
@@ -30,6 +33,12 @@ effects close to its limits}
 \item{initial_method}{The method employed to define the initial river
 corridor geometry. See \code{\link[=initial_corridor]{initial_corridor()}} for the available methods}
 
+\item{initial_buffer}{Buffer region to add to the river geometry to setup the
+initial corridor (only used if \code{initial_method} is \code{"buffer"})}
+
+\item{dem}{Digital elevation model (DEM) of the region (only used if
+\code{initial_method} is \code{"valley"})}
+
 \item{capping_method}{The method employed to connect the corridor edge end
 points (i.e. to "cap" the corridor). See \code{\link[=cap_corridor]{cap_corridor()}} for
 the available methods}
@@ -41,6 +50,10 @@ the available methods}
 \item{segments}{Whether to carry out the corridor segmentation}
 
 \item{riverspace}{Whether to carry out the riverspace delineation}
+
+\item{...}{Additional (optional) input arguments for retrieving the DEM
+dataset (see \code{\link[=get_dem]{get_dem()}}). Only relevant if \code{initial_method} is \code{"valley"}
+and \code{dem} is NULL}
 }
 \value{
 A simple feature geometry
diff --git a/man/dem_to_cog.Rd b/man/dem_to_cog.Rd
index d0d5c24..94862aa 100644
--- a/man/dem_to_cog.Rd
+++ b/man/dem_to_cog.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/valley.R
 \name{dem_to_cog}
 \alias{dem_to_cog}
-\title{Write dem to cloud optimized GeoTiff file as specified location}
+\title{Write DEM to cloud optimized GeoTiff file as specified location}
 \usage{
 dem_to_cog(dem, fpath, output_directory = NULL)
 }
@@ -16,6 +16,10 @@ the output directory}
 \item{output_directory}{where file should be written.
 If specified fpath is treated as filename only.}
 }
+\value{
+The input DEM. This function is used for the side-effect of writing
+values to a file.
+}
 \description{
-Write dem to cloud optimized GeoTiff file as specified location
+Write DEM to cloud optimized GeoTiff file as specified location
 }
diff --git a/man/figures/README-example-1.png b/man/figures/README-example-1.png
index 1bcac86..4281d57 100644
Binary files a/man/figures/README-example-1.png and b/man/figures/README-example-1.png differ
diff --git a/man/filter_network.Rd b/man/filter_network.Rd
index 4577d36..872aac6 100644
--- a/man/filter_network.Rd
+++ b/man/filter_network.Rd
@@ -15,5 +15,6 @@ filter_network(network, target)
 A spatial network object
 }
 \description{
-Subset a network keeping the only nodes that intersect a target geometry.
+If subsetting results in multiple disconnected components, we keep the main
+one.
 }
diff --git a/man/get_corridor.Rd b/man/get_corridor.Rd
index 4980706..a422e0f 100644
--- a/man/get_corridor.Rd
+++ b/man/get_corridor.Rd
@@ -9,7 +9,9 @@ get_corridor(
   river_centerline,
   river_surface,
   bbox,
-  initial_method = "buffer",
+  initial_method = "valley",
+  buffer = NULL,
+  dem = NULL,
   capping_method = "direct"
 )
 }
@@ -26,10 +28,15 @@ centerline}
 \item{initial_method}{The method employed to define the initial river
 corridor geometry. See \code{\link[=initial_corridor]{initial_corridor()}} for the available methods}
 
+\item{buffer}{Buffer region to add to the river geometry to setup the initial
+corridor (only used if \code{initial_method} is \code{"buffer"})}
+
+\item{dem}{Digital elevation model (DEM) of the region (only used if
+\code{initial_method} is \code{"valley"})}
+
 \item{capping_method}{The method employed to connect the corridor edge end
 points (i.e. to "cap" the corridor). See \code{\link[=cap_corridor]{cap_corridor()}} for
-the available methods
-TODO: add input arguments for the initialization method}
+the available methods}
 }
 \value{
 A simple feature geometry representing the river corridor
diff --git a/man/get_dem.Rd b/man/get_dem.Rd
index ed17850..92301b5 100644
--- a/man/get_dem.Rd
+++ b/man/get_dem.Rd
@@ -2,26 +2,41 @@
 % Please edit documentation in R/valley.R
 \name{get_dem}
 \alias{get_dem}
-\title{Load dem from a STAC endpoint}
+\title{Access digital elevation model (DEM) for a given region}
 \usage{
-get_dem(bb, resource = "STAC", ...)
+get_dem(
+  bb,
+  dem_source = "STAC",
+  stac_endpoint = NULL,
+  stac_collection = NULL,
+  crs = NULL
+)
 }
 \arguments{
 \item{bb}{A bounding box, provided either as a matrix (rows for "x", "y",
-columns for "min", "max") or as a vector ("xmin", "ymin", "xmax", "ymax")}
+columns for "min", "max") or as a vector ("xmin", "ymin", "xmax", "ymax"),
+in lat/lon coordinates (WGS84 coordinate referece system)}
 
-\item{resource}{from which to source dem. Can be "STAC".
-if "STAC" the parameters the following parameters
-must be supplied as named parameters. If omitted defaults
-will be used.}
+\item{dem_source}{Source of the DEM:
+\itemize{
+\item If "STAC" (default), DEM tiles are searched on a SpatioTemporal Asset
+Catalog (STAC) end point, then accessed and mosaicked to the area of
+interest
+}}
 
-\item{...}{Additional parameters to be passed depending on the resource.
-In case \code{resource = "STAC"}, the arguments \code{endpoint} and
-\code{collection} need to be passed to \code{get_stac_asset_urls()}.}
+\item{stac_endpoint}{URL of the STAC API endpoint (only used if \code{dem_source}
+is \code{"STAC"}). To be provided together with \code{stac_collection}, or leave
+blank to use defaults (see \code{\link{default_stac_dem}})}
+
+\item{stac_collection}{Identifier of the STAC collection to be queried (only
+used if \code{dem_source} is \code{"STAC"}). To be provided together with
+\code{stac_endpoint}, or leave blank to use defaults (see \code{\link{default_stac_dem}})}
+
+\item{crs}{Coordinate reference system (CRS) which to transform the DEM to}
 }
 \value{
-dem
+DEM as a terra SpatRaster object
 }
 \description{
-Load dem from a STAC endpoint
+Access digital elevation model (DEM) for a given region
 }
diff --git a/man/get_osm_city_boundary.Rd b/man/get_osm_city_boundary.Rd
index eee4d3d..43a63a6 100644
--- a/man/get_osm_city_boundary.Rd
+++ b/man/get_osm_city_boundary.Rd
@@ -4,13 +4,13 @@
 \alias{get_osm_city_boundary}
 \title{Get the city boundary from OpenStreetMap}
 \usage{
-get_osm_city_boundary(city_name, bb, crs, multiple = FALSE)
+get_osm_city_boundary(bb, city_name, crs = NULL, multiple = FALSE)
 }
 \arguments{
-\item{city_name}{A character string with the name of the city}
-
 \item{bb}{Bounding box of class \code{bbox}}
 
+\item{city_name}{A character string with the name of the city}
+
 \item{crs}{Coordinate reference system as EPSG code}
 
 \item{multiple}{A logical indicating if multiple city boundaries should be
@@ -27,5 +27,5 @@ The result is filtered by the city name.
 \examples{
 bb <- get_osm_bb("Bucharest")
 crs <- get_utm_zone(bb)
-get_osm_city_boundary("Bucharest", bb, crs)
+get_osm_city_boundary(bb, "Bucharest", crs)
 }
diff --git a/man/get_osm_railways.Rd b/man/get_osm_railways.Rd
index 8ca5eb9..e105580 100644
--- a/man/get_osm_railways.Rd
+++ b/man/get_osm_railways.Rd
@@ -4,7 +4,7 @@
 \alias{get_osm_railways}
 \title{Get OpenStreetMap railways}
 \usage{
-get_osm_railways(bb, crs)
+get_osm_railways(bb, crs = NULL)
 }
 \arguments{
 \item{bb}{Bounding box of class \code{bbox}}
diff --git a/man/get_osm_river.Rd b/man/get_osm_river.Rd
index 9a3d2f2..ab94404 100644
--- a/man/get_osm_river.Rd
+++ b/man/get_osm_river.Rd
@@ -4,13 +4,13 @@
 \alias{get_osm_river}
 \title{Get the river centreline and surface from OpenStreetMap}
 \usage{
-get_osm_river(river_name, bb, crs)
+get_osm_river(bb, river_name, crs = NULL)
 }
 \arguments{
-\item{river_name}{The name of the river}
-
 \item{bb}{Bounding box of class \code{bbox}}
 
+\item{river_name}{The name of the river}
+
 \item{crs}{Coordinate reference system as EPSG code}
 }
 \value{
@@ -22,5 +22,5 @@ Get the river centreline and surface from OpenStreetMap
 \examples{
 bb <- get_osm_bb("Bucharest")
 crs <- get_utm_zone(bb)
-get_osm_river("Dâmbovița", bb, crs)
+get_osm_river(bb, "Dâmbovița", crs)
 }
diff --git a/man/get_osm_streets.Rd b/man/get_osm_streets.Rd
index 4dd27ac..e6aa49f 100644
--- a/man/get_osm_streets.Rd
+++ b/man/get_osm_streets.Rd
@@ -4,7 +4,7 @@
 \alias{get_osm_streets}
 \title{Get OpenStreetMap streets}
 \usage{
-get_osm_streets(bb, crs, highway_values = NULL)
+get_osm_streets(bb, crs = NULL, highway_values = NULL)
 }
 \arguments{
 \item{bb}{Bounding box of class \code{bbox}}
diff --git a/man/get_osmdata.Rd b/man/get_osmdata.Rd
index 24f538a..82fac66 100644
--- a/man/get_osmdata.Rd
+++ b/man/get_osmdata.Rd
@@ -4,18 +4,17 @@
 \alias{get_osmdata}
 \title{Retrieve OpenStreetMap data for a given location}
 \usage{
-get_osmdata(city_name, river_name, crs = NULL, buffer = NULL)
+get_osmdata(bb, city_name, river_name, crs = NULL)
 }
 \arguments{
+\item{bb}{Bounding box defining the area of interest}
+
 \item{city_name}{A character string with the name of the city.}
 
 \item{river_name}{A character string with the name of the river.}
 
 \item{crs}{An integer with the EPSG code for the projection. If no CRS is
 specified, the default is the UTM zone for the city.}
-
-\item{buffer}{A numeric with the buffer distance in meters. By default,
-no buffer is applied}
 }
 \value{
 An list with the retrieved OpenStreetMap data sets for the
@@ -27,5 +26,7 @@ the city boundary, the river centreline and surface, the streets, and the
 railways.
 }
 \examples{
-get_osmdata("Bucharest", "Dambovita", buffer = 2000)
+bb <- get_osm_bb("Bucharest")
+crs <- get_utm_zone(bb)
+get_osmdata(bb, "Bucharest", "Dambovita", crs)
 }
diff --git a/man/get_slope.Rd b/man/get_slope.Rd
index 14c4585..7d04ad3 100644
--- a/man/get_slope.Rd
+++ b/man/get_slope.Rd
@@ -2,8 +2,7 @@
 % Please edit documentation in R/valley.R
 \name{get_slope}
 \alias{get_slope}
-\title{Derive slope as percentage from dem
-This makes use of the terrain function of the terra package}
+\title{Derive slope as percentage from DEM}
 \usage{
 get_slope(dem)
 }
@@ -14,6 +13,5 @@ get_slope(dem)
 raster of derived slope over dem extent
 }
 \description{
-Derive slope as percentage from dem
 This makes use of the terrain function of the terra package
 }
diff --git a/man/get_stac_asset_urls.Rd b/man/get_stac_asset_urls.Rd
index 2774710..99845ab 100644
--- a/man/get_stac_asset_urls.Rd
+++ b/man/get_stac_asset_urls.Rd
@@ -2,28 +2,23 @@
 % Please edit documentation in R/valley.R
 \name{get_stac_asset_urls}
 \alias{get_stac_asset_urls}
-\title{Retrieve asset urls for the intersection of a bounding box with a
-remote STAC endpoint}
+\title{Retrieve the URLs of all the assets intersecting a bbox from a STAC API}
 \usage{
-get_stac_asset_urls(
-  bb,
-  endpoint = "https://earth-search.aws.element84.com/v1",
-  collection = "cop-dem-glo-30"
-)
+get_stac_asset_urls(bb, endpoint, collection)
 }
 \arguments{
 \item{bb}{A bounding box, provided either as a matrix (rows for "x", "y",
-columns for "min", "max") or as a vector ("xmin", "ymin", "xmax", "ymax")}
+columns for "min", "max") or as a vector ("xmin", "ymin", "xmax", "ymax"),
+in lat/lon coordinates (WGS84 coordinate referece system)}
 
-\item{endpoint}{url of (remote) STAC endpoint}
+\item{endpoint}{URL of the STAC API endpoint}
 
 \item{collection}{STAC collection to be queried}
 }
 \value{
-A list of urls for the assets in the collection
-overlapping with the specified bounding box
+A list of URLs for the assets in the collection overlapping with
+the specified bounding box
 }
 \description{
-Retrieve asset urls for the intersection of a bounding box with a
-remote STAC endpoint
+Retrieve the URLs of all the assets intersecting a bbox from a STAC API
 }
diff --git a/man/get_valley.Rd b/man/get_valley.Rd
index c363b7b..9697dd0 100644
--- a/man/get_valley.Rd
+++ b/man/get_valley.Rd
@@ -2,23 +2,26 @@
 % Please edit documentation in R/valley.R
 \name{get_valley}
 \alias{get_valley}
-\title{Create vector/polygon representation of valley from dem and river polygon
-for a crs}
+\title{Extract the river valley from the DEM}
 \usage{
-get_valley(dem, river, crs)
+get_valley(dem, river, bbox = NULL)
 }
 \arguments{
-\item{dem}{of region}
+\item{dem}{Digital elevation model of the region}
 
-\item{river}{vector/polygon representation of river area}
+\item{river}{A simple feature geometry representing the river}
 
-\item{crs}{coordiante reference system to use}
+\item{bbox}{Bounding box defining the extent of the area of interest}
 }
 \value{
-(multi)polygon representation of valley
-area as st_geometry without holes
+River valley as a simple feature geometry
 }
 \description{
-Create vector/polygon representation of valley from dem and river polygon
-for a crs
+The slope of the digital elevation model (DEM) is used as friction (cost)
+surface to compute the cost distance from any grid cell of the raster to
+the river. A characteristic value (default: the mean) of the cost distance
+distribution in a region surrounding the river (default: a buffer region of
+2 km) is then calculated, and used to threshold the cost-distance surface.
+The resulting area is then "polygonized" to obtain the valley boundary as a
+simple feature geometry.
 }
diff --git a/man/get_valley_mask.Rd b/man/get_valley_mask.Rd
deleted file mode 100644
index 151a3d8..0000000
--- a/man/get_valley_mask.Rd
+++ /dev/null
@@ -1,17 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/valley.R
-\name{get_valley_mask}
-\alias{get_valley_mask}
-\title{Select valley pixels from cost distance based on threshold}
-\usage{
-get_valley_mask(cd, thresh)
-}
-\arguments{
-\item{cd}{cost distance raster}
-
-\item{thresh}{threshold cost distance value below which pixels are assuemd
-to belong to the valley}
-}
-\description{
-Select valley pixels from cost distance based on threshold
-}
diff --git a/man/get_valley_polygon.Rd b/man/get_valley_polygon.Rd
index 14bc3ba..236be60 100644
--- a/man/get_valley_polygon.Rd
+++ b/man/get_valley_polygon.Rd
@@ -10,8 +10,8 @@ get_valley_polygon(valley_mask)
 \item{valley_mask}{raster mask of valley pixels}
 }
 \value{
-(multi)polygon representation of valley area
-as st_geometry without holes
+(multi)polygon representation of valley area as a simple feature
+geometry without holes
 }
 \description{
 Create vector/polygon representation of valley without holes from raster mask
diff --git a/man/initial_corridor.Rd b/man/initial_corridor.Rd
index 7dc4f64..7239330 100644
--- a/man/initial_corridor.Rd
+++ b/man/initial_corridor.Rd
@@ -4,19 +4,30 @@
 \alias{initial_corridor}
 \title{Draw the initial geometry of a river corridor.}
 \usage{
-initial_corridor(river, method = "buffer", ..., bbox = NULL)
+initial_corridor(
+  river,
+  method = "valley",
+  buffer = NULL,
+  dem = NULL,
+  bbox = NULL
+)
 }
 \arguments{
 \item{river}{A simple feature geometry representing the river}
 
 \item{method}{The method employed to draw the initial river corridor:
 \itemize{
-\item "buffer" (default): add a fixed buffer region to the river geometry
-(see \code{\link[=river_buffer]{river_buffer()}})
+\item "buffer": add a fixed buffer region to the river geometry (see
+\code{\link[=river_buffer]{river_buffer()}})
+\item "valley" (default): use the river valley boundary, as estimated from the
+provided digital elevation model (DEM, see \code{\link[=river_valley]{river_valley()}})
 }}
 
-\item{...}{Additional arguments required by the function that implements the
-chosen method (see \code{method})}
+\item{buffer}{Buffer region to add to the river geometry (only used if
+\code{initial_method} is \code{"buffer"})}
+
+\item{dem}{Digital elevation model (DEM) of the region (only used if
+\code{initial_method} is \code{"valley"})}
 
 \item{bbox}{Bounding box defining the extent of the area of interest}
 }
diff --git a/man/reproject.Rd b/man/reproject.Rd
index 5192074..d4c0884 100644
--- a/man/reproject.Rd
+++ b/man/reproject.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/valley.R
+% Please edit documentation in R/utils.R
 \name{reproject}
 \alias{reproject}
 \title{Reproject a raster or vector dataset to the specified
diff --git a/man/river_valley.Rd b/man/river_valley.Rd
new file mode 100644
index 0000000..67d884b
--- /dev/null
+++ b/man/river_valley.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/corridor.R
+\name{river_valley}
+\alias{river_valley}
+\title{Draw a corridor as an estimate of the river valley.}
+\usage{
+river_valley(river, dem, bbox = NULL)
+}
+\arguments{
+\item{river}{A simple feature geometry representing the river}
+
+\item{dem}{The DEM of the area as a terra SpatRaster object}
+
+\item{bbox}{Bounding box defining the extent of the area of interest}
+}
+\value{
+A simple feature geometry
+}
+\description{
+The valley is drawn on a digital elevation model (DEM) of the area,
+see \code{\link[=get_valley]{get_valley()}} for details on the implementation.
+}
diff --git a/man/smooth_dem.Rd b/man/smooth_dem.Rd
index 80a9658..99a06d7 100644
--- a/man/smooth_dem.Rd
+++ b/man/smooth_dem.Rd
@@ -9,8 +9,8 @@ smooth_dem(dem, method = "median", window = 5)
 \arguments{
 \item{dem}{raster data of dem}
 
-\item{method}{smoothing function to be used, e.g. "median".
-As accepted by focal}
+\item{method}{smoothing function to be used, e.g. "median",
+as accepted by \code{\link[terra:focal]{terra::focal()}}}
 
 \item{window}{size of smoothing kernel}
 }
diff --git a/tests/testthat/test-corridor.R b/tests/testthat/test-corridor.R
index 1b5e00d..c6b308d 100644
--- a/tests/testthat/test-corridor.R
+++ b/tests/testthat/test-corridor.R
@@ -1,11 +1,26 @@
-test_that("river_buffer parameters can be configured via initial_corridor", {
+test_that("proper parameters must be provided depending on selected method", {
   river <- bucharest_osm$river_centerline
-  actual <- initial_corridor(river, buffer = 1)
-  expected <- river_buffer(river, buffer = 1)
-  expect_setequal(actual, expected)
+
+  # for "buffer" method, we need the "buffer" parameter
+  expect_error(initial_corridor(river, method = "buffer"),
+               "Buffer should be provided if `method` is `'buffer'`")
+  with_mocked_bindings(river_buffer = function(...) NULL, {
+    expect_no_error(initial_corridor(river, method = "buffer", buffer = 42))
+  })
+
+  # for "valley" method, we need the "dem" parameter
+  expect_error(initial_corridor(river, method = "valley"),
+               "DEM should be provided if `method` is `'valley'`")
+  with_mocked_bindings(river_valley = function(...) NULL, {
+    expect_no_error(initial_corridor(river, method = "valley", dem = 42))
+  })
+
+  # inexistent method raise an error
+  expect_error(initial_corridor(river, method = "crisp"),
+               "Unknown method to initialize river corridor: crisp")
 })
 
-test_that("River buffer properly implements a buffer function", {
+test_that("River buffer implements a buffer function", {
   river <- bucharest_osm$river_centerline
   actual <- river_buffer(river, buffer = 0.5)
   expected <- sf::st_buffer(river, 0.5)
diff --git a/tests/testthat/test-network.R b/tests/testthat/test-network.R
index 0160138..c8c6c64 100644
--- a/tests/testthat/test-network.R
+++ b/tests/testthat/test-network.R
@@ -277,6 +277,18 @@ test_that("Filter network properly splits network across adjacent regions", {
   expect_length(nodes_area_2, 2)
 })
 
+test_that("Filter network drops smallest disconnected components", {
+  # p4 is within the area, but it is left out since it remains disconnected
+  # from the main network component
+  area <- sf::st_as_sfc(sf::st_bbox(c(xmin = -1, xmax = 4,
+                                      ymin = 0, ymax = 2)))
+  network_filtered <- filter_network(network, area)
+  edges_area <- sf::st_geometry(sf::st_as_sf(network_filtered, "edges"))
+  nodes_area <- sf::st_geometry(sf::st_as_sf(network_filtered, "nodes"))
+  expect_length(edges_area, 2)
+  expect_length(nodes_area, 3)
+})
+
 test_that("Network setup with real data", {
   edges <- bucharest_osm$streets
   network <- as_network(edges, clean = FALSE, flatten = FALSE)
diff --git a/tests/testthat/test-osmdata.R b/tests/testthat/test-osmdata.R
index 055ecf6..5bf15b5 100644
--- a/tests/testthat/test-osmdata.R
+++ b/tests/testthat/test-osmdata.R
@@ -1,18 +1,10 @@
 test_that("City boundary of Bucharest is correctly retreived", {
   skip_on_ci()
-  bucharest_boundary <- bucharest_osm$boundary |>
-    sf::st_geometry() |>
-    sf::st_transform(4326)
-
   city_name <- "Bucharest"
   bb <- get_osm_bb(city_name)
-  crs <- get_utm_zone(bb)
-
-  bucharest_boundary_osm <- get_osm_city_boundary(city_name, bb, crs) |>
-    sf::st_geometry() |>
-    sf::st_transform(4326)
+  bucharest_boundary <- get_osm_city_boundary(bb, city_name)
 
-  expect_equal(bucharest_boundary, bucharest_boundary_osm, tolerance = 1e-4)
+  expect_equal(as.numeric(bb), as.numeric(sf::st_bbox(bucharest_boundary)))
 })
 
 test_that("City boundary of Paris is returned without error", {
@@ -22,7 +14,8 @@ test_that("City boundary of Paris is returned without error", {
   bb <- get_osm_bb(city_name)
   crs <- get_utm_zone(bb)
 
-  expect_no_error(get_osm_city_boundary(city_name, bb, crs))
+  boundary <- get_osm_city_boundary(bb, city_name, crs)
+  expect_true(!sf::st_is_empty(boundary))
 })
 
 test_that("Wrong city name throws error", {
@@ -32,7 +25,8 @@ test_that("Wrong city name throws error", {
 
 test_that("OSM data for Bucharest is correctly retreived", {
   skip_on_ci()
-  bucharest <- get_osmdata("Bucharest", "Dâmbovița", buffer = 2000)
+  bb <- bucharest_osm$bb
+  bucharest <- get_osmdata(bb, "Bucharest", "Dâmbovița")
 
   expect_length(bucharest, 6)
   expect_true(all(sapply(bucharest, function(x) length(x) >= 1)))
@@ -46,6 +40,6 @@ test_that("Multiple boundaries are correcly retreived", {
   crs <- get_utm_zone(bb)
 
   expect_true(
-    length(get_osm_city_boundary(city_name, bb, crs, multiple = TRUE)) > 1
+    length(get_osm_city_boundary(bb, city_name, crs, multiple = TRUE)) > 1
   )
 })
diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R
index ef75c12..8ac1229 100644
--- a/tests/testthat/test-utils.R
+++ b/tests/testthat/test-utils.R
@@ -89,3 +89,84 @@ test_that("a bbox object does not change class", {
   expect_true(all(as.vector(bbox) == c(0, 1, 2, 3)))
   expect_equal(sf::st_crs(bbox), sf::st_crs(crs))
 })
+
+test_that("buffering a bbox properly enlarge the area of interest", {
+  # bbox in UTM zone 2N
+  x <- c(xmin = 263554, xmax = 736446, ymin = 4987330, ymax = 5654109)
+  bbox_utm2n <- sf::st_bbox(x, crs = "EPSG:32602")
+
+  bbox_buffer_actual <- buffer_bbox(bbox_utm2n, 1000)
+
+  y <- c(x[c("xmin", "ymin")] - 1000, x[c("xmax", "ymax")] + 1000)
+  bbox_buffer_expected <- sf::st_bbox(y, crs = "EPSG:32602")
+
+  expect_equal(bbox_buffer_actual, bbox_buffer_expected)
+})
+
+test_that("buffering a bbox does not change its CRS", {
+  # bbox in WGS 84
+  x <- c(xmin = -174, xmax = -168, ymin = 45, ymax = 51)
+  bbox_wgs84 <- sf::st_bbox(x, crs = "EPSG:4326")
+
+  bbox_buffer <- buffer_bbox(bbox_wgs84, 1000)
+
+  crs_expected <- sf::st_crs(bbox_wgs84)
+  crs_actual <- sf::st_crs(bbox_buffer)
+  expect_equal(crs_actual, crs_expected)
+})
+
+test_that("reproject works with raster data", {
+  # raster in UTM zone 2 (lon between -174 and -168 deg), northern emisphere
+  x <- terra::rast(xmin = -174, xmax = -168, ymin = 45, ymax = 51, res = 1,
+                   vals = 1, crs = "EPSG:4326")
+
+  # reproject with integer (EPSG code)
+  x_repr_int <- reproject(x, 32602)
+
+  # reproject with string
+  x_repr_str <- reproject(x, "EPSG:32602")
+
+  crs_expected <- terra::crs("EPSG:32602")
+  crs_actual_int <- terra::crs(x_repr_int)
+  expect_equal(crs_actual_int, crs_expected)
+  crs_actual_str <- terra::crs(x_repr_str)
+  expect_equal(crs_actual_str, crs_expected)
+})
+
+test_that("reproject works with vector data", {
+  # polygon in UTM zone 2 (lon between -174 and -168 deg), northern emisphere
+  x <- sf::st_linestring(cbind(c(-174, -174, -168, -168, -174),
+                               c(45, 51, 51, 45, 45)))
+  x <- sf::st_polygon(list(x))
+  x <- sf::st_sfc(x, crs = "EPSG:4326")
+
+  # reproject with integer (EPSG code)
+  x_repr_int <- reproject(x, 32602)
+
+  # reproject with string
+  x_repr_str <- reproject(x, "EPSG:32602")
+
+  crs_expected <- sf::st_crs("EPSG:32602")
+  crs_actual_int <- sf::st_crs(x_repr_int)
+  expect_equal(crs_actual_int, crs_expected)
+  crs_actual_str <- sf::st_crs(x_repr_str)
+  expect_equal(crs_actual_str, crs_expected)
+})
+
+test_that("reproject works with bbox", {
+  # bbox in UTM zone 2 (lon between -174 and -168 deg), northern emisphere
+  x <- c(xmin = -174, xmax = -168, ymin = 45, ymax = 51)
+  x <- sf::st_bbox(x, crs = "EPSG:4326")
+
+  # reproject with integer (EPSG code)
+  x_repr_int <- reproject(x, 32602)
+
+  # reproject with string
+  x_repr_str <- reproject(x, "EPSG:32602")
+
+  crs_expected <- sf::st_crs("EPSG:32602")
+  crs_actual_int <- sf::st_crs(x_repr_int)
+  expect_equal(crs_actual_int, crs_expected)
+  crs_actual_str <- sf::st_crs(x_repr_str)
+  expect_equal(crs_actual_str, crs_expected)
+})
diff --git a/tests/testthat/test-valley.R b/tests/testthat/test-valley.R
index 12ca12c..ecf4a83 100644
--- a/tests/testthat/test-valley.R
+++ b/tests/testthat/test-valley.R
@@ -14,37 +14,56 @@ test_that("STAC asset urls are correctly retrieved", {
 
   asset_urls_retrieved <- get_stac_asset_urls(bb, endpoint = ep,
                                               collection = col)
-  asset_urls_retrieved_default <- get_stac_asset_urls(bb)
   expected_asset_urls <- asset_urls
 
   expect_equal(expected_asset_urls, asset_urls_retrieved)
-  expect_equal(expected_asset_urls, asset_urls_retrieved_default)
 })
 
-test_that("raster data are correctly retrieved and merged", {
+test_that("load_raster correctly retrieve and merge local data", {
+
+  write_local_raster <- function(fname, xmin, xmax, ymin, ymax) {
+    rast <- terra::rast(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax,
+                        res = 1, vals = 1, crs = "EPSG:4326")
+    terra::writeRaster(rast, fname)
+  }
+
+  bbox <- sf::st_bbox(c(xmin = 1, xmax = 4, ymin = 1, ymax = 7),
+                      crs = "EPSG:4326")
+  # create local rasters with adjacent bboxes
+  withr::with_file(list("r1.tif" = write_local_raster("r1.tif", 1, 4, 1, 4),
+                        "r2.tif" = write_local_raster("r2.tif", 1, 4, 4, 7)), {
+      rast <- load_raster(bbox, c("r1.tif", "r2.tif"))
+      # all values should be 1
+      expect_true(all(terra::values(rast) == 1))
+      # 2 rasters with 3x3 pixels -> 18 pixels in total
+      expect_length(terra::values(rast), 18)
+      # expect_equal on the two terra::ext objects somehow fails
+      expect_true(terra::ext(rast) == terra::ext(bbox))
+    }
+  )
+})
+
+test_that("load_raster correctly retrieve and merge remote data", {
   skip_on_ci()
 
-  dem <- load_raster(bb, asset_urls) |> CRiSp::reproject(32635)
-  expected_dem <- terra::unwrap(bucharest_dem)
+  dem <- load_raster(bb, asset_urls)
 
-  expect_equal(terra::values(dem), terra::values(expected_dem))
-  expect_equal(terra::crs(dem), terra::crs(expected_dem))
-  expect_true(all.equal(terra::ext(dem), terra::ext(expected_dem),
-                        tolerance = 1e-4))
+  expect_equal(terra::crs(dem), terra::crs("EPSG:4326"))
+  expect_equal(as.vector(terra::ext(dem)), as.vector(terra::ext(bb)),
+               tolerance = 1.e-5)
 })
 
 test_that("valley polygon is correctly constructed", {
   dem <- terra::unwrap(bucharest_dem)
   river <- bucharest_osm$river_surface
-  crs <- "epsg:32635"
 
-  valley <- get_valley(dem, river, crs)
-  expected_valley <- sf::st_read("./testdata/expected_valley.gpkg",
-                                 quiet = TRUE) |>
-    sf::st_as_sfc()
+  valley <- get_valley(dem, river)
+  expected_valley_path <- testthat::test_path("testdata",
+                                              "expected_valley.gpkg")
+  expected_valley <- sf::st_read(expected_valley_path, quiet = TRUE)
 
-  valley <- st_set_precision(valley, 1e-06)
-  expected_valley <- st_set_precision(expected_valley, 1e-06)
+  valley <- sf::st_set_precision(valley, 1e-06)
+  expected_valley <- sf::st_set_precision(expected_valley, 1e-06)
   expect_true(sf::st_equals_exact(valley, expected_valley,
-                                  par = 1.e-4, sparse = FALSE))
+                                  par = 0, sparse = FALSE))
 })
diff --git a/tests/testthat/testdata/expected_valley.gpkg b/tests/testthat/testdata/expected_valley.gpkg
index d285245..34fd9f5 100644
Binary files a/tests/testthat/testdata/expected_valley.gpkg and b/tests/testthat/testdata/expected_valley.gpkg differ
diff --git a/tests/testthat/testdata/generate_valley.R b/tests/testthat/testdata/generate_valley.R
index 79bbaac..cebcb47 100644
--- a/tests/testthat/testdata/generate_valley.R
+++ b/tests/testthat/testdata/generate_valley.R
@@ -1,8 +1,7 @@
 dem <- terra::unwrap(CRiSp::bucharest_dem)
 river <- CRiSp::bucharest_osm$river_surface
-crs <- "epsg:32635"
 
-valley <- CRiSp::get_valley(dem, river, crs)
+valley <- CRiSp::get_valley(dem, river)
 
 filepath <- testthat::test_path("testdata", "expected_valley.gpkg")
 sf::st_write(valley, filepath, append = FALSE, quiet = TRUE)