sf 方法从通道中描绘中心线

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

我希望用Voronoi方法来划定河道中心线。我遇到了多个来源，但由于已退役的软件包而未能复制该工作（

rgeos
）。

即https://gist.github.com/FloodHydrology/829d3846f0722610ff69d896ad65af3e 识别栅格地图上的线性特征并使用 R

返回线性形状对象 
library(sf)
library(tibble)
library(sf)

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)

st_voronoi(channel) %>% 
  st_collection_extract()

我已经做到了： 但不知道如何使用 Voronoi 对象的边缘来获取中心点，然后获取中心线。

tidyverse gis r-sf
1个回答
0
投票

这里有一个简化的方法，只是为了说明一般工作流程，但对于实际任务，评论中建议的 

centerline
包可能会为您节省大量时间，并且您也可能会得到更好的结果。

Voronoi 曲面细分对中心线检测的适用性很大程度上取决于形状的顶点密度和顶点分布，因此在这里我只是添加更多点，在第一个示例中使用 

st_line_sample()
（均匀分布，适合可视化，可能会错过一些关键的点）点在弯曲处），并且在第二个示例中使用 
st_segmentize()
（可能在中心线末端引入了一些伪影）。

由于我们追求的是 Vornoi 边缘并且并不真正关心多边形，因此我们可以选择 

st_voronoi(... , bOnlyEdges = TRUE)

从生成的边缘中，我们只对 

channel
多边形内的边缘感兴趣，因此我们可以使用空间过滤器来删除大部分（但可能不是全部）杂波。从那里，我们可以将其视为一个图形问题 - 中心线是最长的测地线，即图形直径。

tidygraph
/ 
igarph
和 
sf
整齐地组合在 
sfnetworks
包中，这让我们可以方便地将剩余 Voronoi 边的 
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)
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.