From c1726e8aa08fe02e24700d08cba1c0886790171e Mon Sep 17 00:00:00 2001 From: Josiah Parry Date: Sat, 30 Sep 2023 13:06:19 -0400 Subject: [PATCH] Conditionally use rsgeo::line_segmentize() if available (#542) * conditionally use rsgeo::line_segmentize() if available * Document use_rsgeo * Document --------- Co-authored-by: robinlovelace --- DESCRIPTION | 3 +- NAMESPACE | 2 ++ NEWS.md | 5 +++ R/linefuns.R | 86 +++++++++++++++++++++++++++++++++++++++++---- man/line_segment.Rd | 10 +++++- 5 files changed, 97 insertions(+), 9 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 4cdd972c..1fd96c0a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: stplanr Title: Sustainable Transport Planning -Version: 1.1.2 +Version: 1.1.2.9000 Authors@R: c( person("Robin", "Lovelace", , "rob00x@gmail.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-5679-6536")), @@ -65,6 +65,7 @@ Suggests: osrm, pct, rmarkdown (>= 1.10), + rsgeo (>= 0.1.6), testthat (>= 2.0.0), tmap VignetteBuilder: diff --git a/NAMESPACE b/NAMESPACE index 00efc842..a5160822 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -21,6 +21,8 @@ S3method(line2points,sfc) S3method(line2points,sfg) S3method(line2pointsn,sf) S3method(line2vertices,sf) +S3method(line_segment,sf) +S3method(line_segment,sfc_LINESTRING) S3method(n_vertices,sf) S3method(od2line,sf) S3method(onewaygeo,sf) diff --git a/NEWS.md b/NEWS.md index d620751e..ae633bed 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,8 @@ +# stplanr (development version) + +- `line_segment()` now will use `rsgeo::line_segmentize()` if available and only when the input geometry is not in geographic coordinates or does not have a CRS set as it uses euclidean distance + - `line_segment()` becomes an S3 generic which now has methods for `sf` and `sfc` class objects + # stplanr 1.1.2 (2023-09) - Export S3 methods diff --git a/R/linefuns.R b/R/linefuns.R index 8ca0b159..29df5a5e 100644 --- a/R/linefuns.R +++ b/R/linefuns.R @@ -151,6 +151,8 @@ line_midpoint <- function(l, tolerance = NULL) { #' @inheritParams line2df #' @param n_segments The number of segments to divide the line into #' @param segment_length The approximate length of segments in the output (overides n_segments if set) +#' @param use_rsgeo Should the `rsgeo` package be used? +#' If `rsgeo` is available, this faster implementation is used by default. #' @family lines #' @export #' @examples @@ -167,13 +169,72 @@ line_midpoint <- function(l, tolerance = NULL) { #' l <- routes_fast_sf[2:4, ] #' l_seg_multi = line_segment(l, segment_length = 1000) #' plot(sf::st_geometry(l_seg_multi), col = seq(nrow(l_seg_100)), lwd = 5) -line_segment <- function( +line_segment <- function(l, n_segments = NA, segment_length = NA, + use_rsgeo = rlang::is_installed("rsgeo", version = "0.1.6")) { + UseMethod("line_segment") +} + +#' @export +line_segment.sf <- function( l, n_segments = NA, - segment_length = NA + segment_length = NA, + use_rsgeo = rlang::is_installed("rsgeo", version = "0.1.6") ) { + + if (is.na(n_segments) && is.na(segment_length)) { + rlang::abort( + "`n_segment` or `segment_length` must be set.", + call = rlang::caller_env() + ) + } + + # if rsgeo is available use it + if (use_rsgeo) { + # if CRS is NA then we can continue or if IsGeographic is NA + crs <- sf::st_crs(l) + is_geographic <- crs$IsGeographic + + # if its NA set FALSE, if not keep + is_geographic <- ifelse(is.na(is_geographic), FALSE, is_geographic) + + if (is.na(crs) || !is_geographic) { + # extract geometry and convert to rsgeo + geo <- rsgeo::as_rsgeo(sf::st_geometry(l)) + + # if n_segments is missing it needs to be calculated + if (is.na(n_segments)) { + l_length <- rsgeo::length_euclidean(geo) + n_segments <- max(round(l_length / segment_length), 1) + } + + # segmentize the line strings + res <- rsgeo::line_segmentize(geo, n_segments) + + # make them into sfc_LINESTRING + res <- sf::st_cast(sf::st_as_sfc(res), "LINESTRING") + + # give them them CRS + res <- sf::st_set_crs(res, crs) + + # calculate the number of original geometries + n <- length(geo) + # create index ids to grab rows from + ids <- rep.int(1:n, rep(n_segments, n)) + + # index the original sf object + res_tbl <- sf::st_drop_geometry(l)[ids,] + + # assign the geometry column + res_tbl[[attr(l, "sf_column")]] <- res + + # convert to sf and return + return(sf::st_as_sf(res_tbl)) + } + } + n_row_l = nrow(l) - if(n_row_l > 1) { + if (n_row_l > 1) { res_list = pbapply::pblapply(seq(n_row_l), function(i) { l_segmented = line_segment(l[i, ], n_segments, segment_length) res_names <- names(sf::st_drop_geometry(l_segmented)) @@ -186,12 +247,10 @@ line_segment <- function( }) res = bind_sf(res_list) return(res) - } - if (is.na(n_segments)) { + } else if (is.na(n_segments)) { l_length <- as.numeric(sf::st_length(l)) n_segments <- max(round(l_length / segment_length), 1) - } - if(n_segments == 1) { + } else if (n_segments == 1) { return(l) } from_to_sequence = seq(from = 0, to = 1, length.out = n_segments + 1) @@ -210,6 +269,19 @@ line_segment <- function( } +#' @export +line_segment.sfc_LINESTRING <- function( + l, + n_segments = NA, + segment_length = NA, + use_rsgeo = rlang::is_installed("rsgeo", version = "0.1.6") +) { + l <- sf::st_as_sf(l) + res <- line_segment(l, n_segments, segment) + sf::st_geometry(res) +} + + make_bidirectional <- function(bearing) { is_na_bearings <- is.na(bearing) non_na_bearings <- bearing[!is_na_bearings] diff --git a/man/line_segment.Rd b/man/line_segment.Rd index 1943cdaa..341b134d 100644 --- a/man/line_segment.Rd +++ b/man/line_segment.Rd @@ -4,7 +4,12 @@ \alias{line_segment} \title{Divide an sf object with LINESTRING geometry into regular segments} \usage{ -line_segment(l, n_segments = NA, segment_length = NA) +line_segment( + l, + n_segments = NA, + segment_length = NA, + use_rsgeo = rlang::is_installed("rsgeo", version = "0.1.6") +) } \arguments{ \item{l}{A spatial lines object} @@ -12,6 +17,9 @@ line_segment(l, n_segments = NA, segment_length = NA) \item{n_segments}{The number of segments to divide the line into} \item{segment_length}{The approximate length of segments in the output (overides n_segments if set)} + +\item{use_rsgeo}{Should the \code{rsgeo} package be used? +If \code{rsgeo} is available, this faster implementation is used by default.} } \description{ This function keeps the attributes