Skip to content

Commit

Permalink
Merge pull request #41 from CityRiverSpaces/37-segmentation-fn
Browse files Browse the repository at this point in the history
Implement segmentation as part of the package
  • Loading branch information
fnattino authored Dec 19, 2024
2 parents cfe22b9 + 18d53e4 commit dfa958c
Show file tree
Hide file tree
Showing 25 changed files with 669 additions and 105 deletions.
12 changes: 8 additions & 4 deletions DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -13,12 +13,14 @@ Description: CRiSp (City River Spaces) provides tools to automate the
License: Apache License (>= 2)
URL: https://cityriverspaces.github.io/CRiSp/
BugReports: https://github.com/CityRiverSpaces/crisp/issues
Depends:
Depends:
R (>= 2.10)
Imports:
Imports:
dbscan,
dplyr,
lwgeom,
osmdata,
rcoins,
rlang,
rstac,
sf,
Expand All @@ -29,17 +31,19 @@ Imports:
tibble,
tidygraph,
units
Suggests:
Suggests:
ggplot2,
gridExtra,
knitr,
purrr,
rmarkdown,
testthat (>= 3.0.0)
VignetteBuilder:
VignetteBuilder:
knitr
Config/testthat/edition: 3
Encoding: UTF-8
LazyData: true
Roxygen: list(markdown = TRUE)
RoxygenNote: 7.3.2
Remotes:
github::CityRiverSpaces/rcoins
4 changes: 2 additions & 2 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -4,13 +4,12 @@ export(add_weights)
export(as_bbox)
export(as_network)
export(clean_network)
export(corridor)
export(delineate_corridor)
export(delineate_riverspace)
export(delineate_segments)
export(dem_to_cog)
export(flatten_network)
export(get_cd_char)
export(get_corridor)
export(get_cost_distance)
export(get_dem)
export(get_osm_bb)
Expand All @@ -19,6 +18,7 @@ export(get_osm_railways)
export(get_osm_river)
export(get_osm_streets)
export(get_osmdata)
export(get_segments)
export(get_slope)
export(get_stac_asset_urls)
export(get_utm_zone)
Expand Down
25 changes: 1 addition & 24 deletions R/corridor.R
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
#'
#' @return A simple feature geometry representing the river corridor
#' @export
corridor <- function(
get_corridor <- function(
network, river_centerline, river_surface, bbox, initial_method = "buffer",
capping_method = "direct"
) {
Expand Down Expand Up @@ -129,17 +129,6 @@ split_aoi <- function(bbox, river) {
}
}

#' Split a geometry along a (multi)linestring.
#'
#' @param geometry Geometry to split
#' @param line Dividing (multi)linestring
#'
#' @return A simple feature object
split <- function(geometry, line) {
lwgeom::st_split(geometry, line) |>
sf::st_collection_extract()
}

#' Identify the initial edges of the river corridor
#'
#' These are defined by splitting the initial corridor boundary into the
Expand Down Expand Up @@ -208,15 +197,3 @@ cap_corridor <- function(edges, method = "direct", network = NULL) {
}
as_polygon(c(edges, cap_edge_1, cap_edge_2))
}

as_linestring <- function(points) {
points_union <- sf::st_union(points)
sf::st_cast(points_union, "LINESTRING")
}

