R Raster: Расчет расстояния от каждой ячейки сетки до других растровых ячеек в радиусе

Question

У меня есть растровый объект, который содержит числовые значения и NA, и я хочу знать расстояние от каждой растровой ячейки со значениями до других числовых ячеек в радиусе.Результатом является матрица / растр, который содержит суммированное расстояние для каждой числовой растровой ячейки, которые не являются NA.Клетки, которые не считаются, получают NA.Проблема в том, что вычисление занимает относительно много времени, потому что мне нужно сделать цикл по всем растровым ячейкам, чтобы проверить, являются ли они числовыми, а затем вычислить расстояния.После этого мне нужно выбрать все растровые ячейки, которые являются числовыми, и выбрать все числовые ячейки в пределах определенного радиуса, так как функция accCost не учитывает ячейки с NA для вычисления расстояний.Существует ли более быстрый способ вычисления суммированных расстояний растровых ячеек в пределах определенного радиуса?

Во-первых, у меня есть растр, который я должен изменить, потому что я хочу знать только суммарное расстояние ячеек, которые лежат в определенной области.Поскольку функция accCost не учитывает NA, мне нужно дать им значение.Затем я определяю ядра для функции foreach.Поскольку я использую функцию «accCost» для вычисления расстояний от одной растровой ячейки до других, мне нужно сделать некоторые настройки по умолчанию и рассчитать координаты xy ячеек.Чтобы выбрать расстояния, которые не являются NA, я делаю булевский запрос.Затем я проверяю каждую растровую ячейку, чтобы проверить, имеют ли они определенное значение.Если да, то я вычисляю расстояние до каждой ячейки сетки с помощью функции «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 

 }
}

Результат выглядит так: Результат

Answer 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)