在 R 中，找到离另一点最近的线上点的坐标

Question

在 R 中，我有两个对象，一个 SpatialLinesDataFrame（表示道路网络）和一个 SpatialPointsDataFrame（表示对象位置）。我需要输出距离每个对象位置最近的道路上的坐标。

我所能找到的搜索结果都是其他语言（例如 Python - 如何找到线段上到任意点最近的点？）或者用于查找点与线之间最小距离的方法（例如使用 geosphere::dist2Line() 或 rgeos::gDistance()）。我只对返回线上最近点的坐标感兴趣。

编辑： 这是我的道路网络的一个小子集：

new("SpatialLinesDataFrame"
    , data = structure(list(ID = c(0, 0, 0, 0, 0, 0, 0),
                            ET_ID = c("4", "4", "4", "4", "4", "4", "4"),
                            length = c(0.280848, 0.812133, 0.0402004, 0.209611, 0.0433089, 0.501865, 0.363501)),
                       .Names = c("ID", "ET_ID", "length"),
                       row.names = c(980L, 982L, 983L, 984L, 987L, 988L, 989L),
                       class = "data.frame")
, lines = list(<S4 object of class structure("Lines", package = "sp")>, 
               <S4 object of class structure("Lines", package = "sp")>, 
               <S4 object of class structure("Lines", package = "sp")>, 
               <S4 object of class structure("Lines", package = "sp")>, 
               <S4 object of class structure("Lines", package = "sp")>, 
               <S4 object of class structure("Lines", package = "sp")>, 
               <S4 object of class structure("Lines", package = "sp")>)
, bbox = structure(c(433266.568837884, 547825.73420664, 437050.511867258, 548168.921069476),
                   .Dim = c(2L, 2L),
                   .Dimnames = list(c("x", "y"), c("min", "max")))
, proj4string = new("CRS", projargs = "+proj=utm +zone=37 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
)

我的物体位置：

new("SpatialPointsDataFrame"
    , data = structure(list(x = c(38.4129, 38.41697, 38.41501), y = c(4.95659, 4.95809, 4.96122)),
                       .Names = c("x", "y"), row.names = c(105L, 166L, 185L), class = "data.frame")
    , coords.nrs = numeric(0)
    , coords = structure(c(434912.0166297, 435363.392542353, 435146.398500838, 547894.637850701, 548060.055746692, 548406.25007762),
                         .Dim = c(3L, 2L), .Dimnames = list(c("105", "166", "185"), c("x", "y")))
    , bbox = structure(c(434912.0166297, 547894.637850701, 435363.392542353, 548406.25007762),
                       .Dim = c(2L, 2L), .Dimnames = list(c("x", "y"), c("min", "max")))
    , proj4string = new("CRS", projargs = "+proj=utm +zone=37 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
)

提前感谢！

- Andrew

---

也许你愿意分享一些你的数据的例子吗？ - mtoto

好的，我会尝试将其输出。 - Andrew

@Spacedman 是的，我认为这是我想要的，但对于地理坐标并不自动工作。我猜我可以使用 @coords 来表示点的坐标，但我不确定如何输出线的坐标。 - Andrew

啊，snapPointsToLines适用于Spatial*对象。由于它不适用于纬度和经度坐标，所以需要使用几何（投影）坐标。只有在处理全球范围很大的数据集时才会出现问题... - Spacedman

谢谢- maptools::snapPointsToLines 可以完成任务。我的研究区域很小，所以没有范围的问题。 - Andrew

Answer 1

如果您不介意自己实现，可以从这个伪代码进行移植...

Point snap_to_segment(Point point, Segment segment):
    Point s1 = segment.start;
    Point s2 = segment.end;

    Vector v = s2 - s1;
    Vector w = point - s1;

    double c1 = Vector.dot_product(w,v);
    double c2 = Vector.dot_product(v,v);

    Point snap;

    if c1 <= 0:
        snap = s1

    elif c2 <= c1:
        snap = s2

    else:
        snap = s1 + v * (c1 / c2)

    return snap

除此之外，当我决定使用Shapely时，我在Python中比上述函数做得好多了。事实证明： rgeos是R语言对应于Python的Shapely库。Shapely和rgeos都基于GEOS（即PostGIS引擎）（source）

然而，在Shapely中，我这样找到所需的点：

desiredPoint = road.distance(road.project(point));

但是，rgeos似乎缺少linestring.project(point)和linestring.distance(scalarValue)方法...

- heltonbiker

Answer 2

这是heltonbiker在R中的答案的一个快速hack（2D版本），供那些可能需要的人使用。当然，最好使用一个库，但有时你不想使用它，如果你只需要得到一条线上最近的点，这个方法就可以做到。

findClosestPtToSeg <- function( ptx,pty, sx0,sy0,sx1,sy1 ){

  sdx <- sx1-sx0
  sdy <- sy1-sy0
  wx <- ptx-sx0
  wy <- pty-sy0

  lambda <- (wx*sdx + wy*sdy) / (sdx*sdx + sdy*sdy)
  trunclamb <- min(max(lambda,0),1)

  rv <- list() # return the coordinates in a list
  rv$x <- sx0 + trunclamb*sdx
  rv$y <- sy0 + trunclamb*sdy
  return(rv)
}

以下是一些在几个测试案例中的输出结果：

> findClosestPtToSeg( 1,2, 0,0,1,1 )
$x
[1] 1
$y
[1] 1

> findClosestPtToSeg( -1,-2, 0,0,1,1 )
$x
[1] 0
$y
[1] 0

> findClosestPtToSeg( 0.5,0.6, 0,0,1,1 )
$x
[1] 0.55
$y
[1] 0.55

请注意这不是优美的代码，它没有检查退化区间，并且处理输入点和返回值的方式不同（坐标传递 vs 列表打包）。哦好吧，这只是一个快速的hack，对我来说效果还不错。

如果这对heltonbiker有所帮助，请确保他也获得一些赞吧。

- Mike Wise

Answer 3

这是 heltonbiker 方法的向量化 n 维版本。它还返回从 p4 到 p3 的距离以及从 p4 到 p1 的距离相对于从 p1 到 p2 的距离的分数，其中 p3 是目标点，p4 是找到的最近点。

nrow=1e5
ncol=10

p1=matrix(rnorm(ncol*nrow),nrow)
p2=matrix(rnorm(ncol*nrow),nrow)
p3=rnorm(ncol)

p21=p2-p1
p31=outer(rep(1,nrow(p1)),p3)-p1
c1=rowSums(p31*p21)
c2=rowSums(p21*p21)
frac=c1/c2
p4=p1+p21*frac

# uncomment to restrict points to the line segment between p1 and p2
# o1=c1>c2;u0=c1<0;p4[o1]=p1[o1];frac[o1]=1;p4[u0]=p2[u0];frac[u0]=0

dist=sqrt(outer(rowSums(p4^2),sum(p3^2),"+")-tcrossprod(p4,2*t(p3)))[,1]

outer(rep(1,nrow(m)),v)-m 是一种快速的方法，用于在向量和矩阵的每一行之间应用逐行算术运算（参见/questions/39443055/add-a-vector-to-all-rows-of-a-matrix）。

sqrt(outer(rowSums(m^2),sum(v^2),"+")-tcrossprod(m,2*t(v)))[,1] 是一种快速的方法，用于计算向量v到矩阵m中每一行的欧几里得距离。

非向量化：

p1=c(1,5,2,4)
p2=c(7,3,2,5)
p3=c(5,5,2,3)

p21=p2-p1
p31=p3-p1
c1=(p31%*%p21)[1]
c2=(p21%*%p21)[1]
frac=c1/c2
p4=p1+p21*frac
# if(c1>c2){p4=p1;frac=0};if(c1<0){p4=p1;frac=1}
dist=sqrt(sum((p31-frac*p21)^2))
list(p4,frac,dist)

- nisetama