使用 sf 将多线串合并为一个线串以生成规则间隔的点

问题描述 投票：0回答：1

我正在尝试使用 SF 在线串上生成随机点。但是，数据是多线串，并且 

st_line_sample
需要线串。所以我想将多线串变成只有一个线串。如果你看一下代码，我应该只为线串生成一个点，但它为每个线串生成一个点。我无法连接这些线，因此它变成只有一个线串，可用于在其上生成一个点。

如何从多线串中只获取一根线串来生成线上的随机点数？

library(sf)
library(mapview)

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

数据如下

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"))
r spatial mapview
1个回答
0
投票

以下方法可能不容易推广到其他形状，但恰好适用于这个特定的表示。

投影坐标首先四舍五入（到一米），以确保线端点最终成为图形网络中的相同节点，然后将线串转换为 

sfnetworks
对象（
igraph
/ 
tidygraph
+ 空间）。从那里，我们可以通过欧拉路径获得正确的线段序列（通过每条边一次），将边转换回 
sf
对象以及具有路径中的边索引的子集线段。

还有 

st_line_merge()
，一旦坐标被舍入，它的工作就非常好，但它无法处理那个 3 路连接点并返回 2 个线串，循环+那个短部分。

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)

创建于 2024-04-24，使用 reprex v2.1.0

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.