Combine Multilinestring into one line string to generate regular spaced points using sf

Question

I'm trying to generate random points on a linestring using SF. However, the data is multilinestring and st_line_sample requires linestring. So I'd like to make the multilinestring into only one linestring. If you look at the code, I'm supposed to generate only one point for the linestring, but it generates one point per linestring. I'm not able to join the lines so it become only one linestring which could be used to generate one point on it.

How is it possible to get only one line string from multilinestring to generate the number of random points on the line?

library(sf)
library(mapview)

chemins.simple %>% 
  st_zm(geom) %>% 
  st_cast('LINESTRING') %>% 
  st_line_sample(n = 1, 
            type = "regular") %>% 
  mapview() + mapview(chemins.simple)

The data is below

chemins.simple = structure(list(Name = "Sentier Parc Maisonneuve", description = NA_character_, 
    timestamp = structure(NA_real_, class = c("POSIXct", "POSIXt"
    ), tzone = ""), begin = structure(NA_real_, class = c("POSIXct", 
    "POSIXt"), tzone = ""), end = structure(NA_real_, class = c("POSIXct", 
    "POSIXt"), tzone = ""), altitudeMode = NA_character_, tessellate = NA_integer_, 
    extrude = NA_integer_, visibility = NA_integer_, drawOrder = NA_integer_, 
    icon = NA_character_, sentier = NA_character_, layer = NA_character_, 
    path = NA_character_, geom = structure(list(structure(list(
        structure(c(299582.847886435, 299567.18907952, 299533.785023203, 
        299515.754560643, 299489.877968861, 299474.41195845, 
        299482.409654204, 299490.968980874, 5047541.45303887, 
        5047515.52023998, 5047476.83372278, 5047450.17574194, 
        5047418.25329369, 5047405.40685391, 5047374.54606886, 
        5047360.74761632, 0, 0, 0, 0, 0, 0, 0, 0), dim = c(8L, 
        3L)), structure(c(299490.968980874, 299497.29763185, 
        299517.115667715, 299552.432591057, 5047360.74761632, 
        5047350.76804681, 5047334.25633307, 5047315.1869463, 
        0, 0, 0, 0), dim = 4:3), structure(c(299552.432591057, 
        299572.471975622, 299580.660505604, 299582.021973701, 
        5047315.1869463, 5047290.4503043, 5047268.99526155, 5047206.73936411, 
        0, 0, 0, 0), dim = 4:3), structure(c(299582.021973701, 
        299588.0164882, 299600.947249195, 299630.602987485, 299665.55385166, 
        299701.470512757, 299727.863346363, 299758.996957146, 
        299802.35340954, 299839.140374969, 299876.649347262, 
        299907.772848518, 299930.504446245, 5047206.73936411, 
        5047178.37871541, 5047157.46467836, 5047127.49377699, 
        5047104.97108792, 5047092.39963862, 5047085.8799339, 
        5047078.22045237, 5047075.4591553, 5047067.43201061, 
        5047049.77064173, 5047028.70621312, 5047003.83116685, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), dim = c(13L, 
        3L)), structure(c(299651.605306308, 299634.64341842, 
        299624.972126109, 299582.847886435, 5047607.50989447, 
        5047604.07051358, 5047596.26302026, 5047541.45303887, 
        0, 0, 0, 0), dim = 4:3), structure(c(299903.07988968, 
        299838.094689635, 299772.477232293, 299750.23396301, 
        299735.638185202, 299715.208156948, 299702.078535598, 
        299685.670961566, 299669.997900071, 299657.600720696, 
        299651.605306308, 5047528.38059123, 5047560.33007056, 
        5047598.91464364, 5047604.02168543, 5047594.94570617, 
        5047587.32852712, 5047587.33920063, 5047592.8052447, 
        5047604.45040378, 5047608.09567279, 5047607.50989447, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), dim = c(11L, 3L)), 
        structure(c(300348.696267149, 300327.322202229, 300298.695288719, 
        300263.694927326, 300246.383397842, 300232.079546476, 
        300194.889733763, 300177.248808883, 300154.270060961, 
        300141.323210388, 300128.611331444, 300115.132109253, 
        300099.688617737, 300072.709345899, 299903.07988968, 
        5047307.57255003, 5047318.94799577, 5047324.78492774, 
        5047343.34999599, 5047361.17527199, 5047377.0898541, 
        5047393.47561699, 5047396.94207427, 5047395.32322611, 
        5047396.42340946, 5047406.92957364, 5047427.66034101, 
        5047443.0760719, 5047456.36482701, 5047528.38059123, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), dim = c(15L, 
        3L)), structure(c(299903.127127851, 299903.320414223, 
        299915.378648026, 299923.777440411, 299930.504446245, 
        5046960.75487254, 5046961.28337721, 5046988.90350206, 
        5047001.25750258, 5047003.83116685, 0, 0, 0, 0, 0), dim = c(5L, 
        3L)), structure(c(299930.504446245, 299949.583123368, 
        299959.545929842, 299980.866628314, 300003.611260665, 
        300023.310726376, 300048.444105637, 300077.234279445, 
        300100.919930205, 300104.481131745, 5047003.83116685, 
        5046971.68800182, 5046943.00547969, 5046920.72166176, 
        5046911.97900684, 5046915.78123032, 5046930.48587239, 
        5046957.6847821, 5046987.79590053, 5046994.35977285, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0), dim = c(10L, 3L)), structure(c(300104.481131745, 
        300119.138934739, 300134.109164134, 300150.529356202, 
        300177.155842816, 300212.510325105, 300224.537076035, 
        300224.346205566, 300237.013058666, 300263.691579113, 
        300284.285763285, 300390.750589685, 300420.305984624, 
        5046994.35977285, 5047024.27306488, 5047045.34739571, 
        5047054.4237419, 5047054.76757684, 5047019.84146242, 
        5047006.38152575, 5046994.74822331, 5046983.65083931, 
        5046961.0006817, 5046941.44510252, 5046881.74703166, 
        5046896.63190805, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0), dim = c(13L, 3L)), structure(c(300457.455775157, 
        300447.44086401, 300426.293881923, 300391.700426951, 
        300375.480661379, 300360.17836955, 300348.696267149, 
        5047223.3217804, 5047233.97932101, 5047244.71822396, 
        5047259.60178559, 5047274.83602678, 5047297.9309186, 
        5047307.57255003, 0, 0, 0, 0, 0, 0, 0), dim = c(7L, 3L
        )), structure(c(300420.305984624, 300449.031711306, 300489.528913656, 
        300498.295918156, 300511.803110218, 300514.373762732, 
        300511.842867013, 300499.106689568, 300493.293484818, 
        300480.192280261, 300468.535068035, 300460.197952614, 
        5046896.63190805, 5046898.79319606, 5046912.94366797, 
        5046931.84221716, 5046948.5561874, 5046974.00273835, 
        5047006.63268926, 5047116.61272606, 5047149.5168537, 
        5047192.05949512, 5047212.74350708, 5047221.64446475, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), dim = c(12L, 3L)), 
        structure(c(300460.197952613, 300457.455775157, 5047221.64446475, 
        5047223.3217804, 0, 0), dim = 2:3)), class = c("XYZ", 
    "MULTILINESTRING", "sfg"))), n_empty = 0L, crs = structure(list(
        input = "EPSG:32188", wkt = "PROJCRS[\"NAD83 / MTM zone 8\",\n    BASEGEOGCRS[\"NAD83\",\n        DATUM[\"North American Datum 1983\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4269]],\n    CONVERSION[\"MTM zone 8\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-73.5,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9999,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",304800,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"easting (E(X))\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"northing (N(Y))\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Canada - Quebec between 75°W and 72°W.; Canada - Ontario - east of 75°W.\"],\n        BBOX[44.98,-75,62.53,-72]],\n    ID[\"EPSG\",32188]]"), class = "crs"), class = c("sfc_MULTILINESTRING", 
    "sfc"), precision = 0, bbox = structure(c(xmin = 299474.41195845, 
    ymin = 5046881.74703166, xmax = 300514.373762732, ymax = 5047608.09567279
    ), class = "bbox"), z_range = structure(c(zmin = 0, zmax = 0
    ), class = "z_range"))), row.names = 1L, sf_column = "geom", agr = structure(c(Name = NA_integer_, 
description = NA_integer_, timestamp = NA_integer_, begin = NA_integer_, 
end = NA_integer_, altitudeMode = NA_integer_, tessellate = NA_integer_, 
extrude = NA_integer_, visibility = NA_integer_, drawOrder = NA_integer_, 
icon = NA_integer_, sentier = NA_integer_, layer = NA_integer_, 
path = NA_integer_), levels = c("constant", "aggregate", "identity"
), class = "factor"), class = c("sf", "data.frame"))

