Sf-методы для отделения осевой линии от канала

Вот упрощенный подход, просто чтобы проиллюстрировать общий рабочий процесс, хотя для реальной задачи пакет centerline, предложенный в комментариях, может сэкономить вам массу времени, и вы, вероятно, также получите лучшие результаты.

Пригодность тесселяции Вороного для обнаружения осевой линии сильно зависит от плотности вершин формы и распределения вершин, поэтому здесь я просто добавляю больше точек, в первом примере с st_line_sample() (равномерное распределение, удобное для визуализации, может пропустить некоторые критические точки на изгибах ), а во втором примере — с st_segmentize() (вероятно, вносит некоторые артефакты на концах центральной линии).

Поскольку нас интересуют края Ворного и нас не особо волнуют полигоны, мы можем пойти st_voronoi(... , bOnlyEdges = TRUE)

Из получившихся ребер нас интересуют только те, которые находятся внутри многоугольника channel, поэтому мы можем использовать пространственный фильтр, чтобы удалить большую часть (но, скорее всего, не весь) беспорядка. Отсюда мы можем думать об этом как о задаче о графе: осевая линия — это самая длинная геодезическая, диаметр графа.

tidygraph/igarph и sf аккуратно собраны в пакете sfnetworks, который позволяет нам удобно конвертировать наш sf объект оставшихся ребер Вороного в граф. И после обработки графа мы можем преобразовать его обратно в sf.

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(dplyr)
library(sfnetworks)
library(tidygraph)
library(ggplot2)

# uniform point distribution along polygon line 
# for simplification & visualization
edge_points <-
  channel %>% 
  st_cast("LINESTRING") %>% 
  st_line_sample(density = 0.02) %>% 
  st_union()
edge_points
#> Geometry set for 1 feature 
#> Geometry type: MULTIPOINT
#> Dimension:     XY
#> Bounding box:  xmin: 656774.6 ymin: 5552713 xmax: 657642.6 ymax: 5553838
#> CRS:           NA
#> MULTIPOINT ((656930.3 5553621), (656880.9 55536...

# we can set bOnlyEdges to get edges instead of polygons
vor_edges <- 
  st_voronoi(edge_points, bOnlyEdges = TRUE) %>% 
  # from single MULTILINESTRING to a set of LINESTRINGs
  st_cast("LINESTRING")
vor_edges
#> Geometry set for 183 features 
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 655649.7 ymin: 5551588 xmax: 658767.4 ymax: 5554963
#> CRS:           NA
#> First 5 geometries:
#> LINESTRING (655649.7 5554589, 656873.7 5553711)
#> LINESTRING (656873.7 5553711, 656873.7 5553711)
#> LINESTRING (656873.7 5553711, 656873.7 5553711)
#> LINESTRING (656873.7 5553711, 655649.7 5553857)
#> LINESTRING (656873.7 5553711, 656873.6 5553711)

# keep only edges within channel
edge_within <- 
  vor_edges %>% 
  # for convenience of st_filter we need to convert sfc to sf
  st_sf() %>% 
  st_filter(channel, .predicate = st_within) 
edge_within
#> Simple feature collection with 68 features and 0 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 656873.6 ymin: 5552816 xmax: 657543.1 ymax: 5553731
#> CRS:           NA
#> First 10 features:
#>                          geometry
#> 1  LINESTRING (656873.7 555371...
#> 2  LINESTRING (656873.7 555371...
#> 3  LINESTRING (656873.7 555371...
#> 4  LINESTRING (656996.3 555330...
#> 5  LINESTRING (656997.3 555329...
#> 6  LINESTRING (656873.8 555371...
#> 7  LINESTRING (656873.8 555371...
#> 8  LINESTRING (656996.6 555330...
#> 9  LINESTRING (656900.8 555371...
#> 10 LINESTRING (656900.8 555371...

# to remove edges that are not part of the centre line, 
# yet are too short to intersect with channel polygon outline, 
# we can convert the remaining list of edges to a graph, 
# keep only nodes (and edges) that are part of a graph diameter (longest geodesic),  
# remove pseudo nodes (nodes with 2 incident edges) with to_spatial_smooth morpher and 
# convert resulting graph to a single LINESTRING
channel_cline <- 
  edge_within %>% 
  as_sfnetwork(directed = FALSE) %>% 
  activate(nodes) %>% 
  # get_diameter returns a list node tbl indices, so we can use it with slice
  slice(igraph::get_diameter(.) %>% as.vector()) %>% 
  convert(to_spatial_smooth) %>%
  st_as_sf(active = "edges")
channel_cline
#> Simple feature collection with 1 feature and 3 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 656873.6 ymin: 5552831 xmax: 657543.1 ymax: 5553725
#> CRS:           NA
#> # A tibble: 1 × 4
#>    from    to .tidygraph_edge_index                                     geometry
#>   <int> <int> <list>                                                <LINESTRING>
#> 1     1     2 <int [63]>            (656873.6 5553711, 656873.7 5553711, 656873…


# visualize previous steps ---------------------------------------------------
# cropped variant for plot
vor_crop <- st_crop(vor_edges, st_buffer(channel, 100))

p1 <- ggplot() + 
  geom_sf(data = channel, fill = "lightblue") +
  geom_sf(data = edge_points) +
  # Calssify edges (within / not within channel polygon)
  geom_sf(data = vor_crop, aes(color = st_within(vor_crop, channel, sparse = FALSE)), show.legend = FALSE) +
  ggtitle("Vor. edges from \nunif. sampled points") +
  theme_void() 

