我有2个SHP文件,一个有许多字符串(铁路),另一个有许多多边形(美国的县).这两个文件之间根本没有交集.我的目标是将线捕捉到最近的多边形边界,以便进行交叉并计算每个多边形中线的长度.

我try 使用sf包中的st_snap来实现这一目标,但没有得到预期的结果. 这让我在一些简单的数据上判断这个函数,并得到了这些结果:

library(dplyr)
library(tidyverse)
library(sf)
library(units)
library(ggplot2)

poly = st_polygon(list(cbind(c(0, 0, 1, 1, 0), c(0, 1, 1, 0, 0))))
lines = st_multilinestring(list(
 cbind(c(0, 1), c(1.01, 1.05)),
 cbind(c(0, 1), c(-0.01, -.05)),
 cbind(c(1.025, 1.05, 1.025), c(1.05, .5, -.05))
))
t <- st_intersection(poly, lines)
st_is_empty(t)
# [1] TRUE

ggplot()+
    geom_sf(data = lines, col = "green")+
    geom_sf(data = poly, col = "red", fill = NA)

Lines (green) & polygon (red)

拍照后我的预期结果是:

但我得到的是:

tolerance=0.5

snapped1 = st_snap(lines, poly, tolerance=0.5)
print(snapped1)
# MULTILINESTRING ((0 1, 1 1), (0 0, 1 0), (1 1, 1.05 0.5, 1 0))

ggplot()+
    geom_sf(data = lines, col = "green")+
    geom_sf(data = poly, col = "red", fill = NA)+
    geom_sf(data = snapped1, col = "blue", alpha = 0.5)

Snpped lines in blue

tolerance=1.005

snapped2 = st_snap(lines, poly, tolerance=1.005)
print(snapped2)
# MULTILINESTRING ((1 1, 0 1, 0 0, 1 0), (0 1, 1 1, 1 0, 0 0), (0 1, 1 1, 1.05 0.5, 1 0, 0 0))

其中包含预期结果:(0 1, 1 1, 1 0, 0 0),但也包含不需要的行.

ggplot()+
    geom_sf(data = lines, col = "green")+
    geom_sf(data = poly, col = "red", fill = NA)+
    geom_sf(data = snapped2, col = "blue", alpha = 0.5)

Snpped lines in blue

tolerance=1.1

snapped3 = st_snap(lines, poly, tolerance=1.1)
print(snapped3)
# MULTILINESTRING ((1 1, 0 1, 0 0, 1 0), (0 1, 1 1, 1 0, 0 0), (1 1, 1 0, 0 0, 0 1))

其中还包含预期结果:(0 1, 1 1, 1 0, 0 0),但也包含不需要的行.

ggplot()+
    geom_sf(data = lines, col = "green")+
    geom_sf(data = poly, col = "red", fill = NA)+
    geom_sf(data = snapped3, col = "blue", alpha = 0.5)

Snpped lines in blue

对于任何更大的容忍度,结果与第三个相同. 显然,我对这个功能的工作方式有些不理解.我试图查看文档,但找不到任何有用的东西.我很感激对此有任何见解,以及如何通过该功能或其他功能获得我的预期结果.

谢谢.

推荐答案

st_snap snaps the vertices and segments of a geometry to another geometry's vertices.
[ ?sf::st_snap ]

如果您的多边形有长(类似)的直部分,则您很可能需要添加更多点(顶点),以便线要素的顶点和分段有一些可以捕捉的东西.而且,在咬合公差和所需的点密度之间找到最佳平衡可能有些棘手.

对于此reprex,st_segmentize(..., dfMaxLength = .2)st_snap(..., tolerance = .18)(即略小于dfMaxLength以避免错误快照)应该做:

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(ggplot2)
library(dplyr, warn.conflicts = FALSE)

poly <- st_polygon(list(cbind(c(0, 0, 1, 1, 0), c(0, 1, 1, 0, 0))))
st_as_text(poly)
#> [1] "POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))"

# add points to straight lines
poly_2 <- st_segmentize(poly, dfMaxLength = .2)
st_as_text(poly_2)
#> [1] "POLYGON ((0 0, 0 0.2, 0 0.4, 0 0.6, 0 0.8, 0 1, 0.2 1, 0.4 1, 0.6 1, 0.8 1, 1 1, 1 0.8, 1 0.6, 1 0.4, 1 0.2, 1 0, 0.8 0, 0.6 0, 0.4 0, 0.2 0, 0 0))"

lines <- st_multilinestring(list(
  cbind(c(0, 1), c(1.01, 1.05)),
  cbind(c(0, 1), c(-0.01, -.05)),
  cbind(c(1.025, 1.05, 1.025), c(1.05, .5, -.05))
))

