foreach 不会改变 R 中栅格单元的值
foreach doesn't change value of raster cell in R
我正在尝试在 R 中模拟羊群行为。
这是代码
library(raster)
library(sp)
library(foreach)
K=100
sig=0.2
G=0.3
x <- raster(ncol=2000,nrow=2000)
values(x) <- sign(rnorm(4000000,mean=0,sd=0.3))
y <- raster(ncol=2000,nrow=2000)
values(y) <- sign(rnorm(4000000,mean=0,sd=0.3))
#plot(x)
ei <- rnorm(4000000)
j=0
while(j < 30) {
for(i in 1:4000000){
ad <- adjacent(x,cell=c(i))[,2]
y[i] <- sign(K*sum(x[ad])+sig*ei[i]+G)
}
x <- y
plot(x)
j = j+1
}
经典的循环方法太慢了。
如果我使用 foreach 循环而不是经典的 for 循环,它不会在每次迭代中更改 y 的值。
我根本无法修复它。
有人可以帮忙吗?
谢谢
您有一个动态模型,其中每个(时间)步骤的输出都是下一步的输入。并行执行此操作是不可能的。但这并不意味着您不能使模型 运行 更快。
在 R 中遍历栅格单元格总是很慢,所以我们需要避免这种情况。通常这样的问题可以用 focal
解决(见底部代码)---但在这种情况下很难,因为你有效地使用了两个栅格(x
和 ei
) --- 我将着眼于在 terra
包中实现多层焦点操作。
这是 getFocalValues
的一种方法。它 快 很多(我使用 Sys.sleep
来减慢它的速度)。
library(raster)
set.seed(0)
x <- raster(ncol=200, nrow=200)
values(x) <- sign(rnorm(ncell(x),mean=0,sd=0.3))
y <- raster(x)
values(y) <- sign(rnorm(ncell(x),mean=0,sd=0.3))
ei <- rnorm(ncell(x))
K=100
sig=0.2
G=0.3
for (j in 1:29) {
# with large rasters, you may need to do the below in chunks
v <- getValuesFocal(x, 1, nrow(x), c(3,3))
# only keep the rook neighbors
v <- v[, c(2,4,6,8)]
v <- rowSums(v, na.rm=TRUE)
values(x) <- sign(K*v+sig*ei+G)
plot(x)
Sys.sleep(0.1)
}
这就是您在类似情况下如何使用 focal
w <- matrix(c(0,1,0,1,0,1,0,1,0), 3, 3)
y <- focal(x, w, fun=function(i)sign(K*sum(i)+sig+G))
另请参阅 ?focal
中的元胞自动机示例
我正在尝试在 R 中模拟羊群行为。
这是代码
library(raster)
library(sp)
library(foreach)
K=100
sig=0.2
G=0.3
x <- raster(ncol=2000,nrow=2000)
values(x) <- sign(rnorm(4000000,mean=0,sd=0.3))
y <- raster(ncol=2000,nrow=2000)
values(y) <- sign(rnorm(4000000,mean=0,sd=0.3))
#plot(x)
ei <- rnorm(4000000)
j=0
while(j < 30) {
for(i in 1:4000000){
ad <- adjacent(x,cell=c(i))[,2]
y[i] <- sign(K*sum(x[ad])+sig*ei[i]+G)
}
x <- y
plot(x)
j = j+1
}
经典的循环方法太慢了。 如果我使用 foreach 循环而不是经典的 for 循环,它不会在每次迭代中更改 y 的值。
我根本无法修复它。 有人可以帮忙吗?
谢谢
您有一个动态模型,其中每个(时间)步骤的输出都是下一步的输入。并行执行此操作是不可能的。但这并不意味着您不能使模型 运行 更快。
在 R 中遍历栅格单元格总是很慢,所以我们需要避免这种情况。通常这样的问题可以用 focal
解决(见底部代码)---但在这种情况下很难,因为你有效地使用了两个栅格(x
和 ei
) --- 我将着眼于在 terra
包中实现多层焦点操作。
这是 getFocalValues
的一种方法。它 快 很多(我使用 Sys.sleep
来减慢它的速度)。
library(raster)
set.seed(0)
x <- raster(ncol=200, nrow=200)
values(x) <- sign(rnorm(ncell(x),mean=0,sd=0.3))
y <- raster(x)
values(y) <- sign(rnorm(ncell(x),mean=0,sd=0.3))
ei <- rnorm(ncell(x))
K=100
sig=0.2
G=0.3
for (j in 1:29) {
# with large rasters, you may need to do the below in chunks
v <- getValuesFocal(x, 1, nrow(x), c(3,3))
# only keep the rook neighbors
v <- v[, c(2,4,6,8)]
v <- rowSums(v, na.rm=TRUE)
values(x) <- sign(K*v+sig*ei+G)
plot(x)
Sys.sleep(0.1)
}
这就是您在类似情况下如何使用 focal
w <- matrix(c(0,1,0,1,0,1,0,1,0), 3, 3)
y <- focal(x, w, fun=function(i)sign(K*sum(i)+sig+G))
另请参阅 ?focal
中的元胞自动机示例