p2 <- ggplot() +
  geom_sf(data = channel, fill = "lightblue") +
  geom_sf(data = edge_within) +
  ggtitle("Vor. edges within \nchannel poly") +
  theme_void() +
  theme(legend.position = "bottom", legend.title.position = "top")

p3 <- ggplot() +
  geom_sf(data = channel, fill = "lightblue") +
  geom_sf(data = channel_cline) +
  ggtitle("Longest geodesic\n") +
  theme_void() +
  theme(legend.position = "bottom", legend.title.position = "top")

patchwork::wrap_plots(p1, p2, p3)

Поскольку st_voronoi работает с многоугольником, мы можем просто добавить некоторую плотность вершин к краю этого канала с помощью st_segmentize(), поэтому в одном конвейере это может выглядеть примерно так:

cline_sfc <- 
  channel %>% 
  # we can add more points to the polygon outline and pass it directly to st_voronoi
  st_segmentize(dfMaxLength = 10) %>%
  st_voronoi(bOnlyEdges = TRUE) %>% 
  st_cast("LINESTRING") %>% 
  st_sf() %>% 
  st_filter(channel, .predicate = st_within) %>%
  as_sfnetwork(directed = FALSE) %>% 
  activate(nodes) %>% 
  slice(igraph::get_diameter(.) %>% as.vector()) %>% 
  # there's a good chance that there are few artefacts at the beginning and end,
  # naive but perhaps good enough approach would be just to remove few nodes
  # from both ends of the diameter list:
  # slice(igraph::get_diameter(.) %>% as.vector() %>% head(-2) %>% tail(-2)) %>% 
  convert(to_spatial_smooth) %>%
  st_as_sf(active = "edges") %>%
  st_geometry() 
  
cline_sfc
#> Geometry set for 1 feature 
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 656873.7 ymin: 5552826 xmax: 657553.8 ymax: 5553729
#> CRS:           NA
#> LINESTRING (656873.7 5553711, 656873.7 5553711,...

ggplot() +
  geom_sf(data = channel, fill = "lightblue") +
  geom_sf(data = cline_sfc) +
  ggtitle("cline_sfc") +
  theme_void() 

Форма ввода:

channel <- st_as_sf(tribble(
  ~x, ~y,
  656873.7, 5553711,
  656993.1, 5553738,
  656993.6, 5553738,
  656994.1, 5553738,
  656994.6, 5553738,
  656995.1, 5553738,
  656995.7, 5553738,
  656996.2, 5553738,
  656996.7, 5553738,
  656997.2, 5553738,
  656997.7, 5553738,
  656998.2, 5553738,
  656998.7, 5553738,
  656999.2, 5553738,
  656999.7, 5553737,
  657000.2, 5553737,
  657000.6, 5553737,
  657001.0, 5553737,
  657001.5, 5553736,
  657001.9, 5553736,
  657002.3, 5553736,
  657002.6, 5553735,
  657003.0, 5553735,
  657003.3, 5553734,
  657003.6, 5553734,
  657126.8, 5553550,
  657127.0, 5553550,
  657127.3, 5553549,
  657127.5, 5553549,
  657127.7, 5553548,
  657127.9, 5553548,
  657128.1, 5553548,
  657128.2, 5553547,
  657128.3, 5553546,
  657128.4, 5553546,
  657128.4, 5553545,
  657128.4, 5553545,
  657129.8, 5553466,
  657129.8, 5553465,
  657129.7, 5553465,
  657129.7, 5553464,
  657129.6, 5553464,
  657129.5, 5553463,
  657129.4, 5553463,
  657129.2, 5553462,
  657129.0, 5553462,
  657128.8, 5553461,
  657128.6, 5553461,
  657128.3, 5553460,
  657128.1, 5553460,
  657127.8, 5553459,
  657127.4, 5553459,
  657127.1, 5553459,
  657126.7, 5553458,
  657126.4, 5553458,
  657126.0, 5553458,
  656995.2, 5553355,
  656971.7, 5553303,
  656975.2, 5553281,
  657061.9, 5553229,
  657233.8, 5553179,
  657234.3, 5553179,
  657234.8, 5553179,
  657235.3, 5553179,
  657235.9, 5553178,
  657236.3, 5553178,
  657312.0, 5553130,
  657312.5, 5553130,
  657312.9, 5553130,
  657313.3, 5553129,
  657313.7, 5553129,
  657314.1, 5553128,
  657314.5, 5553128,
  657314.8, 5553128,
  657315.1, 5553127,
  657315.4, 5553127,
  657315.6, 5553126,
  657315.9, 5553126,
  657316.1, 5553125,
  657316.2, 5553125,
  657406.2, 5552837,
  657467.5, 5552813,
  657534.9, 5552855,
  657535.4, 5552855,
  657535.9, 5552855,
  657536.3, 5552855,
  657536.8, 5552855,
  657537.3, 5552856,
  657537.8, 5552856,
  657538.3, 5552856,
  657538.8, 5552856,
  657539.4, 5552856,
  657539.9, 5552856,
  657540.4, 5552856,
  657540.9, 5552856,
  657541.4, 5552856,
  657542.0, 5552856,
  657542.5, 5552856,
  657543.0, 5552856),
  coords = c("x", "y")) %>% 
  st_combine() %>% 
  st_cast("LINESTRING") 

channel <- st_buffer(channel, 100)