diff --git a/DESCRIPTION b/DESCRIPTION
index 3b64c8d6..86e5f57d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -54,8 +54,7 @@ Imports:
     Rcpp (>= 0.12.1),
     rlang (>= 0.2.2),
     sf (>= 0.6.3),
-    sfheaders,
-    vctrs (>= 0.4.0)
+    sfheaders
 Suggests:
     cyclestreets,
     dodgr (>= 0.2.15),
@@ -68,9 +67,11 @@ Suggests:
     osrm,
     pct,
     rmarkdown (>= 1.10),
-    rsgeo (>= 0.1.6),
+    rsgeo (>= 0.1.6.9000),
     testthat (>= 2.0.0),
     tmap
+Remotes:
+    JosiahParry/rsgeo
 VignetteBuilder: 
     knitr
 Encoding: UTF-8
diff --git a/NAMESPACE b/NAMESPACE
index a5160822..96e821ac 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -23,6 +23,8 @@ S3method(line2pointsn,sf)
 S3method(line2vertices,sf)
 S3method(line_segment,sf)
 S3method(line_segment,sfc_LINESTRING)
+S3method(line_segment1,sf)
+S3method(line_segment1,sfc_LINESTRING)
 S3method(n_vertices,sf)
 S3method(od2line,sf)
 S3method(onewaygeo,sf)
@@ -60,6 +62,7 @@ export(line_breakup)
 export(line_cast)
 export(line_midpoint)
 export(line_segment)
+export(line_segment1)
 export(line_via)
 export(mats2line)
 export(n_vertices)
diff --git a/R/geo_projected.R b/R/geo_projected.R
index a64868cd..82e79ee0 100644
--- a/R/geo_projected.R
+++ b/R/geo_projected.R
@@ -74,6 +74,14 @@ geo_projected.sf <- function(shp, fun, crs = geo_select_aeq(shp), silent = TRUE,
     sf::st_crs(shp) <- sf::st_crs(4326)
   }
   crs_orig <- sf::st_crs(shp)
