如何根据 R 中未定义和未知的点空间序列创建 SpatialLines 对象？

Question

我有 lon/lat 点，问题是我没有行 ID 或行内的排序 ID。所以我想根据点之间的距离标准从点制作多线。例如在图像 here 中，至少应创建 4 行。

我怎样才能做到这一点？

Answer 1

我认为这个问题没有简单的通用解决方案。这是使用 spatstat 实现目标的一些灵感：

library(spatstat)

生成测试数据：

d1 <- data.frame(x=1:9, y=-4)
d2 <- data.frame(x=1:9, y=0)
d3 <- data.frame(x=1:9, y=4)
X <- as.ppp(rbind(d1,d2,d3), W = owin(c(0,10), c(-5,5)))
set.seed(42)
X <- rjitter(X, 0.5)
plot(X)

找到连接的组件并在每个组件内将每个点连接到最近的两个邻居：

Xcomp <- connected.ppp(X, R = 2)
Xcomp <- split(Xcomp)
neighbours <- list()
line_list <- list()
for(i in seq_along(Xcomp)){
  pts <- Xcomp[[i]]
  nn <- nnwhich(pts, k=1:2)
  x0 <- c(pts$x, pts$x)
  y0 <- c(pts$y, pts$y)
  x1 <- c(pts$x[nn[,1]], pts$x[nn[,2]])
  y1 <- c(pts$y[nn[,1]], pts$y[nn[,2]])
  line_list[[i]] <- psp(x0, y0, x1, y1, window = Window(X))
}

将组件放回一起并转换为线性网络 (linnet) 基本上是一个无向图，其中节点在欧几里德 space 而不是抽象。

L <- Reduce(superimpose.psp, line_list)
L <- as.linnet(L)
#> Warning: Network is not connected
plot(L)

剩下的任务是找到所有三角形并删除最长的边哪个更繁琐。您可以使用 edges2triangles 来查找所有三角形：

tri <- edges2triangles(L$from, L$to)
tri
#>       i  j  k
#> [1,]  1  2  3
#> [2,]  4  5  6
#> [3,]  7  8  9
#> [4,] 10 11 12
#> [5,] 13 14 15
#> [6,] 16 17 18
#> [7,] 19 20 21
#> [8,] 25 26 27

例如顶点 25,26,27 形成一个三角形

i <- as.numeric(tri[8,])
Li <- thinNetwork(L, retainvertices = i)
plot(Li)

三角形从 i 到 j 有 3 条边：

j <- i[c(2,3,1)]
i
#> [1] 25 26 27
j
#> [1] 26 27 25

所有顶点之间的距离矩阵（矫枉过正，但很容易计算和应该只做一次——避免大数据集）

D <- pairdist(vertices(L))

最长距离索引：

long <- which.max(diag(D[i,j]))
long
#> [1] 1

因此应该删除从i[long]到j[long]的边

plot(L)
edge <- which(paste(L$from,L$to)==paste(sort(c(i[long],j[long])), collapse = " "))
plot(thinNetwork(L, retainedges = edge), add = TRUE, col = 2, lwd = 1.5)

我们应该将此代码应用于循环中的所有三角形：

edge <- numeric(nrow(tri))
for(k in seq_len(nrow(tri))){
  i <- tri[k,]
  j <- i[c(2,3,1)]
  long <- which.max(diag(D[i,j]))
  edge[k] <- which(paste(L$from,L$to)==paste(sort(c(i[long],j[long])), collapse = " "))
}
Lfinal <- thinNetwork(L, retainedges = -edge)
plot(Lfinal)

如果单独需要这些行，我们可以使用 connected:

Lfinal_list <- connected.linnet(Lfinal, what = "components")
Lfinal_list
#> [[1]]
#> Linear network with 9 vertices and 8 lines
#> Enclosing window: rectangle = [0, 10] x [-5, 5] units
#> 
#> [[2]]
#> Linear network with 9 vertices and 8 lines
#> Enclosing window: rectangle = [0, 10] x [-5, 5] units
#> 
#> [[3]]
#> Linear network with 9 vertices and 8 lines
#> Enclosing window: rectangle = [0, 10] x [-5, 5] units

查找和删除三角形可以很容易地为每个组件完成在收集所有行时构建行而不是在最后。这在大型数据集上会更有效，但这很好用作为概念证明。 谨防黑客攻击，例如上面的粘贴技巧找到边缘数——这可能不是很稳健，我不确定是否它适用于所有情况。

如何根据 R 中未定义和未知的点空间序列创建 SpatialLines 对象？

How to make a SpatialLines object from undefined and unknowon spatial sequence of points in R?

r

spatial

sp