EDIT:

With this data (redrawn in QGIS), it works just fine. However, I don't know why this can't be done easily (to go from the data above, to the data below...)

chemins.simple %>% 
          filter(Name %in% "sentier_velo") %>% 
  st_transform(crs = 32188) %>% 
  st_zm(geom) %>%
  st_sample(size = 25, type = "regular") %>% 
  mapview() + 
  mapview(chemins.simple)

chemins.simple = structure(list(Name = "sentier_velo", description = NA_character_, 
    timestamp = structure(NA_real_, tzone = "", class = c("POSIXct", 
    "POSIXt")), begin = structure(NA_real_, tzone = "", class = c("POSIXct", 
    "POSIXt")), end = structure(NA_real_, tzone = "", class = c("POSIXct", 
    "POSIXt")), altitudeMode = NA_character_, tessellate = NA_integer_, 
    extrude = NA_integer_, visibility = NA_integer_, drawOrder = NA_integer_, 
    icon = NA_character_, sentier = NA_character_, layer = NA_character_, 
    path = NA_character_, geom = structure(list(structure(list(
        structure(c(299903.127127852, 299915.378648026, 299923.777440411, 
        299930.504446245, 299949.583123368, 299959.545929842, 
        299980.866628314, 300003.611260665, 300023.310726376, 
        300048.444105637, 300077.234279445, 300100.919930205, 
        300119.138934739, 300134.109164134, 300150.529356202, 
        300177.155842816, 300224.537076035, 300224.346205566, 
        300284.285763285, 300390.750589685, 300420.305984624, 
        300449.031711306, 300489.528913656, 300498.295918157, 
        300511.803110218, 300514.373762732, 300511.842867013, 
        300499.106689568, 300493.293484818, 300480.192280261, 
        300468.535068035, 300460.197952613, 300457.455775157, 
        300447.44086401, 300426.293881923, 300391.700426951, 
        300375.480661379, 300360.17836955, 300348.696267149, 
        300327.322202229, 300298.695288719, 300263.694927326, 
        300246.383397842, 300232.079546476, 300194.889733763, 
        300177.248808883, 300154.270060961, 300141.323210388, 
        300128.611331444, 300115.132109253, 300099.688617737, 
        300072.709345899, 299903.07988968, 299838.094689635, 
        299772.477232293, 299750.23396301, 299735.638185202, 
        299715.208156948, 299702.078535598, 299685.670961566, 
        299669.997900071, 299657.600720696, 299651.605306308, 
        299634.64341842, 299624.972126109, 299582.847886435, 
        299567.18907952, 299533.785023203, 299515.754560643, 
        299489.877968861, 299474.41195845, 299482.409654204, 
        299490.968980874, 299497.29763185, 299517.115667715, 
        299552.432591057, 299572.471975622, 299580.660505604, 
        299582.021973701, 299588.0164882, 299600.947249195, 299630.602987485, 
        299665.55385166, 299701.470512757, 299758.996957146, 
        299802.35340954, 299839.140374969, 299876.649347262, 
        299907.772848518, 299930.504446245, 5046960.75487254, 
        5046988.90350206, 5047001.25750258, 5047003.83116685, 
        5046971.68800182, 5046943.00547969, 5046920.72166176, 
        5046911.97900684, 5046915.78123032, 5046930.48587239, 
        5046957.6847821, 5046987.79590053, 5047024.27306488, 
        5047045.34739571, 5047054.4237419, 5047054.76757684, 
        5047006.38152575, 5046994.74822331, 5046941.44510252, 
        5046881.74703166, 5046896.63190805, 5046898.79319606, 
        5046912.94366797, 5046931.84221716, 5046948.5561874, 
        5046974.00273835, 5047006.63268926, 5047116.61272606, 
        5047149.5168537, 5047192.05949512, 5047212.74350708, 
        5047221.64446475, 5047223.3217804, 5047233.97932101, 
        5047244.71822396, 5047259.60178558, 5047274.83602677, 
        5047297.9309186, 5047307.57255003, 5047318.94799577, 
        5047324.78492774, 5047343.34999599, 5047361.17527199, 
        5047377.0898541, 5047393.47561699, 5047396.94207427, 
        5047395.3232261, 5047396.42340945, 5047406.92957363, 
        5047427.66034101, 5047443.0760719, 5047456.36482701, 
        5047528.38059123, 5047560.33007056, 5047598.91464364, 
        5047604.02168543, 5047594.94570617, 5047587.32852712, 
        5047587.33920063, 5047592.8052447, 5047604.45040378, 
        5047608.09567278, 5047607.50989447, 5047604.07051358, 
        5047596.26302026, 5047541.45303887, 5047515.52023998, 
        5047476.83372278, 5047450.17574194, 5047418.25329369, 
        5047405.40685391, 5047374.54606886, 5047360.74761632, 
        5047350.76804681, 5047334.25633307, 5047315.1869463, 
        5047290.4503043, 5047268.99526155, 5047206.73936411, 
        5047178.37871541, 5047157.46467836, 5047127.49377699, 
        5047104.97108792, 5047092.39963862, 5047078.22045236, 
        5047075.4591553, 5047067.43201061, 5047049.77064173, 
        5047028.70621312, 5047003.83116685), dim = c(90L, 2L))), class = c("XY", 
    "MULTILINESTRING", "sfg"))), class = c("sfc_MULTILINESTRING", 
    "sfc"), precision = 0, bbox = structure(c(xmin = 299474.41195845, 
    ymin = 5046881.74703166, xmax = 300514.373762732, ymax = 5047608.09567278
    ), class = "bbox"), crs = structure(list(input = "EPSG:32188", 
        wkt = "PROJCRS[\"NAD83 / MTM zone 8\",\n    BASEGEOGCRS[\"NAD83\",\n        DATUM[\"North American Datum 1983\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4269]],\n    CONVERSION[\"MTM zone 8\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-73.5,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9999,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",304800,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"easting (E(X))\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"northing (N(Y))\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Canada - Quebec between 75°W and 72°W.; Canada - Ontario - east of 75°W.\"],\n        BBOX[44.98,-75,62.53,-72]],\n    ID[\"EPSG\",32188]]"), class = "crs"), n_empty = 0L)), row.names = 1L, sf_column = "geom", agr = structure(c(Name = NA_integer_, 
description = NA_integer_, timestamp = NA_integer_, begin = NA_integer_, 
end = NA_integer_, altitudeMode = NA_integer_, tessellate = NA_integer_, 
extrude = NA_integer_, visibility = NA_integer_, drawOrder = NA_integer_, 
icon = NA_integer_, sentier = NA_integer_, layer = NA_integer_, 
path = NA_integer_), levels = c("constant", "aggregate", "identity"
), class = "factor"), class = c("sf", "data.frame"))

