Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement valley as initial method for corridor delineation #59

Merged
merged 27 commits into from
Jan 20, 2025
Merged
Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
27 commits
Select commit Hold shift + click to select a range
1499d5f
implement valley initial method
fnattino Jan 13, 2025
f9901c8
fix example
fnattino Jan 13, 2025
91d64a0
update man page
fnattino Jan 13, 2025
a38da70
fix syntax
fnattino Jan 14, 2025
1a1b463
add and improve tests
fnattino Jan 14, 2025
021bedf
add dependency
fnattino Jan 14, 2025
f870c83
fix for failing tests on windows
fnattino Jan 14, 2025
e5c8659
rollback changes
fnattino Jan 14, 2025
2b38343
make consistent use of the same library for reprojection
fnattino Jan 14, 2025
6725526
drop function
fnattino Jan 15, 2025
8f12f84
add return statement, fixing r cmd check
fnattino Jan 15, 2025
34b50e1
fix issue with bboxes intersecting a single dem tile
fnattino Jan 15, 2025
845745d
update documentation and exports
fnattino Jan 15, 2025
0b0329e
Merge remote-tracking branch 'origin/main' into 52-corridor-valley-fn
fnattino Jan 15, 2025
fdd6f2a
fix default for bbox buffer
fnattino Jan 15, 2025
75b6e71
reintroduce initial corridor being built on both river centerline and…
fnattino Jan 16, 2025
82a4b94
update default initial_method
fnattino Jan 17, 2025
80d7708
update code snippet in readme
fnattino Jan 17, 2025
39cbcc1
update man pages
fnattino Jan 17, 2025
529d1dc
Update exmaple plot in README
cforgaci Jan 17, 2025
7245ae7
remove default for initial buffer
fnattino Jan 17, 2025
c1c6c70
Update R/corridor.R
fnattino Jan 17, 2025
1bd21e0
Merge branch '52-corridor-valley-fn' of github.com:CityRiverSpaces/CR…
fnattino Jan 17, 2025
28cc75b
fix line length
fnattino Jan 17, 2025
645d4af
update man pages
fnattino Jan 17, 2025
94a650c
Update R/valley.R
fnattino Jan 20, 2025
5438189
Update R/valley.R
fnattino Jan 20, 2025
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
49 changes: 39 additions & 10 deletions R/corridor.R
Original file line number Diff line number Diff line change
Expand Up @@ -7,22 +7,24 @@
#' @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"
buffer = 1000, dem = NULL, capping_method = "direct"
) {

# Draw the initial corridor geometry within the area of interest
river <- c(river_centerline, river_surface)
fnattino marked this conversation as resolved.
Show resolved Hide resolved
# TODO: remove hardcoded buffer value here
corridor_init <- initial_corridor(river, initial_method, bbox, buffer = 1000)
corridor_init <- initial_corridor(river_surface, method = initial_method,
buffer = buffer, dem = dem, bbox = bbox)

# Find the corridor end points
end_points <- corridor_end_points(river_centerline, bbox)
Expand Down Expand Up @@ -54,15 +56,28 @@ get_corridor <- function(
#' @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`)
#' - "valley": use the river valley boundary, as estimated from the provided
fnattino marked this conversation as resolved.
Show resolved Hide resolved
#' 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 = "buffer", 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)
Expand All @@ -87,6 +102,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,
Expand Down
33 changes: 25 additions & 8 deletions R/delineate.R
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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 = "buffer", initial_buffer = 1000, 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)
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) {
Expand Down
8 changes: 7 additions & 1 deletion R/network.R
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
68 changes: 32 additions & 36 deletions R/osmdata.R
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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.
Expand All @@ -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()
}

Expand All @@ -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.")
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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))
}
Expand All @@ -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")
}
Expand All @@ -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)
}
Expand All @@ -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)
}
Loading
Loading