Add functions from edge delineation prototype

DESCRIPTION

@@ -8,17 +8,25 @@ Authors@R: c(
            comment = c(ORCID = "0000-0003-3286-0139")),
     person("Netherlands eScience Center", , , "[email protected]", role = "fnd")
   )
-Description: CRiSp (City River Spaces) provides tools to automate the 
+Description: CRiSp (City River Spaces) provides tools to automate the
     morphological delineation of riverside urban areas.
 License: Apache License (>= 2)
 URL: https://cityriverspaces.github.io/CRiSp/
 BugReports: https://github.com/CityRiverSpaces/crisp/issues
+Imports: 
+    dplyr,
+    lwgeom,
+    osmdata,
+    sf,
+    sfnetworks,
+    tidygraph
 Suggests: 
     knitr,
     rmarkdown,
     testthat (>= 3.0.0)
+VignetteBuilder: 
+    knitr
 Config/testthat/edition: 3
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.3.1
-VignetteBuilder: knitr

NAMESPACE

@@ -1,2 +1,20 @@
 # Generated by roxygen2: do not edit by hand
 
+export(calc_weights)
+export(clean_network)
+export(create_network)
+export(define_aoi)
+export(delineate_corridor)
+export(get_corridor_edge)
+export(get_latlon)
+export(get_osmdata)
+export(get_osmdata_river_corridor)
+export(get_target_points)
+export(get_vertices)
+export(merge_streets)
+export(not_intersects)
+export(osm_bb)
+export(osmdata_as_sf)
+export(simplify_network)
+export(split_aoi)
+export(trim_network)

R/delineate-corridor.R