Answer 1

Following approach is probably not easily generalizable for other shapes, but happens to work with this particular reprex.

Projected coordinates are first rounded (to a meter) to make sure line endpoints would end up as same nodes in graph network, then linestrings are converted to sfnetworks object (igraph / tidygraph + spatial). From there, we can get the correct segment sequence through Eulerian path (pass through every edge exactly once), convert edges back to sf object and subset line segments with edge indices from the path.

There's also st_line_merge(), and once coordinates are rounded, it does a pretty good job, but it's not able to handle that 3-way junction and returns 2 linestrings, loop + that short section.

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(sfnetworks)
library(tidygraph, warn.conflicts = FALSE)
library(dplyr, warn.conflicts = FALSE)
library(mapview)

line_sfn <-
  chemins.simple %>% 
  st_zm(geom) %>% 
  st_geometry() %>%
  st_cast("LINESTRING") %>% 
  # round coordinates so line endpoints would end up in same graph nodes
  lapply(function(x) round(x, 0)) %>% 
  # back to sfc
  st_sfc(crs = st_crs(chemins.simple)) %>% 
  # directed graph
  as_sfnetwork(directed = TRUE)

line_sfn
#> # A sfnetwork with 13 nodes and 13 edges
#> #
#> # CRS:  EPSG:32188 
#> #
#> # A directed simple graph with 1 component with spatially explicit edges
#> #
#> # A tibble: 13 × 1
#>                  x
#>        <POINT [m]>
#> 1 (299583 5047541)
#> 2 (299491 5047361)
#> 3 (299552 5047315)
#> 4 (299582 5047207)
#> 5 (299931 5047004)
#> 6 (299652 5047608)
#> # ℹ 7 more rows
#> #
#> # A tibble: 13 × 3
#>    from    to                                                                  x
#>   <int> <int>                                                   <LINESTRING [m]>
#> 1     1     2 (299583 5047541, 299567 5047516, 299534 5047477, 299516 5047450, …
#> 2     2     3   (299491 5047361, 299497 5047351, 299517 5047334, 299552 5047315)
#> 3     3     4   (299552 5047315, 299572 5047290, 299581 5047269, 299582 5047207)
#> # ℹ 10 more rows

