Menu
I am trying to split a large line shape by pairs of points along the line. Based on previous questions asked mainly by @mbcaradima here, here, here and here, I have put together some code which works for single points, but does not for multiple points.
See reproducible example below:
# libraries
library(tidyverse)
library(sf)
library(rgdal)
library(sp)
river <- structure(list(FWK_Code = structure(1L, .Label = "1_F461", class = "factor"),
FWK_Name = structure(1L, .Label = "Glonn von Odelzhausen bis Mündung in die Amper", class = "factor"),
Achse_SWK = structure(NA_integer_, .Label = character(0), class = "factor"),
GEWKZ_11 = structure(1L, .Label = "1648", class = "factor"),
NAM11_SEG = structure(1L, .Label = "Glonn", class = "factor"),
GWO11_SEG = structure(1L, .Label = "1330", class = "factor"),
FKT11_SEG = structure(1L, .Label = "2403", class = "factor"),
ZUST_WWA = structure(1L, .Label = "M", class = "factor"),
ZUST_REG = structure(1L, .Label = "OB", class = "factor"),
FF_LAND = structure(1L, .Label = "DEBY", class = "factor"),
GEW_TYP_K = structure(1L, .Label = "F2.2", class = "factor"),
GEW_TYP_L = structure(1L, .Label = "Kleine Flüsse des Alpenvorlandes", class = "factor"),
HMWB = structure(1L, .Label = "HMWB", class = "factor"),
length_m = 27210, code_vis = structure(1L, .Label = "1_F461", class = "factor"),
geometry = structure(list(structure(c(4454482.82717063, 4454477.69,
4454464.67, 4454447.61, 4454437.73, 4454426.51, 4454409.9,
4454391.94, 4454373.98, 4454359.62, 4454344.8, 4454335.82,
4454317.21, 4454313.45, 4454312.77, 4454309.34, 4454300.59,
4454291.38, 4454284.64, 4454268.93, 4454253.67, 4454234.36,
4454211.47, 4454186.33, 4454158.04, 4454140.98, 4454133.83,
4454125.77, 4454120.53, 4454117.07, 4454097.87, 4454073.83,
4454011.25, 4453877.38, 4453783.85, 4453742.74, 4453742.32,
4453727.02, 4453716.29, 4453710.08, 4453698.84, 4453688.08,
4453682.1, 4453664.37, 4453653.96, 4453641.58, 4453618.07,
4453588.82, 4453561.8, 4453548.29, 4453537.32, 4453529.3,
4453523.71, 4453522.55, 4453517.06, 4453514.53, 4453514.95,
4453523.81, 4453525.5, 4453522.97, 4453510.31, 4453495.53,
4453463.67, 4453452.21, 4453432.45, 4453423.25, 4453395.29,
4453350.73, 4453259.43, 4453181.41, 4453163.58, 4453158.68,
4453154.58, 4453146.63, 4453134.59, 4453098.66, 4453000.7,
4452909.68, 4452627.91, 4452620.51, 4452562.96, 4452508.19,
4452448.4, 4452400.29, 4452376.72, 4452311.09, 4452225.2,
4452203.98, 4452176.98, 4452166.66, 4452135.14, 4452101.77,
4452052.37, 4452014.28, 4451988.32, 4451899.24, 4451893.85,
4451871.58, 4451817.72, 4451785.2, 4451761.21, 4451731.24,
4451712.01, 4451693.25, 4451682.14, 4451657.13, 4451636.02,
4451630.51, 4451622.89, 4451615.39, 4451591.93, 4451581.61,
4451560.51, 4451537.06, 4451523.92, 4451511.73, 4451483.78,
4451476.25, 4451471.86, 4451469.3, 4451438.56, 4451384.62,
4451284.5, 4451259.71, 4451221.42, 4451119.19, 4451110.21,
4451107.87, 4451099.18, 4451091.28, 4451046.03, 4450945.32,
4450884.24, 4450846.94, 4450819.81, 4450790.1, 4450758.62,
4450714.54, 4450696.54, 4450669.05, 4450656.15, 4450619.21,
4450587.43, 4450583.75, 4450567.55, 4450507.21, 4450483.95,
4450465.2, 4450375.63, 4450354.81, 4450328, 4450252.04, 4450237.08,
4450220.65, 4450187.34, 4450152.66, 4450131.38, 4450124.66,
4450110.43, 4450085.02, 4450056.18, 4450033.86, 4450018.41,
4449976.87, 4449966.73, 4449964.81, 4449954.93, 4449938.61,
4449917.72, 4449902.2, 4449887.72, 4449884.86, 4449859, 4449853.42,
4449849, 4449844.8, 4449803.42, 4449777.71, 4449752.73, 4449744.70318071,
5358599.95603184, 5358596.63, 5358590.34, 5358582.71, 5358576.42,
5358568.79, 5358556.22, 5358541.86, 5358527.04, 5358514.47,
5358500.11, 5358493.82, 5358477.42, 5358473.28, 5358472.52,
5358468.68, 5358456.59, 5358443.54, 5358432.76, 5358410.32,
5358389.67, 5358362.73, 5358328.61, 5358290.45, 5358242.86,
5358215.93, 5358210, 5358204.52, 5358201.19, 5358198.95,
5358187.14, 5358181.75, 5358176.36, 5358173.05, 5358175.29,
5358176.19, 5358176.2, 5358172.41, 5358157.12, 5358147.5,
5358123.63, 5358104, 5358093.08, 5358067.75, 5358055.79,
5358041.58, 5358016, 5357984.17, 5357955.89, 5357933.94,
5357912.84, 5357896.38, 5357876.81, 5357872.74, 5357846.99,
5357828.42, 5357811.11, 5357733.87, 5357711.07, 5357689.55,
5357645.65, 5357603.02, 5357523.83, 5357493.98, 5357445.35,
5357424.83, 5357379.53, 5357334.4, 5357245.54, 5357170.31,
5357153.66, 5357149.82, 5357146.6, 5357140.36, 5357130.91,
5357110.43, 5357055.29, 5357002.09, 5356939.82, 5356937.7,
5356920.74, 5356906.82, 5356890.51, 5356880.48, 5356875.56,
5356856.4, 5356828.51, 5356822.82, 5356812.72, 5356808.62,
5356788.72, 5356764.09, 5356726.68, 5356706.5, 5356696.34,
5356684.31, 5356683.6, 5356679.74, 5356656.68, 5356636.56,
5356613.5, 5356583.86, 5356560.41, 5356531.79, 5356511.3,
5356451.12, 5356398.12, 5356387.4, 5356381.23, 5356374.67,
5356363.88, 5356362, 5356357.31, 5356355.43, 5356355.43,
5356357.31, 5356362.31, 5356363.6, 5356364.35, 5356364.6,
5356367.63, 5356373.26, 5356382.3, 5356382.62, 5356379.07,
5356353.88, 5356351.95, 5356351.44, 5356348.99, 5356347.05,
5356335.95, 5356307.71, 5356291.47, 5356283.7, 5356278.35,
5356278.2, 5356289.69, 5356311.72, 5356321.69, 5356330.54,
5356332.74, 5356335.88, 5356329.43, 5356328.02, 5356321.84,
5356294.15, 5356285.09, 5356277.79, 5356250.55, 5356243.77,
5356235.05, 5356200.51, 5356193.71, 5356184.67, 5356163.04,
5356144.16, 5356135.23, 5356134.24, 5356132.14, 5356131.8,
5356136.6, 5356145.19, 5356151.37, 5356173.49, 5356177.49,
5356178.04, 5356178.04, 5356176.53, 5356173.65, 5356168.12,
5356161.8, 5356159.97, 5356138.51, 5356131.75, 5356126.26,
5356121.43, 5356074.91, 5356049.16, 5356016.55, 5356007.68509448
), .Dim = c(180L, 2L), class = c("XY", "LINESTRING", "sfg"
))), class = c("sfc_LINESTRING", "sfc"), precision = 0, bbox = structure(c(xmin = 4449744.70318071,
ymin = 5356007.68509448, xmax = 4454482.82717063, ymax = 5358599.95603184
), class = "bbox"), crs = structure(list(epsg = NA_integer_,
proj4string = "+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +units=m +no_defs"), class = "crs"), n_empty = 0L)), row.names = 0L, class = c("sf",
"data.frame"), sf_column = "geometry", agr = structure(c(FWK_Code = NA_integer_,
FWK_Name = NA_integer_, Achse_SWK = NA_integer_, GEWKZ_11 = NA_integer_,
NAM11_SEG = NA_integer_, GWO11_SEG = NA_integer_, FKT11_SEG = NA_integer_,
ZUST_WWA = NA_integer_, ZUST_REG = NA_integer_, FF_LAND = NA_integer_,
GEW_TYP_K = NA_integer_, GEW_TYP_L = NA_integer_, HMWB = NA_integer_,
length_m = NA_integer_, code_vis = NA_integer_), class = "factor", .Label = c("constant",
"aggregate", "identity")))
sites <- st_line_sample(river, density = 1/1500) %>%
st_jitter(500) %>%
st_sf() %>%
st_cast('POINT') %>%
mutate(group = c(1,1,2,2)) %>%
mutate(type = rep(c("start","stop"),2))
mapview(sites) + river
Using only a single site, the following code works. As the spapped points do not exactly intersect the line, I have added a small buffer. This is hardly elegant, but will do for my situation.
nrst = st_nearest_points(sites, river)
on_line = st_cast(nrst, "POINT")[2]
buf <- st_buffer(on_line,0.1)
parts = st_collection_extract(lwgeom::st_split(river$geometry, buf),"LINESTRING")
mapview(parts[1],color="red",lwd=5) + mapview(parts[2],color="green",lwd=5) + mapview(parts[3],lwd=5)
Using multiple sites produces Error in st_split.sfc(river$geometry, buf_all) :
length(y) == 1 is not TRUE
on_line_all = st_cast(nrst, "POINT")
buf_all <- st_buffer(on_line_all,0.1)
parts_all = st_collection_extract(lwgeom::st_split(river$geometry, buf_all),"LINESTRING")
Unfortunately, in this case, I'm not sure what "y" is. Also, it would be ideal, if all line segments between pairs of start and stop points could be saved to a separate object automatically, rather than having to select them by hand in e.g. QGIS.
Any suggestions will be well appreciated!
808
Turns out there is an error in the documentation for lwgeom::st_split. It states that y is the
object split with (blade); if y contains more than one feature geometry, the geometries are st_combined
However, they are not st_combined when supplying a multi-feature sfc (as in your case).
To get your st_split to work, you'd need to combine them beforehand:
on_line_all = st_cast(nrst, "POINT")
buf_all <- st_combine(st_buffer(on_line_all,0.1))
parts_all = st_collection_extract(lwgeom::st_split(river$geometry, buf_all),"LINESTRING")
parts_all = st_as_sf(
data.frame(
id = 1:length(parts_all),
geometry = parts_all
)
)
mapview(parts_all, burst = "id")
Not sure if it is what you are looking for, but this function breaks a LINESTRING or a POLYGON into segments defined by the coordinates of the geometry itself:
stdh_cast_substring <- function(x, to = "MULTILINESTRING") {
ggg <- st_geometry(x)
if (!unique(st_geometry_type(ggg)) %in% c("POLYGON", "LINESTRING")) {
stop("Input should be LINESTRING or POLYGON")
}
for (k in 1:length(st_geometry(ggg))) {
sub <- ggg[k]
geom <- lapply(
1:(length(st_coordinates(sub)[, 1]) - 1),
function(i)
rbind(
as.numeric(st_coordinates(sub)[i, 1:2]),
as.numeric(st_coordinates(sub)[i + 1, 1:2])
)
) %>%
st_multilinestring() %>%
st_sfc()
if (k == 1) {
endgeom <- geom
}
else {
endgeom <- rbind(endgeom, geom)
}
}
endgeom <- endgeom %>% st_sfc(crs = st_crs(x))
if (class(x)[1] == "sf") {
endgeom <- st_set_geometry(x, endgeom)
}
if (to == "LINESTRING") {
endgeom <- endgeom %>% st_cast("LINESTRING")
}
return(endgeom)
}
You may want to check the post where I described the function Cast a line to subsegments in R.
32
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