+  # If the original CRS is already projected, run the fun() on the original:
+  if (!is.na(sf::st_crs(shp)) && !sf::st_is_longlat(shp)) {
+    if (!silent) {
+      message("Running function on original projection")
+    }
+    res <- fun(shp, ...)
+    return(res)
+  }
   shp_projected <- sf::st_transform(shp, crs)
   if (!silent) {
     message(paste0("Running function on a temporary projection: ", crs$proj4string))
@@ -86,20 +94,9 @@ geo_projected.sf <- function(shp, fun, crs = geo_select_aeq(shp), silent = TRUE,
 }
 #' @export
 geo_projected.sfc <- function(shp, fun, crs = geo_select_aeq(shp), silent = TRUE, ...) {
-  # assume it's not projected  (i.e. lat/lon) if there is no CRS
-  if (is.na(sf::st_crs(shp))) {
-    sf::st_crs(shp) <- sf::st_crs(4326)
-  }
-  crs_orig <- sf::st_crs(shp)
-  shp_projected <- sf::st_transform(shp, crs)
-  if (!silent) {
-    message(paste0("Running function on a temporary projection: ", crs$proj4string))
-  }
-  res <- fun(shp_projected, ...)
-  if (grepl("sf", x = class(res)[1])) {
-    res <- sf::st_transform(res, crs_orig)
-  }
-  res
+  shp_sf <- sf::st_as_sf(shp)
+  res <- geo_projected.sf(shp_sf, fun, crs, silent, ...)
+  sf::st_geometry(res)
 }
 #' Perform a buffer operation on a temporary projected CRS
 #'
diff --git a/R/linefuns.R b/R/linefuns.R
index 29df5a5e..854fcab4 100644
--- a/R/linefuns.R
+++ b/R/linefuns.R
@@ -69,9 +69,11 @@ is_linepoint <- function(l) {
 #' }
 line_bearing <- function(l, bidirectional = FALSE) {
   p <- sf::st_geometry(line2points(l))
-  i_s <- 1:length(sf::st_geometry(l)) * 2 - 1
+  i_s <- seq_along(sf::st_geometry(l)) * 2 - 1
+
   bearing_radians <- sapply(i_s, function(i) lwgeom::st_geod_azimuth(p[i:(i + 1)]))
-  bearing <- bearing_radians * 180 / (pi)
+
+  bearing <- bearing_radians * 180 / pi
   if (bidirectional) {
     bearing <- make_bidirectional(bearing)
   }
@@ -149,94 +151,52 @@ line_midpoint <- function(l, tolerance = NULL) {
 #' This function keeps the attributes
 #'
 #' @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.
+#'  If `rsgeo` is not available, the `lwgeom` package is used.
+#' @param debug_mode Should debug messages be printed? Default is FALSE.
 #' @family lines
 #' @export
 #' @examples
-#' l <- routes_fast_sf[2, ]
-#' l_seg2 <- line_segment(l = l, n_segments = 2)
-#' l_seg3 <- line_segment(l = l, n_segments = 3)
-#' l_seg_100 <- line_segment(l = l, segment_length = 100)
-#' l_seg_2000 <- line_segment(l = l, segment_length = 2000)
-#' plot(sf::st_geometry(l_seg2), col = 1:2, lwd = 5)
-#' plot(sf::st_geometry(l_seg3), col = 1:3, lwd = 5)
-#' plot(sf::st_geometry(l_seg_100), col = seq(nrow(l_seg_100)), lwd = 5)
-#' plot(sf::st_geometry(l_seg_2000), col = seq(nrow(l_seg_100)), lwd = 5)
-#' # Multiple lines
 #' 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(l, n_segments = NA, segment_length = NA,
-                         use_rsgeo = rlang::is_installed("rsgeo", version = "0.1.6")) {
+#' plot(l_seg_multi, col = seq_along(l_seg_multi), lwd = 5)
+#' # Test rsgeo implementation:
+#' # rsmulti = line_segment(l, segment_length = 1000, use_rsgeo = TRUE)
+#' # waldo::compare(l_seg_multi, rsmulti)
+line_segment <- function(
+  l, 
+  segment_length = NA,
+  use_rsgeo = NULL,
+  debug_mode = FALSE
+) {
   UseMethod("line_segment")
 }
-
 #' @export
 line_segment.sf <- function(
-    l,
-    n_segments = NA,
-    segment_length = NA,
-    use_rsgeo = rlang::is_installed("rsgeo", version = "0.1.6")
-    ) {
-
-  if (is.na(n_segments) && is.na(segment_length)) {
+  l,
+  segment_length = NA,
+  use_rsgeo = NULL,
+  debug_mode = FALSE
+) {
+  if (is.na(segment_length)) {
     rlang::abort(
-      "`n_segment` or `segment_length` must be set.",
+      "`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)
+  # browser()
   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)
+      if (debug_mode) {
+        message(paste0("Processing row ", i, " of ", n_row_l))
+      }
+      # if( i == 108) {
+      #   browser()
+      # }
+      l_segmented = line_segment1(l[i, ], n_segments = NA, segment_length = segment_length, use_rsgeo)
       res_names <- names(sf::st_drop_geometry(l_segmented))
       # Work-around for https://github.com/ropensci/stplanr/issues/531
       if (i == 1) {
@@ -246,41 +206,98 @@ line_segment.sf <- function(
       l_segmented
       })
     res = bind_sf(res_list)
-    return(res)
-  } else if (is.na(n_segments)) {
-    l_length <- as.numeric(sf::st_length(l))
-    n_segments <- max(round(l_length / segment_length), 1)
-  } else if (n_segments == 1) {
-    return(l)
+  } else {
+    # If there's only one row:
+    res = line_segment1(l, n_segments = NA, segment_length = segment_length, use_rsgeo)
   }
-  from_to_sequence = seq(from = 0, to = 1, length.out = n_segments + 1)
-  suppressWarnings({
-    line_segment_list = lapply(seq(n_segments), function(i) {
-      lwgeom::st_linesubstring(
-        x = l,
-        from = from_to_sequence[i],
-        to = from_to_sequence[i + 1]
-      )
-  })
-  })
-  # first_linestring = lwgeom::st_linesubstring(x = l, from = 0, to = 0.2)
-  res <- bind_sf(line_segment_list)
   res
 }
 
 
 #' @export
 line_segment.sfc_LINESTRING <- function(
-    l,
-    n_segments = NA,
-    segment_length = NA,
-    use_rsgeo = rlang::is_installed("rsgeo", version = "0.1.6")
+  l, 
+  segment_length = NA,
+  use_rsgeo = NULL,
+  debug_mode = FALSE
 ) {
   l <- sf::st_as_sf(l)
-  res <- line_segment(l, n_segments, segment)
+  res = line_segment(l, segment_length = segment_length, use_rsgeo, debug_mode)
   sf::st_geometry(res)
 }
 
+#' Segment a single line, using lwgeom or rsgeo
+#' 
+#' @inheritParams line_segment
+#' @param n_segments The number of segments to divide the line into
+#' @family lines
+#' @export
+#' @examples
+#' l <- routes_fast_sf[2, ]
+#' l_seg2 <- line_segment1(l = l, n_segments = 2)
+#' # Test with rsgeo (must be installed):
+#' # l_seg2_rsgeo = line_segment1(l = l, n_segments = 2, use_rsgeo = TRUE)
+#' # waldo::compare(l_seg2, l_seg2_rsgeo)
+#' l_seg3 <- line_segment1(l = l, n_segments = 3)
+#' l_seg_100 <- line_segment1(l = l, segment_length = 100)
+#' l_seg_1000 <- line_segment1(l = l, segment_length = 1000)
+#' plot(sf::st_geometry(l_seg2), col = 1:2, lwd = 5)
+#' plot(sf::st_geometry(l_seg3), col = 1:3, lwd = 5)
+#' plot(sf::st_geometry(l_seg_100), col = seq(nrow(l_seg_100)), lwd = 5)
+#' plot(sf::st_geometry(l_seg_1000), col = seq(nrow(l_seg_1000)), lwd = 5)
+line_segment1 <- function(
+  l, 
+  n_segments = NA,
+  segment_length = NA,
+  use_rsgeo = NULL
+) {
+  UseMethod("line_segment1")
+}
+
+#' @export
+line_segment1.sf <- function(
+  l, 
+  n_segments = NA,
+  segment_length = NA,
+  use_rsgeo = NULL
+) {
+  if (is.na(n_segments) && is.na(segment_length)) {
+    rlang::abort(
+      "`n_segment` or `segment_length` must be set.",
+      call = rlang::caller_env()
+    )
+  }
+  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) {
+    return(l)
+  }
+
+  # Decide whether to use rsgeo or lwgeom, if not set:
+  if (is.null(use_rsgeo)) {
+    use_rsgeo <- use_rsgeo(l)
+  }
+
+  if (use_rsgeo) {
+    res <- line_segment_rsgeo(l, n_segments)
+  } else {
+    res <- line_segment_lwgeom(l, n_segments)
+  }
+  res
+}
+#' @export
+line_segment1.sfc_LINESTRING <- function(
+  l,
+  n_segments = NA,
+  segment_length = NA,
+  use_rsgeo = NULL
+) {
+  l <- sf::st_as_sf(l)
+  res <- line_segment1(l, n_segments, segment_length = segment_length, use_rsgeo)
+  sf::st_geometry(res)
+}
 
 make_bidirectional <- function(bearing) {
   is_na_bearings <- is.na(bearing)
@@ -306,3 +323,67 @@ bind_sf = function(x) {
   x = sf::st_as_sf(x)
   x
 }
+
+use_rsgeo <- function(shp) {
+  rsgeo_installed <- rlang::is_installed("rsgeo", version = "0.1.6")
+  crs <- sf::st_crs(shp)
+  is_projected <- !crs$IsGeographic
+  # if its NA set false
+  if (is.na(is_projected)) {
+    warning("CRS is NA, assuming projected")
+    is_projected <- TRUE
+  }
+  should_use_rsgeo <- rsgeo_installed & is_projected
+  should_use_rsgeo
+}
+
+line_segment_rsgeo <- function(l, n_segments, segment_length) {
+  crs <- sf::st_crs(l)
+  # 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
+  res_sf <- sf::st_as_sf(res_tbl)
+  res_sf
+}
+
+line_segment_lwgeom <- function(l, n_segments) {
+  from_to_sequence = seq(from = 0, to = 1, length.out = n_segments + 1)
+  suppressWarnings({
+    line_segment_list = lapply(seq(n_segments), function(i) {
+      lwgeom::st_linesubstring(
+        x = l,
+        from = from_to_sequence[i],
+        to = from_to_sequence[i + 1]
+      )
+    })
+  })
+  res <- bind_sf(line_segment_list)
+  res
+}
diff --git a/data-raw/test-line_segment.R b/data-raw/test-line_segment.R
new file mode 100644
index 00000000..f8d71101
--- /dev/null
+++ b/data-raw/test-line_segment.R
@@ -0,0 +1,61 @@
+# Test locally:
+# setwd("~/github/ropensci/stplanr")
+# devtools::load_all()
+
+# Test on lastest version
+remotes::install_dev("stplanr")
+library(stplanr)
+
+rnet_y = sf::read_sf("https://github.com/ropensci/stplanr/releases/download/v1.0.2/rnet_y_ed.geojson")
+rnet_y_projected = sf::st_transform(rnet_y, "EPSG:27700")
+summary(sf::st_length(rnet_y_projected))
+# rnet_y_projected_seg = line_segment(rnet_y_projected, segment_length = 10, debug_mode = TRUE)
+
+rnet_y_projected_seg = line_segment(rnet_y_projected, segment_length = 20, debug_mode = TRUE)
+# 108 is to blame
+failing_line = rnet_y_projected[108,]
+plot(failing_line)
+n_vertices(failing_line)
+sf::st_length(failing_line)
+line_segment1(failing_line, segment_length = 20, use_rsgeo = TRUE)
+line_segment(failing_line, segment_length = 20)
+
+# Try again after removing failing line:
+rnet_y_projected_seg = line_segment(rnet_y_projected[-108,], segment_length = 20, debug_mode = TRUE)
+# Still failed at line 108
+rnet_y_projected_seg = line_segment(rnet_y_projected[1:109,], segment_length = 20, debug_mode = TRUE)
+rnet_y_projected_seg = line_segment(rnet_y_projected[-109,], segment_length = 20, debug_mode = TRUE)
+
+# Was actually 109 that was failing...
+failing_line = rnet_y_projected[109,]
+plot(failing_line)
+n_vertices(failing_line)
+sf::st_length(failing_line)
+line_segment1(failing_line, segment_length = 20, use_rsgeo = TRUE)
+line_segment(failing_line, segment_length = 20)
+line_segment1(failing_line$geometry, segment_length = 20, use_rsgeo = TRUE)
+dput(failing_line$geometry)
+
+library(stplanr)
+
+failing_line = structure(list(structure(c(324957.69921197, 324957.873557727, 
+324959.863123514, 324961.852683597, 324963.822867622, 324969.636546456, 
+324976.718443977, 324996.443964294, 673670.123131518, 673680.139281405, 
+673686.784106964, 673693.428933452, 673698.960855279, 673709.992098018, 
+673722.114520549, 673742.922904206), dim = c(8L, 2L), class = c("XY", 
+"LINESTRING", "sfg"))), class = c("sfc_LINESTRING", "sfc"), precision = 0, bbox = structure(c(xmin = 324957.69921197, 
+ymin = 673670.123131518, xmax = 324996.443964294, ymax = 673742.922904206
+), class = "bbox"), n_empty = 0L)
+sf::st_crs(failing_line) = "EPSG:27700"
+line_segment(failing_line, segment_length = 20)
+
+# Try with rsgeo: 
+geo <- rsgeo::as_rsgeo(sf::st_geometry(failing_line))
+
+
+  # segmentize the line strings
+res <- rsgeo::line_segmentize(geo, n = 4)
+res <- sf::st_cast(sf::st_as_sfc(res), "LINESTRING")
+res
+length(res) # should be 4
+sf::st_length(res)
diff --git a/man/angle_diff.Rd b/man/angle_diff.Rd
index dd818a71..98fc6e9d 100644
--- a/man/angle_diff.Rd
+++ b/man/angle_diff.Rd
@@ -48,6 +48,7 @@ Other lines:
 \code{\link{line_bearing}()},
 \code{\link{line_breakup}()},
 \code{\link{line_midpoint}()},
+\code{\link{line_segment1}()},
 \code{\link{line_segment}()},
 \code{\link{line_via}()},
 \code{\link{mats2line}()},
diff --git a/man/geo_toptail.Rd b/man/geo_toptail.Rd
index af543aee..3f03eff3 100644
--- a/man/geo_toptail.Rd
+++ b/man/geo_toptail.Rd
@@ -47,6 +47,7 @@ Other lines:
 \code{\link{line_bearing}()},
 \code{\link{line_breakup}()},
 \code{\link{line_midpoint}()},
+\code{\link{line_segment1}()},
 \code{\link{line_segment}()},
 \code{\link{line_via}()},
 \code{\link{mats2line}()},
diff --git a/man/is_linepoint.Rd b/man/is_linepoint.Rd
index 9b11a005..7e461837 100644
--- a/man/is_linepoint.Rd
+++ b/man/is_linepoint.Rd
@@ -34,6 +34,7 @@ Other lines:
 \code{\link{line_bearing}()},
 \code{\link{line_breakup}()},
 \code{\link{line_midpoint}()},
+\code{\link{line_segment1}()},
 \code{\link{line_segment}()},
 \code{\link{line_via}()},
 \code{\link{mats2line}()},
diff --git a/man/line2df.Rd b/man/line2df.Rd
index 3e5a70b5..7265a1e9 100644
--- a/man/line2df.Rd
+++ b/man/line2df.Rd
@@ -25,6 +25,7 @@ Other lines:
 \code{\link{line_bearing}()},
 \code{\link{line_breakup}()},
 \code{\link{line_midpoint}()},
+\code{\link{line_segment1}()},
 \code{\link{line_segment}()},
 \code{\link{line_via}()},
 \code{\link{mats2line}()},
diff --git a/man/line2points.Rd b/man/line2points.Rd
index 04cc7c6e..660a3057 100644
--- a/man/line2points.Rd
+++ b/man/line2points.Rd
@@ -50,6 +50,7 @@ Other lines:
 \code{\link{line_bearing}()},
 \code{\link{line_breakup}()},
 \code{\link{line_midpoint}()},
+\code{\link{line_segment1}()},
 \code{\link{line_segment}()},
 \code{\link{line_via}()},
 \code{\link{mats2line}()},
diff --git a/man/line_bearing.Rd b/man/line_bearing.Rd
index 36ddad11..2dba9eee 100644
--- a/man/line_bearing.Rd
+++ b/man/line_bearing.Rd
@@ -39,6 +39,7 @@ Other lines:
 \code{\link{line2points}()},
 \code{\link{line_breakup}()},
 \code{\link{line_midpoint}()},
+\code{\link{line_segment1}()},
 \code{\link{line_segment}()},
 \code{\link{line_via}()},
 \code{\link{mats2line}()},
diff --git a/man/line_breakup.Rd b/man/line_breakup.Rd
index 98f7f41e..17c2efc8 100644
--- a/man/line_breakup.Rd
+++ b/man/line_breakup.Rd
@@ -41,6 +41,7 @@ Other lines:
 \code{\link{line2points}()},
 \code{\link{line_bearing}()},
 \code{\link{line_midpoint}()},
+\code{\link{line_segment1}()},
 \code{\link{line_segment}()},
 \code{\link{line_via}()},
 \code{\link{mats2line}()},
diff --git a/man/line_midpoint.Rd b/man/line_midpoint.Rd
index f8c74bd6..60651a4d 100644
--- a/man/line_midpoint.Rd
+++ b/man/line_midpoint.Rd
@@ -30,6 +30,7 @@ Other lines:
 \code{\link{line2points}()},
 \code{\link{line_bearing}()},
 \code{\link{line_breakup}()},
+\code{\link{line_segment1}()},
 \code{\link{line_segment}()},
 \code{\link{line_via}()},
 \code{\link{mats2line}()},
diff --git a/man/line_segment.Rd b/man/line_segment.Rd
index 341b134d..1ef34894 100644
--- a/man/line_segment.Rd
+++ b/man/line_segment.Rd
@@ -4,40 +4,29 @@
 \alias{line_segment}
 \title{Divide an sf object with LINESTRING geometry into regular segments}
 \usage{
-line_segment(
-  l,
-  n_segments = NA,
-  segment_length = NA,
-  use_rsgeo = rlang::is_installed("rsgeo", version = "0.1.6")
-)
+line_segment(l, segment_length = NA, use_rsgeo = NULL, debug_mode = FALSE)
 }
 \arguments{
 \item{l}{A spatial lines object}
 
-\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.}
+If \code{rsgeo} is available, this faster implementation is used by default.
+If \code{rsgeo} is not available, the \code{lwgeom} package is used.}
+
+\item{debug_mode}{Should debug messages be printed? Default is FALSE.}
 }
 \description{
 This function keeps the attributes
 }
 \examples{
-l <- routes_fast_sf[2, ]
-l_seg2 <- line_segment(l = l, n_segments = 2)
-l_seg3 <- line_segment(l = l, n_segments = 3)
-l_seg_100 <- line_segment(l = l, segment_length = 100)
-l_seg_2000 <- line_segment(l = l, segment_length = 2000)
-plot(sf::st_geometry(l_seg2), col = 1:2, lwd = 5)
-plot(sf::st_geometry(l_seg3), col = 1:3, lwd = 5)
-plot(sf::st_geometry(l_seg_100), col = seq(nrow(l_seg_100)), lwd = 5)
-plot(sf::st_geometry(l_seg_2000), col = seq(nrow(l_seg_100)), lwd = 5)
-# Multiple lines
 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)
+plot(l_seg_multi, col = seq_along(l_seg_multi), lwd = 5)
+# Test rsgeo implementation:
+# rsmulti = line_segment(l, segment_length = 1000, use_rsgeo = TRUE)
+# waldo::compare(l_seg_multi, rsmulti)
 }
 \seealso{
 Other lines: 
@@ -49,6 +38,7 @@ Other lines:
 \code{\link{line_bearing}()},
 \code{\link{line_breakup}()},
 \code{\link{line_midpoint}()},
+\code{\link{line_segment1}()},
 \code{\link{line_via}()},
 \code{\link{mats2line}()},
 \code{\link{n_vertices}()},
diff --git a/man/line_segment1.Rd b/man/line_segment1.Rd
new file mode 100644
index 00000000..80094d97
--- /dev/null
+++ b/man/line_segment1.Rd
@@ -0,0 +1,55 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/linefuns.R
+\name{line_segment1}
+\alias{line_segment1}
+\title{Segment a single line, using lwgeom or rsgeo}
+\usage{
+line_segment1(l, n_segments = NA, segment_length = NA, use_rsgeo = NULL)
+}
+\arguments{
+\item{l}{A spatial lines object}
+
+\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.
+If \code{rsgeo} is not available, the \code{lwgeom} package is used.}
+}
+\description{
+Segment a single line, using lwgeom or rsgeo
+}
+\examples{
+l <- routes_fast_sf[2, ]
+l_seg2 <- line_segment1(l = l, n_segments = 2)
+# Test with rsgeo (must be installed):
+# l_seg2_rsgeo = line_segment1(l = l, n_segments = 2, use_rsgeo = TRUE)
+# waldo::compare(l_seg2, l_seg2_rsgeo)
+l_seg3 <- line_segment1(l = l, n_segments = 3)
+l_seg_100 <- line_segment1(l = l, segment_length = 100)
+l_seg_1000 <- line_segment1(l = l, segment_length = 1000)
+plot(sf::st_geometry(l_seg2), col = 1:2, lwd = 5)
+plot(sf::st_geometry(l_seg3), col = 1:3, lwd = 5)
+plot(sf::st_geometry(l_seg_100), col = seq(nrow(l_seg_100)), lwd = 5)
+plot(sf::st_geometry(l_seg_1000), col = seq(nrow(l_seg_1000)), lwd = 5)
+}
+\seealso{
+Other lines: 
+\code{\link{angle_diff}()},
+\code{\link{geo_toptail}()},
+\code{\link{is_linepoint}()},
+\code{\link{line2df}()},
+\code{\link{line2points}()},
+\code{\link{line_bearing}()},
+\code{\link{line_breakup}()},
+\code{\link{line_midpoint}()},
+\code{\link{line_segment}()},
+\code{\link{line_via}()},
+\code{\link{mats2line}()},
+\code{\link{n_vertices}()},
+\code{\link{onewaygeo}()},
+\code{\link{points2line}()},
+\code{\link{toptail_buff}()}
+}
+\concept{lines}
diff --git a/man/line_via.Rd b/man/line_via.Rd
index 8f4b19ba..a483eca0 100644
--- a/man/line_via.Rd
+++ b/man/line_via.Rd
@@ -46,6 +46,7 @@ Other lines:
 \code{\link{line_bearing}()},
 \code{\link{line_breakup}()},
 \code{\link{line_midpoint}()},
+\code{\link{line_segment1}()},
 \code{\link{line_segment}()},
 \code{\link{mats2line}()},
 \code{\link{n_vertices}()},
diff --git a/man/mats2line.Rd b/man/mats2line.Rd
index 996ed059..2fa174b7 100644
--- a/man/mats2line.Rd
+++ b/man/mats2line.Rd
@@ -37,6 +37,7 @@ Other lines:
 \code{\link{line_bearing}()},
 \code{\link{line_breakup}()},
 \code{\link{line_midpoint}()},
+\code{\link{line_segment1}()},
 \code{\link{line_segment}()},
 \code{\link{line_via}()},
 \code{\link{n_vertices}()},
diff --git a/man/n_vertices.Rd b/man/n_vertices.Rd
index b2a77bae..570a8a58 100644
--- a/man/n_vertices.Rd
+++ b/man/n_vertices.Rd
@@ -27,6 +27,7 @@ Other lines:
 \code{\link{line_bearing}()},
 \code{\link{line_breakup}()},
 \code{\link{line_midpoint}()},
+\code{\link{line_segment1}()},
 \code{\link{line_segment}()},
 \code{\link{line_via}()},
 \code{\link{mats2line}()},
diff --git a/man/onewaygeo.Rd b/man/onewaygeo.Rd
index 9e629298..8ddfd784 100644
--- a/man/onewaygeo.Rd
+++ b/man/onewaygeo.Rd
@@ -40,6 +40,7 @@ Other lines:
 \code{\link{line_bearing}()},
 \code{\link{line_breakup}()},
 \code{\link{line_midpoint}()},
+\code{\link{line_segment1}()},
 \code{\link{line_segment}()},
 \code{\link{line_via}()},
 \code{\link{mats2line}()},
diff --git a/man/points2line.Rd b/man/points2line.Rd
index 30e3918f..52c8f4db 100644
--- a/man/points2line.Rd
+++ b/man/points2line.Rd
@@ -27,6 +27,7 @@ Other lines:
 \code{\link{line_bearing}()},
 \code{\link{line_breakup}()},
 \code{\link{line_midpoint}()},
+\code{\link{line_segment1}()},
 \code{\link{line_segment}()},
 \code{\link{line_via}()},
 \code{\link{mats2line}()},
diff --git a/man/toptail_buff.Rd b/man/toptail_buff.Rd
index 7d8a5440..2102c20e 100644
--- a/man/toptail_buff.Rd
+++ b/man/toptail_buff.Rd
@@ -37,6 +37,7 @@ Other lines:
 \code{\link{line_bearing}()},
 \code{\link{line_breakup}()},
 \code{\link{line_midpoint}()},
+\code{\link{line_segment1}()},
 \code{\link{line_segment}()},
 \code{\link{line_via}()},
 \code{\link{mats2line}()},
diff --git a/vignettes/merging-route-networks.Rmd b/vignettes/merging-route-networks.Rmd
index 215101de..5b1d7f9c 100644
--- a/vignettes/merging-route-networks.Rmd
+++ b/vignettes/merging-route-networks.Rmd
@@ -62,13 +62,29 @@ The initial merged result was as follows (original data on left)
 ```{r}
 funs = list(value = sum, Quietness = mean)
 brks = c(0, 100, 500, 1000, 5000)
-rnet_merged = rnet_merge(rnet_x, rnet_y, dist = 10, segment_length = 20, funs = funs)
+system.time({
+  rnet_merged = rnet_merge(rnet_x, rnet_y, dist = 10, segment_length = 20, funs = funs)
+})
 m1 = tm_shape(rnet_y) + tm_lines("value", palette = "viridis", lwd = 5, breaks = brks)  +
   tm_scale_bar()
 m2 = tm_shape(rnet_merged) + tm_lines("value", palette = "viridis", lwd = 5, breaks = brks)
 tmap_arrange(m1, m2, sync = TRUE, nrow = 1)
 ```
 
+Speed-up the results by transforming to a projected coordinate system:
+
+```{r}
+rnet_x = sf::st_transform(rnet_x, 27700)
+rnet_y = sf::st_transform(rnet_y, 27700)
+```
+
+```{r}
+rnet_y_segmented = line_segment(rnet_y, segment_length = 20, use_rsgeo = TRUE)
+system.time({
+  rnet_merged2 = rnet_merge(rnet_x, rnet_y, dist = 10, segment_length = 20, funs = funs)
+})
+```
+
 Let's check the results:
 
 ```{r}