plot(line_sfn)

as.igraph(line_sfn) %>%  plot()

eulerian_path <- 
  line_sfn %>% 
  # edges back to regular sf
  st_as_sf(active = "edges") %>% 
  # find Eulerian path (pass through every edge exactly once),
  # use edge indices to subset edges in sf object in correct order
  slice(igraph::eulerian_path(line_sfn)$epath %>% as.vector()) %>% 
  # lines to multipoints
  st_cast("MULTIPOINT") %>% 
  # to a single multipoint feature, now correctly ordered
  summarise(do_union = FALSE) %>% 
  # points to linestring
  st_cast("LINESTRING")
eulerian_path
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 299474 ymin: 5046882 xmax: 300514 ymax: 5047608
#> Projected CRS: NAD83 / MTM zone 8
#> # A tibble: 1 × 1
#>                                                                                x
#>                                                                 <LINESTRING [m]>
#> 1 (299903 5046961, 299903 5046961, 299915 5046989, 299924 5047001, 299931 50470…

pt_samples <- st_line_sample(eulerian_path, n = 20)
mapview(eulerian_path) + mapview(pt_samples)

---

-e- to a LINESTRING with just sf

sf::st_line_merge() does bit more than I expected, at least with this reprex we can get a valid result when casting st_line_merge() MULTILINESTRING result to POINTs and then to a LINESTRING:

merged_ <- 
  chemins.simple %>% 
  st_zm(geom) %>% 
  # 1 multiline
  st_geometry() %>%
  # to 13 line features
  st_cast("LINESTRING") %>% 
  lapply(function(x) round(x, 0)) %>% 
  st_sfc(crs = st_crs(chemins.simple)) %>% 
  # to 1 multiline with 13 lines
  st_union() %>% 
  # to 1 multiline with 2 lines
  st_line_merge() %>% 
  # to 96 point features
  st_cast("POINT") %>% 
  # to 1 multipoint
  st_combine() %>% 
  st_cast("LINESTRING")
merged_
#> Geometry set for 1 feature 
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 299474 ymin: 5046882 xmax: 300514 ymax: 5047608
#> Projected CRS: NAD83 / MTM zone 8
#> LINESTRING (299903 5046961, 299915 5046989, 299...
plot(merged_)

^{Created on 2024-04-25 with reprex v2.1.0}