R Raster：计算从每个网格单元到半径内其他栅格单元的距离答案

【问题标题】：R Raster: Calculation Distance from each grid cell to other raster cells within a radiusR Raster：计算从每个网格单元到半径内其他栅格单元的距离
【发布时间】：2019-10-02 10:46:35
【问题描述】：

我有一个包含数值和 NA 的栅格对象，我想知道从每个带有值的栅格像元到半径内其他数值像元的距离。结果是一个矩阵/栅格，其中包含每个非 NA 数值栅格像元的总距离。未考虑的单元格得到 NA。现在的问题是计算需要相对较长的时间，因为我必须对所有栅格单元进行循环以测试它们是否为数字，然后计算距离。之后，我必须选择所有数字栅格单元格，并选择一定半径内的所有数字单元格，因为 accCost 函数不考虑具有 NA 的单元格来计算距离。有没有更快的方法来计算一定半径内栅格单元的总距离？

首先，我有一个栅格，我必须对其进行修改，因为我只想知道位于特定区域内的单元格的总距离。由于 accCost 函数不考虑 NA，我需要给它们一个值。然后我定义了“foreach”函数的核心。当我使用“accCost”函数计算一个栅格像元到其他像元的距离时，我需要进行一些默认设置并计算像元的 xy 坐标。为了选择不是 NA 的距离，我进行了 boolena 查询。然后我循环遍历每个栅格单元以测试它们是否具有特定值。如果是，那么我使用“accCost”函数计算到每个网格单元的距离。然后我对生成的栅格进行子集化。否则，栅格单元将获得 NA。

#load library
library(doParallel)
library(foreach)
library(gdistance)
library(raster)

#Create a raster
mraster<-matrix(nrow=15,ncol = 15)
mraster[c(5,10,9,15,13,5),c(4,8,9,7,7,15)]<-1
builtupraster<-raster(mraster, xmn=1, xmx=1500, ymn=1, ymx=1500)
proj4string(builtupraster) <- CRS("+proj=somerc +lat_0=46.952405555556 +lon_0=7.439583333333 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs")
builtupraster[is.na(builtupraster)]<-0

#define parallel function
cores<-detectCores()
cl<-makeCluster(cores-2)
doParallel::registerDoParallel(cl)

#calculate the distance between cells fast
r <- builtupraster
r2 <- transition(r, transitionFunction = function(x){1}, directions = 16)
r2 <- geoCorrection(r2,  scl=FALSE)
siedlungsfl<-as.matrix(r>0)

xcord<-xFromCol(r,1:ncol(r))
ycord<-yFromRow(r,1:nrow(r))

testrow<-1:nrow(r)
testcol<-1:ncol(r)

distmat<-foreach(row= testrow, .combine = "rbind" ,.packages = "gdistance") %:%
 foreach(col= testcol, .combine="c",.packages = "gdistance") %dopar%{
 if(r[row,col]>0){
   d <- accCost(r2,c(xcord[row],ycord[col]))
   d2 <- d[which(siedlungsfl)]
   d3 <- d2[d2<=2000]
   d4 <- sum(sqrt(2*d3+1)-1)/length(d3)+(sqrt(0.97428*30+1.46)-0.996249)

 } else{

   d4 <- NA 

 }
}

结果是这样的：Result

【问题讨论】：

您应该提供一个小而尽可能简单的可重现示例，其中包含一些输入数据（使用代码生成）和正确的输出数据。
@RobertHijmans 我现在已经用一个例子调整了代码。我将输出结果添加为图像文件。谢谢你的建议！

标签： r loops foreach distance raster

【解决方案1】：

我已经简化了您的脚本 --- 删除了并行化，以便更容易优化。我删除了一个循环和一个 if 语句，并进行了一些其他更改；这可能会有所帮助，但可能不是很大。

library(gdistance)

r <- raster(nrow=15,ncol = 15, xmn=1, xmx=1500, ymn=1, ymx=1500, vals=0, crs = "+proj=somerc +lat_0=46.952405555556 +lon_0=7.439583333333 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m"
)
r[c(5,10,9,15,13,5),c(4,8,9,7,7,15)] <- 1

tr <- transition(r, transitionFunction = function(x){1}, directions = 16)
tr <- geoCorrection(tr,  scl=FALSE)

cells <- Which(r > 0, cells=TRUE)
xy <- xyFromCell(r, cells)

constant <- sqrt(0.97428 * 30 + 1.46) - 0.996249
d <- rep(NA, length(cells))
for(i in 1:length(cells)) {
    dst <- accCost(tr, xy[i,])[cells]
    dst <- dst[dst <= 1000] # 2000 is too high for this example
    d[i] <- mean(sqrt(2 * dst + 1) - 1)
}

result <- raster(r)
result[cells] <- d + constant
plot(result)
text(result)

另一种可能较慢的方法可能是这样的：

library(raster)
r <- raster(nrow=15,ncol = 15, xmn=1, xmx=1500, ymn=1, ymx=1500, crs = "+proj=somerc +lat_0=46.952405555556 +lon_0=7.439583333333 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m"
)
r[c(5,10,9,15,13,5),c(4,8,9,7,7,15)] <- 1
d <- distance(r)

cells <- Which(r > 0, cells=TRUE)
b <- buffer(SpatialPoints(xyFromCell(r, cells)), 1000, dissolve=FALSE)
e <- extract(d, b, fun=sum)

result2 <- raster(r)
result2[cells] <- e
plot(result2)
text(result2)

结果不一样，因为我（根据你的描述）对距离求和，而你做constant + mean(sqrt(2 * dst + 1) - 1)

【讨论】：