snapped <- st_snap(lines, poly_2, tolerance = .18)
st_as_text(snapped)
#> [1] "MULTILINESTRING ((0 1, 0.2 1, 0.4 1, 0.6 1, 0.8 1, 1 1), (0 0, 0.2 0, 0.4 0, 0.6 0, 0.8 0, 1 0), (1 1, 1 0.8, 1 0.6, 1 0.4, 1 0.2, 1 0))"

ggplot() +
  geom_sf(data = st_cast(poly_2,"LINESTRING"), aes(color = "poly"), fill = NA, linewidth = 3) +
  geom_sf(data = st_cast(poly_2, "MULTIPOINT"), aes(color = "poly"), size = 6) +
  geom_sf(data = lines, aes(color = "lines"), linewidth = 2, alpha = .8) +
  geom_sf(data = st_cast(lines, "MULTIPOINT"), aes(color = "lines"), size = 4, alpha = .8) +
  geom_sf(data = snapped, aes(color = "snapped"), linewidth = 1, alpha = .8) +
  geom_sf(data = st_cast(snapped, "MULTIPOINT"), aes(color = "snapped"), size = 2, alpha = .8) +
  scale_color_viridis_d() +
  theme_minimal()


-e- Spatial join with nearest feature predicate

根据铁路网络数据集中要素的组织方式,与nearest feature个断言(st_join(... , join = st_nearest_feature))进行空间连接可能是解决原始问题的更好方法.然而,请留意有些模棱两可的情况,例如下面例子中的B D.

poly_sf <- 
  st_make_grid(poly, n = 2) |>
  st_sf(geometry = _) |>
  mutate(name = LETTERS[row_number()] |> as.factor())
poly_sf
#> Simple feature collection with 4 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 0 ymin: 0 xmax: 1 ymax: 1
#> CRS:           NA
#>                         geometry name
#> 1 POLYGON ((0 0, 0.5 0, 0.5 0...    A
#> 2 POLYGON ((0.5 0, 1 0, 1 0.5...    B
#> 3 POLYGON ((0 0.5, 0.5 0.5, 0...    C
#> 4 POLYGON ((0.5 0.5, 1 0.5, 1...    D

lines_sf <- 
  st_sfc(lines) |> 
  st_sf(geometry = _) |>
  st_cast("LINESTRING")
lines_sf
#> Simple feature collection with 3 features and 0 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 0 ymin: -0.05 xmax: 1.05 ymax: 1.05
#> CRS:           NA
#>                         geometry
#> 1    LINESTRING (0 1.01, 1 1.05)
#> 2  LINESTRING (0 -0.01, 1 -0.05)
#> 3 LINESTRING (1.025 1.05, 1.0...

lines_w_names <- 
  st_join(lines_sf, poly_sf, join = st_nearest_feature)
lines_w_names
#> Simple feature collection with 3 features and 1 field
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 0 ymin: -0.05 xmax: 1.05 ymax: 1.05
#> CRS:           NA
#>   name                       geometry
#> 1    C    LINESTRING (0 1.01, 1 1.05)
#> 2    A  LINESTRING (0 -0.01, 1 -0.05)
#> 3    D LINESTRING (1.025 1.05, 1.0...

ggplot() +
  geom_sf(data = poly_sf, aes(fill = name), show.legend = FALSE) +
  geom_sf_label(data = poly_sf, aes(label = name), alpha = .8) +
  geom_sf(data = lines_w_names, aes(color = name),  show.legend = FALSE, linewidth = 2) +
  geom_sf_label(data = lines_w_names, aes(label = name), alpha = .8) +
  scale_fill_viridis_d() +
  scale_color_viridis_d() +
  theme_minimal()

创建于2024年4月21日,共有reprex v2.1.0

R相关问答推荐

R:随机抽取所有可能排列的样本

R -列表元素中所有命名项的总和

在R中,将一个函数作为输入传递给另一个函数时进行参数判断

使用facet_wrap()时如何将面板标题转换为脚注?

在Julia中调用R函数

如何使用ggplot重新绘制LASO回归图?

高质量地将R格式的图表从Word中输出

如何将具有重复名称的收件箱合并到R中的另一列中,而结果不同?

使用ggplot 2根据R中的类别排列Likert比例gplot

ggplot geom_smooth()用于线性回归虚拟变量-没有回归线

如何在xyplot中 for each 面板打印R^2

迭代通过1个长度的字符串长字符R

如何在geom_col中反转条

在另存为PNG之前隐藏htmlwidget绘图元素

如何自定义3D散点图的图例顺序?

使用不同的定性属性定制主成分分析中点的 colored颜色 和形状

从数据创建数字的命名列表.R中的框

在散点图中使用geom_point放置线图例

R:如何在数据集中使用Apply

如果极点中存在部分匹配,则替换整个字符串