@@ -0,0 +1,169 @@
+#' Define an area of interest (AoI) from a bounding box and a coordinate reference system (CRS).
+#'
+#' @param bb A bounding box as a matrix with 4 elements: xmin, ymin, xmax, ymax
+#' @param crs A coordinate reference system as an epsg code, e.g. 4326 for WGS84
+#'
+#' @return An area of interest as a simple feature geometry
+#' @export
+define_aoi <- function(bb, crs) {
+  bbox <- as.vector(bb)
+  names(bbox) <- c("xmin", "ymin", "xmax", "ymax")
+  aoi <- sf::st_bbox(bbox, crs = sf::st_crs(4326)) |>
+    sf::st_as_sfc() |>
+    sf::st_transform(crs)
+}
+
+#' Split the area of interest (AoI) by a river.
+#'
+#' @param aoi Area of interest as a simple feature
+#' @param river River as a simple feature
+#'
+#' @return A simple feature geometry set of two areas of interest
+#' @export
+split_aoi <- function(aoi, river) {
+  areas <- aoi |>
+    lwgeom::st_split(river) |>
+    sf::st_collection_extract()
+}
+
+#' Trim a network to an area of interest (AoI) and a river corridor.
+#'
+#' @param net A network object
+#' @param area An area of interest as a simple feature
+#' @param river_corridor A river corridor as a simple feature
+#'
+#' @return A network object
+#' @export
+trim_network <- function(net, area, river_corridor){
+  net |>
+    sfnetworks::activate("nodes") |>
+    sf::st_filter(area, .predicate = sf::st_intersects) |>
+    sf::st_filter(river_corridor, .predicate = CRiSp::not_intersects)
+}
+
+#' Simplify a street network by removing multiple edges and loops.
+#'
+#' @param net A network object
+#'
+#' @return A simplifed network object
+#' @export
+simplify_network <- function(net) {
+  net|>
+    sfnetworks::activate("edges") |>
+    # TODO incorporate this comment in the function description
+    # reorder the edges so that the shortest is kept
+    dplyr::arrange(sfnetworks::edge_length()) |>
+    dplyr::filter(!tidygraph::edge_is_multiple()) |>
+    dplyr::filter(!tidygraph::edge_is_loop())
+}
+
+#' Clean a street network by subdividing edges and removing pseudo-nodes.
+#'
+#' @param net A network object
+#'
+#' @return A cleaned network object
+#' @export
+clean_network <- function(net) {
+  net |>
+    CRiSp::simplify_network() |>
+    tidygraph::convert(sfnetworks::to_spatial_subdivision) |>
+    tidygraph::convert(sfnetworks::to_spatial_smooth)
+}
+
+#' Determine the end vertices of the initial river corridor.
+#'
+#' Determine the "vertices" of the initial river corridor as the intersections
+#' of the initial river corridor with the AoI boundary.
+#'
+#' @param aoi Area of interest as a simple feature
+#' @param corridor_initial Initial river corridor as a simple feature
+#'
+#' @return A simple feature geometry set of two points
+#' @export
+get_vertices <- function(aoi, corridor_initial) {
+  aoi |>
+    sf::st_boundary() |>
+    sf::st_intersection(corridor_initial) |>
+    # this should consists of two linestring components, determine the endpoints
+    sf::st_cast("POINT")
+}
+
+#' Get start and end points of corridor edge on the network.
+#'
+#' TODO add description
+#'
+#' @param vertices A simple feature geometry set of two points
+#' @param area An area of interest as a simple feature
+#' @param threshold A numeric value
+#'
+#' @return A simple feature geometry set of two points
+#' @export
+get_target_points <- function(vertices, area, threshold = 0.001){
+  vertices |>
+    sf::st_as_sf() |>
+    # TODO incorporate this comment into the function description
+    # keep threshold to check which points intersect the polygons
+    sf::st_filter(area, .predicate = sf::st_is_within_distance, dist = threshold) |>
+    sf::st_geometry()
+}
+
+#' Determine the corridor edge on the network.
+#'
+#' Find the corridor edge on one side of the river by using a shortest path algorithm.
+#'
+#' @param net A network object
+#' @param area An area of interest as a simple feature
+#' @param vertices A simple feature geometry set of two points
+#'
+#' @return A simple feature geometry
+#' @export
+get_corridor_edge <- function(net, area, vertices){
+  target_points <- CRiSp::get_target_points(vertices, area)
+
+  paths <- sfnetworks::st_network_paths(
+    net,
+    from = target_points[1],
+    to = target_points[2],
+    weights = "weight",
+    type = "shortest"
+  )
+
+  edges <- net |> sfnetworks::activate("edges") |> sf::st_geometry()
+  edge_path <- paths |> dplyr::pull("edge_paths") |> unlist()
+  edges[edge_path]
+}
+
+#' Delineate a corridor around a river.
+#'
+#' @param place A place name as a string
+#' @param river A river name as a string
+#' @param crs A coordinate reference system as an epsg code, e.g. 4326 for WGS84
+#'
+#' @return A simple feature geometry
+#' @export
+delineate_corridor <- function(place, river, crs = NULL) {
+  bb <- CRiSp::osm_bb(place)
+  cap <- CRiSp::osmdata_as_sf("place", "city", bb)
+  cap <- cap$osm_multipolygons |>
+    sf::st_geometry()
+  if (!is.null(crs)) cap <- cap |> sf::st_transform(crs)
+
+  highways_value <- c("motorway", "primary", "secondary", "tertiary")
+  net <- CRiSp::osmdata_as_sf("highway", highways_value, bb) |>
+    CRiSp::merge_streets() |>
+    CRiSp::create_network()
+
+  areas <- CRiSp::define_aoi(bb, crs) |> CRiSp::split_aoi(river)
+
+  vertices <- CRiSp::get_vertices(areas, river)
+
+  corridor_edge_1 <- CRiSp::get_corridor_edge(net, areas[1], vertices)
+  corridor_edge_2 <- CRiSp::get_corridor_edge(net, areas[2], vertices)
+  corridor_edges <- sf::st_union(corridor_edge_1, corridor_edge_2)
+
+  corridor <- cap |>
+    lwgeom::st_split(corridor_edges) |>
+    sf::st_collection_extract("POLYGON") |>
+    sf::st_as_sf() |>
+    sf::st_filter(river, .predicate = sf::st_intersects)
+}

R/preprocess.R

@@ -0,0 +1,30 @@
+# TODO if we can avoid adding lines from polygons, we should remove this
+#' Merge streets from polygons and lines
+#'
+#' @param highways A list object containing streets from OpenStreetMap
+#'
+#' @return A data frame with a column named 'highway' containing line strings
+#' @export
+merge_streets <- function(highways) {
+  poly_to_lines <- highways$osm_polygons |>
+    sf::st_cast("LINESTRING")
+  highways_lines <- highways$osm_lines |>
+    dplyr::bind_rows(poly_to_lines)
+}
+
+#' Create a network from a line strings
+#'
+#' @param data A data frame with a column named 'highway' containing line strings
+#' @param crs A coordinate reference system as an epsg code, e.g. 4326 for WGS84
+#'
+#' @return A network object
+#' @export
+create_network <- function(data, crs = NULL) {
+  net <- data |>
+    dplyr::select("highway") |>
+    sfnetworks::as_sfnetwork(directed = FALSE)
+
+  if (!is.null(crs)) {
+    net <- net |> sf::st_transform(crs)
+  }
+}

