Combine Multilinestring into one line string to generate regular spaced points using sf
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"))
Solution
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
- grid.arrange's relative height argument didn't works well in RMarkdown
- How apply a function to every element of a matrix?
- Polar plot in R using ggplot2
- Find total length of time in multiple values of a dataset that overlaps with a given period of time
- Faster than R's `simplify2array`?
- What is the most simple way to label a data.table frequency chart in R?
- How to show dates in a secondary x-axis while the primary x-axis is numerical?
- online compiler for r-markdown document
- Unecessarily large column widths in HTML tables published by R Quarto when the table header contains special characters
- Bootstrapping a vector of results, by group in R
- Speeding Up Monte Carlo Simulations for Weather-Dependent Claims
- Is there an R class which allows to list 2 elements but shows only the first of them by default when called?
- Convert day of week number to weekday name in R
- How do you use "<<-" (scoping assignment) in R?
- Using reactivePoll to monitor and process file changes
- How to find out which graph is most similar to another graph?
- ggplot2: rectangular background behind axis text in circular bar plot
- How to simulate button click using rvest
- CMIP6, automating by month
- Covert dummy variables to single categorical in R?
- Animated reducing bar chart with gganimate
- Create an object class matrix/array from data frame
- How to add a vector of names to `sort_by()` in base R
- Among many subfolders, find CSVs starting with a given string and within ZIPs or not, and merge them while adding their names in a new column with R
- How to add a greek letter to an axis label to a 3D plot in plotly?
- bigrquery + dplyr error when using functions and sourcing from separate R script
- Shiny App Error: /v1/applications/ 400 - Validation Error Execution halted
- Rupee Unicode replacement not working in R
- Mutate multiple variable to create multiple new variables
- How to use the Stata equivalent of a foreach loop in R with mutate()