沿 SpatialLine 的 R 栅格提取：将提取的值与实际距离相关联答案

【问题标题】：R-raster extraction along SpatialLine: relate extracted values to actual distance沿 SpatialLine 的 R 栅格提取：将提取的值与实际距离相关联
【发布时间】：2015-05-08 07:25:07
【问题描述】：

在R中沿SpatialLine提取栅格值时，如何将这些值与沿这条线的实际距离相关联？

假设我想沿以下行提取 R 标志的值：

library(raster)
r <- raster(system.file("external/rlogo.grd", package="raster"))
x=c(5, 95)
y=c(20, 50)
line = SpatialLines(list(Lines(Line(cbind(x,y)), ID="a")))
plot(r)
plot(line, add=TRUE)

我可以提取值并绘制它们 - 但是如何用实际距离（例如从线左侧的 0 开始）替换 x 值（下面的1:length(vals)）？

vals <- extract(r, line)[[1]]
plot(1:length(vals), vals, type='o')

我可以将单元格的提取与xyFromCell 结合起来，以按照here 的建议获取提取的单元格的坐标，但我不清楚如何走得更远。

【问题讨论】：

标签： r extraction r-raster

【解决方案1】：

我不确定您到底要问什么，但是如果您要查找线段的最左侧坐标与线所通过的单元格中心之间的距离，那么您可以找到这样的距离：

x <- extract(r, l, cellnumbers=TRUE)[[1]]
xy <- xyFromCell(r, x[,1]) # get cell coordinates where the line passes
start <- xy[which.min(xy[,1]),] # leftmost coordinate of the line
d <- apply(xy, 1, function(x, start) sqrt(sum((x-start)^2)), start=start) # find distances between the line segment start and the cells
plot(1:length(d), d, type='o')

【讨论】：

谢谢@jvj。假设提取的值与线本身相关联，我更多的是考虑线本身的距离（我知道这不完全正确，并且还取决于用于提取的method，但不应该太大如果网格分辨率与线的长度相比较高，则会出现问题）。直观地说，我更多的是考虑线的最左边坐标与线上单元格中心的投影之间的距离（我想这很容易在 R 的其他地方存在） - 但也许以更高效/严谨/内置的方式...
顺便说一句@jvj，你上面的start 不是实际上不是（或者很可能不是，除非线路穿过单元中心）在线上？此外，start 不能保证定义为最左边的和最低单元格，是吗？还是我哪里错了？
要获取线路本身的距离，您也许可以使用 geosphere::dist2Line

【解决方案2】：

这是一个解决方案（部分基于@jvj 的输入），通过尝试计算raster::extract 提供的单元格中心在线上的正交投影，然后计算沿线的距离。

（这是一个 R 初学者脚本，可能很容易改进，但似乎可以工作（当然仅适用于投影距离相关的栅格））

vals <- extract(r, line, cellnumbers=TRUE)[[1]]
cellsxy <- xyFromCell(r, vals[,1]) # coordinates of intersected cells (likely not ON the line)
linexy = spsample(line, 1000, "regular") # get the line as points
linexy <- matrix(cbind(linexy$x, linexy$y), ncol=2) # easier than Spatial object for later
orthoproj <- c() # to store the orthogonal projections of cells centres on the line
for (i in 1:nrow(cellsxy)) {
  xypt = cellsxy[i,]
  min.index <- which.min(spDistsN1(linexy, xypt)) 
  orthopt <- linexy[min.index, ] # orthogonal projections = smaller distance 
  orthoproj <- c(orthoproj, c(orthopt[1], orthopt[2]))
}
orthoproj <- matrix(orthoproj, ncol=2, byrow=T)
orthoproj <- data.frame(x=orthoproj[,1], y=orthoproj[,2])
orthoproj <- orthoproj[order(orthoproj[,1]),] # reorder with increasing distance
orthoproj <- data.frame(x=orthoproj$x, y=orthoproj$y)
start <- linexy[which.min(linexy[,1]),] # leftmost coordinate of the line
dists <- apply(orthoproj, 1, 
               function(xy, start) sqrt(sum((xy-start)^2)), 
               start=start) # distances between 'start' and  the orthogonal projections

plot(dists, rev(vals[,2]), type='o') # !! beware: order of 'vals' and 'dists' 
                    # depending on the order in which cellnumbers are returned 
                    # in raster::extract and the shape of your line !!

【讨论】：