as_polygon <- function(lines) {
lines_union <- sf::st_union(lines)
sf::st_line_merge(lines_union) |>
sf::st_polygonize() |>
sf::st_collection_extract()
}
28 changes: 15 additions & 13 deletions R/delineate.R
Original file line number Diff line number Diff line change
Expand Up @@ -11,15 +11,18 @@
#' @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
#' @param angle_threshold Only network edges forming angles above this threshold
#' (in degrees) are considered when forming segment edges. See
#' [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
#'
#' @return A simple feature geometry
#' @export
delineate_corridor <- function(
city_name, river_name, crs = NULL, bbox_buffer = NULL,
initial_method = "buffer", capping_method = "direct", segments = FALSE,
riverspace = FALSE
initial_method = "buffer", 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)
Expand All @@ -32,25 +35,24 @@ delineate_corridor <- function(
network <- as_network(network_edges)

# Run the corridor delineation on the spatial network
corridor <- corridor(
corridor <- get_corridor(
network, osm_data$river_centerline, osm_data$river_surface, bbox,
initial_method, capping_method
)

if (segments) delineate_segments()
if (segments) {
# Select the relevant part of the network
buffer_corridor <- 100 # TODO should this be an additional input parameter?
corridor_buffer <- sf::st_buffer(corridor, buffer_corridor)
network_filtered <- filter_network(network, corridor_buffer)

corridor <- get_segments(corridor, network_filtered,
osm_data$river_centerline, angle_threshold)
}
if (riverspace) delineate_riverspace()
return(corridor)
}


#' Delineate segments of a river corridor.
#'
#' @return A simple feature geometry
#' @export
delineate_segments <- function() {
stop("Segmentation not yet implemented.")
}

#' Delinate the riverspace.
#'
#' @return A simple feature geometry
Expand Down
18 changes: 16 additions & 2 deletions R/network.R
Original file line number Diff line number Diff line change
Expand Up @@ -308,6 +308,20 @@ filter_network <- function(network, target) {
tidygraph::filter(sfnetworks::node_intersects(target))
}

strokes <- function() {
stop("`strokes` not yet implemented.")
#' Identify network edges that are intersecting a geometry
#'
#' @param network A spatial network object
#' @param geometry A simple feature geometry
#' @param index Whether to return the indices of the matchin edges or the
#' geometries
#'
#' @return Indices or geometries of the edges intersecting the given geometry
get_intersecting_edges <- function(network, geometry, index = FALSE) {
edges <- sf::st_as_sf(network, "edges")
intersects <- sf::st_intersects(edges, geometry, sparse = FALSE)
if (index) {
return(which(intersects))
} else {
return(edges[intersects, ])
}
}
211 changes: 211 additions & 0 deletions R/segments.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,211 @@
#' Split a river corridor into segments
#'
#' Segments are defined as corridor subregions separated by river-crossing
#' transversal lines that form continuous strokes in the network.
#'
#' @param corridor The river corridor as a simple feature geometry
#' @param network The spatial network to be used for the segmentation
#' @param river_centerline The river centerline as a simple feature geometry
#' @param angle_threshold Only consider angles above this threshold (in degrees)
#' to form continuous strokes in the network. See [`rcoins::stroke()`] for
#' more details.
#'
#' @return Segment polygons as a simple feature geometry
#' @export
get_segments <- function(corridor, network, river_centerline,
angle_threshold = 90) {

# Find river crossings in the network and build continuous strokes from them
crossings <- get_intersecting_edges(network, river_centerline, index = TRUE)
crossing_strokes <- rcoins::stroke(network, from_edge = crossings,
angle_threshold = angle_threshold)

# Clip strokes and select the ones that could be used as segment boundaries
block_edges <- clip_and_filter(crossing_strokes, corridor, river_centerline)

# Split the corridor into candidate segments ("blocks")
blocks <- split(corridor, block_edges)

# Refine the blocks to make sure that all segments touch the river and cross
# the corridor from side to side
refine_segments(blocks, river_centerline, corridor)
}

#' Clip lines to the extent of the corridor, and filter valid segment edges
#'
#' Lines that intersect the river only once and that cross the corridor from
#' side to side are considered valid segment edges. We group valid segment edges
#' that cross the river in nearby locations, and select the shortest line per
#' cluster.
#'
#' @param lines Candidate segment edges as a simple feature geometry
#' @param corridor The river corridor as a simple feature geometry
#' @param river_centerline The river centerline as a simple feature geometry
#'
#' @return Candidate segment edges as a simple feature geometry
#' @importFrom rlang .data
clip_and_filter <- function(lines, corridor, river_centerline) {

# Split corridor along the river centerline to find edges on the two sides
corridor_edges <- split(corridor, river_centerline, boundary = TRUE)

# Clip the lines, keeping the only fragments that intersect the river
lines_clipped <- sf::st_intersection(lines, corridor) |>
sf::st_as_sf() |>
dplyr::filter(sf::st_is(.data$x, c("MULTILINESTRING", "LINESTRING"))) |>
sfheaders::sf_cast("LINESTRING") |>
sf::st_filter(river_centerline, .predicate = sf::st_intersects) |>
sf::st_geometry()

# Select the fragments that cross the river only once and intersect both
# sides of the corridor
river_intersections <- sf::st_intersection(lines_clipped, river_centerline)
# TODO: we could generalize the following to allow for more complex river
# geometries (e.g. for river islands)
intersects_river <- sf::st_is(river_intersections, "POINT")
intersects_side_1 <- sf::st_intersects(lines_clipped, corridor_edges[[1]],
sparse = FALSE)
intersects_side_2 <- sf::st_intersects(lines_clipped, corridor_edges[[2]],
sparse = FALSE)
is_valid <- intersects_river & intersects_side_1 & intersects_side_2
lines_valid <- lines_clipped[is_valid]

# Cluster valid segment edges and select the shortest line per cluster
filter_clusters(lines_valid, river_centerline)
}

#' Cluster the river crossings and select the shortest crossing per cluster
#'
#' Create groups of edges that are crossing the river in nearby locations,
#' using a density-based clustering method (DBSCAN). This is to make sure that
#' edges representing e.g. different lanes of the same street are treated as
#' part of the same crossing. For each cluster, select the shortest edge.
#'
#' @param crossings Crossing edge geometries as a simple feature object
#' @param river The river geometry as a simple feature object
#' @param eps DBSCAN parameter referring to the size (radius) distance of the
#' neighborhood. Should approximate the distance between edges that we want
#' to consider as a single river crossing
#'
#' @return A simple feature geometry including the shortest edge per cluster
filter_clusters <- function(crossings, river, eps = 100) {
intersections <- sf::st_intersection(crossings, river)
# By computing centroids we make sure we only have POINT geometries here
intersections_centroids <- sf::st_centroid(intersections)
intersections_coords <- sf::st_coordinates(intersections_centroids)
# We should not enforce a min mumber of elements - one-element clusters are OK
db <- dbscan::dbscan(intersections_coords, eps = eps, minPts = 1)

crossings_clustered <- sf::st_as_sf(crossings)
crossings_clustered$cluster <- db$cluster
crossings_clustered$length <- sf::st_length(crossings_clustered)
crossings_clustered |>
dplyr::group_by(.data$cluster) |>
dplyr::filter(length == min(length) & !duplicated(length)) |>
sf::st_geometry()
}

#' Refine candidate segments via recursive merging
#'
#' Recursively merge the candidate segments provided ("blocks"), until they all
#' intersect the river centerline and both sides of the corridor.
#'
#' @param blocks Candidate segments as a simple feature geometry
#' @param river_centerline The river centerline as a simple feature geometry
#' @param corridor The river corridor as a simple feature geometry
#'
#' @return Refined corridor segments as a simple feature geometry
refine_segments <- function(blocks, river_centerline, corridor) {

# Split corridor along the river centerline to find edges on the two sides
corridor_edges <- split(corridor, river_centerline, boundary = TRUE)

# Recursively merge blocks until all blocks intersect the river
not_intersect_river <- function(blocks) {
index_instersect_river <- find_intersects(blocks, river_centerline)
index <- seq_along(blocks)
index[!index %in% index_instersect_river]
}
while (TRUE) {
index_not_instersects_river <- not_intersect_river(blocks = blocks)
if (length(index_not_instersects_river) == 0) break
blocks <- merge_blocks(blocks, index_not_instersects_river,
method = "longest-intersection")
}

# Recursively merge blocks until all blocks intersect both edges
not_intersect_both_edges <- function(blocks) {
index_intersects_edge_1 <- find_intersects(blocks, corridor_edges[1])
index_intersects_edge_2 <- find_intersects(blocks, corridor_edges[2])
index <- seq_along(blocks)
index[!(index %in% index_intersects_edge_1 &
index %in% index_intersects_edge_2)]
}
while (TRUE) {
index_not_intersect_both_edges <- not_intersect_both_edges(blocks)
if (length(index_not_intersect_both_edges) == 0) break
blocks <- merge_blocks(blocks, index_not_intersect_both_edges,
method = "smallest")
}
return(blocks)
}

#' Merge a set of blocks to adjacent ones
#'
#' Adjacent blocks are defined as the blocks that are neighbours to the blocks
#' that need to be merged, and that intersect them via a (Multi)LineString. We
#' consider the blocks to merge one by one, from the smallest to the largest,
#' merging them to the other blocks recursively.
#'
#' @param blocks Simple feature geometry representing all the blocks
#' @param to_merge Indices of the blocks to merge
#' @param method Strategy for merging, see [merge_block()]
#'
#' @return Blocks merged to the specified ones as a simple feature geometry
merge_blocks <- function(blocks, to_merge, method = "longest-intersection") {
if (length(to_merge) == 0) {
return(blocks)
}
# Pick the first block to merge, i.e. the smallest one, and the targets
index_smallest <- find_smallest(blocks[to_merge])
index_current <- to_merge[index_smallest]
current <- blocks[index_current]
targets <- blocks[!seq_along(blocks) %in% index_current]
# Keep track of the geometries of the blocks that still need merging: their
# position in the list of blocks might change after merging the current one!
index_others <- to_merge[!seq_along(to_merge) %in% index_smallest]
others <- blocks[index_others]
# Merge current block with one of the targets
merged <- merge_block(targets, current, method = method)
# Determine the new indices of the other blocks that need to be merged
is_equal <- sf::st_equals(merged, others, sparse = TRUE)
index_others <- which(apply(is_equal, any, MARGIN = 1))
# Continue merging other blocks, recursively
merge_blocks(blocks = merged, to_merge = index_others, method = method)
}

#' Merge a block to one of the target geometries
#'
#' @param targets Sequence of target blocks as a simple feature geometry
#' @param block Block to merge as a simple feature geometry
#' @param method Strategy for merging, choose between "smallest" (merge to
#' smallest adjacent block) and "longest-intersection" (merge to block which
#' it shares the longest intersection with)
#'
#' @return Blocks merged to the specified one as a simple feature geometry
merge_block <- function(targets, block, method = "longest-intersection") {
index_adjacent <- find_adjacent(targets, block)
if (method == "longest-intersection") {
intersections <- sf::st_intersection(targets[index_adjacent], block)
index_longest_intersection <- find_longest(intersections)
index_to_merge <- index_adjacent[index_longest_intersection]
} else if (method == "smallest") {
index_smallest <- find_smallest(targets[index_adjacent])
index_to_merge <- index_adjacent[index_smallest]
} else {
stop(sprintf("Method '%s' unknown", method))
}
merged <- sf::st_union(targets[index_to_merge], block)
others <- targets[!seq_along(targets) %in% index_to_merge]
return(c(others, merged))
}
Loading

0 comments on commit dfa958c

Please sign in to comment.