r - 连接线串

Question

如何连接两个线串？

可以在线路周围放置一个非常细的缓冲区，然后像这样连接它们:

one_line <- lines %>% 
  st_buffer(0.05) %>% 
  st_intersection() %>% 
  st_union() %>% 
  st_cast('LINESTRING')

这有两个问题:

a) 下面是我的数据的一个非常小的子集，其中包含一个这样的断开的线段 - 如果我在一小部分上使用上述方法，它会形成一个完整的多边形，当转换为线串时只会形成一个非常窄的环

b) 如果我使用整个数据集，它可以工作，但会在我的原始行周围缓冲区的近似距离处创建行。见下图:

蓝色和红色是边缘线，黑色是原始线。

我想简单地将它们取平均，但是当我将这 2 条线转换为坐标 (st_coordinates()) 时，生成的表格具有不同的长度并且没有按顺序排列。

我环顾四周，但没有真正找到任何有用的答案。

这是几何数据的输入:

lines <- structure(list(structure(list(structure(c(2880, 2880.92, 2881.72, 
                                                2882.47, 2883.17, 2883.84, 2884.5, 2894.05, 2894.69, 2895.29393034826, 
                                                340255.362641509, 340257.22, 340259.03, 340260.85, 340262.69, 
                                                340264.55, 340266.4, 340293.7, 340295.61, 340297.500995024), .Dim = c(10L, 
                                                                                                                      2L)), structure(c(2907.22402724177, 2914.21353757771, 340330.886392736, 
                                                                                                                                        340350.2), .Dim = c(2L, 2L))), class = c("XY", "MULTILINESTRING", 
                                                                                                                                                                                 "sfg")), structure(c(2895.3, 2896.82, 2897.26, 2897.72, 2907.2, 
                                                                                                                                                                                                      340297.52, 340302.26, 340303.58, 340304.89, 340330.82), .Dim = c(5L, 
                                                                                                                                                                                                                                                                       2L), class = c("XY", "LINESTRING", "sfg"))), n_empty = 0L, crs = structure(list(
                                                                                                                                                                                                                                                                         input = "EPSG:31256", wkt = "PROJCRS[\"MGI / Austria GK East\",\n    BASEGEOGCRS[\"MGI\",\n        DATUM[\"Militar-Geographische Institut\",\n            ELLIPSOID[\"Bessel 1841\",6377397.155,299.1528128,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4312]],\n    CONVERSION[\"Austria Gauss-Kruger East\",\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\",16.3333333333333,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",1,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",-5000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"northing (X)\",north,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"easting (Y)\",east,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Austria east of 14°50'E of Greenwich (32°30'E of Ferro).\"],\n        BBOX[46.56,14.83,49.02,17.17]],\n    ID[\"EPSG\",31256]]"), class = "crs"), idx = structure(c(1, 
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             2, 1, 1), .Dim = c(2L, 2L)), class = c("sfc_GEOMETRY", "sfc"), precision = 0, bbox = structure(c(xmin = 2880, 
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              ymin = 340255.362641509, xmax = 2914.21353757771, ymax = 340350.2
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             ), class = "bbox"), classes = c("MULTILINESTRING", "LINESTRING"
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             ))

Answer 1

对于上面的示例，您可以转换为 MULTIPOINT，然后联合，然后转换为 LINESTRING。

``` r
library(tidyverse) #overkill, but easier
library(sf)
library(patchwork) #to plot side-by-side

# load data from above 
# lines <- 

single_line <- lines %>%
  st_as_sf() %>%
  st_cast('MULTIPOINT') %>%
  st_union() %>%
  st_cast('LINESTRING')

head(single_line)
#> Geometry set for 1 feature 
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 2880 ymin: 340255.4 xmax: 2914.214 ymax: 340350.2
#> Projected CRS: MGI / Austria GK East
#> LINESTRING (2880 340255.4, 2880.92 340257.2, 28...

p1 <- ggplot() +
  geom_sf(data = st_as_sf(lines), col = c('red', 'blue')) + ggtitle('lines')
          
p2 <- ggplot() + geom_sf(data = single_line, col = 'black') + ggtitle('lines cast & unioned')

p3 <- p1 + p2

p3

^{由 reprex package 创建于 2022-03-08 (v0.3.0)}