R/retreive-osm.R

@@ -0,0 +1,76 @@
+#' Retrieve bounding box from OpenStreetMap
+#'
+#' @param name A character string with the name of the place to retrieve the bounding box
+#'
+#' @return A list with the bounding box
+#' @export
+#'
+#' @examples
+#' osm_bb("Bucharest")
+osm_bb <- function(name) {
+  bb <- osmdata::getbb(name)
+}
+
+#' Retrieve OpenStreetMap data as sf object
+#'
+#' @param key A character string with the key to filter the data
+#' @param value A character string with the value to filter the data
+#' @param bb A list with the bounding box
+#'
+#' @return An sf object with the retrieved OpenStreetMap data
+#' @export
+osmdata_as_sf <- function(key, value, bb) {
+  bb |>
+    osmdata::opq() |>
+    osmdata::add_osm_feature(key = key, value = value) |>
+    osmdata::osmdata_sf()
+}
+
+#' Retrieve OpenStreetMap data for a given location
+#'
+#' @param name A character string with the name of the place to retrieve the bounding box
+#' @param key A character string with the key to filter the data
+#' @param value A character string with the value to filter the data
+#'
+#' @return An sf object with the retrieved OpenStreetMap data for the given location
+#' @export
+#'
+#' @examples
+#' get_osmdata("Bucharest", "waterway", "river")
+get_osmdata <- function(name, key, value) {
+  bb <- CRiSp::osm_bb(name)
+  CRiSp::osmdata_as_sf(key, value, bb)
+}
+
+#' Get OpenStreetMap data for a river corridor
+#'
+#' @param name A character string with the name of the place to retrieve the bounding box
+#' @param river_name A character string with the name of the river
+#' @param epsg_code An integer with the EPSG code for the projection
+#' @param buffer A numeric with the buffer distance in meters from the water stream
+#'
+#' @return An sf object with the river corridor
+#' @export
+get_osmdata_river_corridor <- function(name, river_name, epsg_code, buffer) {
+  key = "waterway"
+  value = "river"
+  waterways <- CRiSp::get_osmdata(name, key, value)
+  waterway <- waterways$osm_multilines |>
+    dplyr::filter(name == river_name) |>
+    sf::st_transform(epsg_code) |>
+    sf::st_geometry()
+
+  key = "natural"
+  value = "water"
+  water <- CRiSp::get_osmdata(name, key, value)
+
+  waterbody <- dplyr::bind_rows(water$osm_polygons, water$osm_multipolygons) |>
+    sf::st_transform(epsg_code) |>
+    sf::st_filter(waterway, .predicate = sf::st_intersects) |>
+    sf::st_geometry() |>
+    sf::st_union()
+
+  corridor_initial <- c(waterway, waterbody) |>
+    sf::st_buffer(buffer) |>
+    sf::st_union()
+}

R/utils.R

@@ -0,0 +1,33 @@
+#' Get the latitude and longitude of an sf object
+#'
+#' @param x Object of sf class sf, sfc or sfg
+#'
+#' @return An sf object with the latitude and longitude
+#' @export
+get_latlon <- function(x) {
+  sf::st_transform(x, 4326) |>
+    sf::st_geometry()
+}
+
+#' Check if two sf objects do not intersect
+#'
+#' @param x sf object
+#' @param y sf object
+#'
+#' @return A logical vector
+#' @export
+not_intersects <- function(x, y) {
+  !sf::st_intersects(x, y)
+}
+
+#' Calculate the weights of the edges of a network based on their length
+#'
+#' @param net A network object
+#'
+#' @return A network object with a new column 'weight' containing the length of the edges
+#' @export
+calc_weights <- function(net){
+  net |>
+    sfnetworks::activate("edges") |>
+    dplyr::mutate(weight = sfnetworks::edge_length())
+}