R / GIS：查找位置与最近线之间的正交距离答案

【问题标题】：R/GIS: Find orthogonal distance between a location and nearest lineR / GIS：查找位置与最近线之间的正交距离
【发布时间】：2012-04-04 14:16:47
【问题描述】：

我试图找出一组位置坐标和一组线（道路或河流）之间的正交距离。点集采用纬度/经度对的形式，线条位于 shapefile (.shp) 中。使用maptools 或PBSmapping 在地图上绘制它们不是问题。但我的基本问题是找到一个人从某个位置到达道路或河流的最小距离。有没有办法在 R 中做到这一点？

【问题讨论】：

标签： r gis

【解决方案1】：

如果我理解正确，您可以使用 gDistance 包中的 gDistance 来完成此操作。

将行读为SpatialLines/DataFrame，将点读为SpatialPoints/DataFrame，然后遍历每个点，每次计算距离：

require(rgeos)
## untested code
shortest.dists <- numeric(nrow(sp.pts))
for (i in seq_len(nrow(sp.pts)) {
    shortest.dists[i] <- gDistance(sp.pts[i,], sp.lns)
}

这里sp.pts 是空间点对象，sp.lns 是空间线对象。

您必须循环，以便仅将sp.pts 中的单个坐标与sp.lns 中的所有线条几何图形进行比较，否则您将获得所有点的聚合值的距离。

由于您的数据是纬度/经度，您应该将线和点都转换为合适的投影，因为gDistance 函数假定笛卡尔距离。

更多讨论和示例（编辑）

获得直线上最近的点而不是距离会很整洁，但这会打开另一个选项，即您是否需要沿直线最近的坐标，或实际交点与比任何现有顶点更近的线段。如果您的顶点足够密集以至于差异无关紧要，则使用 spDistsN1 包中的 sp 。您必须从集合中的每一行中提取所有坐标（不难，但有点难看），然后遍历每个兴趣点，计算到线顶点的距离 - 然后您可以找到最短的并选择该坐标来自一组顶点，因此您可以轻松获得距离和坐标。也不需要投影，因为该函数可以使用带有longlat = TRUE 参数的椭球距离。

library(maptools)

## simple global data set, which we coerce to Lines
data(wrld_simpl)

wrld_lines <- as(wrld_simpl, "SpatialLinesDataFrame")

## get every coordinate as a simple matrix (scary but quick)
wrld_coords <- do.call("rbind", lapply(wrld_lines@lines, function(x1) do.call("rbind", lapply(x1@Lines, function(x2) x2@coords[-nrow(x2@coords), ]))))

以交互方式检查它，您必须对其进行修改以保存坐标或最小距离。这将绘制线条并等待您单击绘图中的任意位置，然后它会从您的单击到线条上最近的 vertex 绘制一条线。

## no out of bounds clicking . . .
par(mar = c(0, 0, 0, 0), xaxs = "i", yaxs = "i") 

plot(wrld_lines, asp = "")

n <- 5

for (i in seq_len(n)) {
xy <- matrix(unlist(locator(1)), ncol = 2)
    all.dists <- spDistsN1(wrld_coords, xy, longlat = TRUE)
    min.index <- which.min(all.dists)
    points(xy, pch = "X")
lines(rbind(xy, wrld_coords[min.index, , drop = FALSE]), col = "green", lwd = 2)
}

【讨论】：

不错且内容丰富的答案。感谢那。看起来rgeos 没有任何返回（或）最近位置的函数。你碰巧有一个功能的建议吗？看起来 project2segment() in spatstat 这样做，但不幸的是它不使用 sp 类对象，这会更方便......
我会更新的，我正在考虑更多，有一些微妙之处。
spatstat和sp之间的转换并不难，maptools中有强制转换工具，比如spst
“更多讨论”部分中的精彩演示。我希望对此感兴趣的每个人都花一点时间来运行该代码。另外，感谢您提供指向 maptools 包中的全套 sp -> spatstat -> sp 转换工具的指针。
mdsumner... 真的，非常巧妙。非常感谢。这一定是 stackoverflow 上最原始的想法之一。

【解决方案2】：

geosphere 包具有 dist2line 函数，可以对 lon/lat 数据执行此操作。它可以使用 Spatial* 对象或矩阵。

line <- rbind(c(-180,-20), c(-150,-10), c(-140,55), c(10, 0), c(-140,-60))
pnts <- rbind(c(-170,0), c(-75,0), c(-70,-10), c(-80,20), c(-100,-50), 
         c(-100,-60), c(-100,-40), c(-100,-20), c(-100,-10), c(-100,0))
d <- dist2Line(pnts, line)
d

结果说明

plot( makeLine(line), type='l')
points(line)
points(pnts, col='blue', pch=20)
points(d[,2], d[,3], col='red', pch='x')
for (i in 1:nrow(d)) lines(gcIntermediate(pnts[i,], d[i,2:3], 10), lwd=2)

【讨论】：

【解决方案3】：

看起来这可以在 sf 包中使用 st_distance 函数完成。

您将两个 sf 对象传递给函数。与其他解决方案相同的问题是您需要迭代您的点，以便该函数计算每个点到道路上每个点之间的距离。然后取最短距离的结果向量的最小值。

# Solution for one point
min(st_distance(roads_sf, points_sf[1, ]))

# Iterate over all points using sapply
sapply(1:nrow(points_sf), function(x) min(st_distance(roads_sf, points_sf[x, ])))

【讨论】：