使用其他栅格作为指标对栅格堆栈的每个网格单元求和的功能 [英] Function to sum each grid cells of raster stack using other rasters as an indicator

问题描述

## 输入栅格s <- stack(list.files("~/dailyraster", full.names=TRUE)) # 每日栅格堆栈r_start <- raster("~/stackSumSTART.asc") # 这个栅格包含开始儒略日r_end <- raster("~/stackSumEND.asc") # 这个光栅包含结束儒略日noNAcells <- which(!is.na(r[])) # 包含值的单元格编号##虚拟光栅x<-rx[] <- 不适用## 环形对于(我在 noNAcells 中){x[i] <- sum(s[[r_start[i]:r_end[i]]][i])}

我想创建一个类似于 stackApply() 的函数，但我希望它在单元格基础上工作.

上面是一个for()循环版本，效果很好，但是太费时间了.

重点是每个单元格从两个栅格层获取sum()的范围，r_start，r_end在上面的脚本中.

现在我正在努力使用 apply() 系列来转换这段代码.

有没有可能用 for() 循环来提高速度?或者请给我一些在 apply()

中编写此代码的提示

任何意见都会帮助我，谢谢.

解决方案

你的方法

x <- s$layer.1系统时间(for (i in 1:ncell(x)) {x[i] <- sum(s[[r_start[i]:r_end[i]]][i], na.rm = T)})

<块引用>

 用户系统已用完0.708 0.000 0.710

我的提议

您可以在堆栈末尾添加用作索引的栅格，然后使用 calc 来大大加快进程速度(~30-50 倍).

s2 <- stack(s, r_start, r_end)sum_time <- function(x) {sum(x[x[6]:x[7]], na.rm = T)}系统时间(输出 <- calc(s2, fun = sum_time))

<块引用>

 用户系统已用完0.016 0.000 0.015

all.equal(x, output)

<块引用>

[1] 真

示例数据

库(光栅)# 生成随机值的栅格r1 <- r2 <- r3 <- r4 <- r5 <- r_start <- r_end <- raster(ncol=10, nrow=10)r1[] <- rnorm(ncell(r1), 1, 0.2)r2[] <- rnorm(ncell(r2), 1, 0.2)r3[] <- rnorm(ncell(r3), 1, 0.2)r4[] <- rnorm(ncell(r4), 1, 0.2)r5[] <- rnorm(ncell(r5), 1, 0.2)s <- 堆栈(r1，r2，r3，r4，r5)r_start[] <- 样本(1:2，ncell(r_start)，替换 = T)r_end[] <- 样本(3:5，ncell(r_end)，替换 = T)

## input raster
s <- stack(list.files("~/dailyraster", full.names=TRUE)) # daily raster stack
r_start <- raster("~/stackSumSTART.asc") # this raster contain starting Julian day
r_end <- raster("~/stackSumEND.asc") # this raster contain ending Julian day
noNAcells <- which(!is.na(r[])) # cell numbers which contain values

## dummy raster
x <- r
x[] <- NA 

## loop
for (i in noNAcells) {      
  x[i] <- sum(s[[r_start[i]:r_end[i]]][i])
}

I would like to create a function like stackApply(), but I want it to work on a cell basis.

Above is a for() loop version and it works well, but it takes too much time.

The point is that each cell gets the range of sum() from two raster layers, r_start, r_end in above script.

Now I am struggling to transform this code using apply() family.

Is there any possibility to improve the speed with for() loop? or please give me some tips to write this code in apply()

Any comments will help me, thank you.

解决方案

Your approach

x <- s$layer.1

system.time(
for (i in 1:ncell(x)) {
     x[i] <- sum(s[[r_start[i]:r_end[i]]][i], na.rm = T)
   }
)

   user  system elapsed 
  0.708   0.000   0.710

My proposal

You can add the rasters used as indices at the end of your stack and then use calc to highly speed up the process (~30-50x).

s2 <- stack(s, r_start, r_end)
sum_time <- function(x) {sum(x[x[6]:x[7]], na.rm = T)}

system.time(
   output <- calc(s2, fun = sum_time)
)

   user  system elapsed 
  0.016   0.000   0.015

all.equal(x, output)

[1] TRUE

Sample Data

library(raster)

# Generate rasters of random values
r1 <- r2 <- r3 <- r4 <- r5 <- r_start <- r_end <- raster(ncol=10, nrow=10)

r1[] <- rnorm(ncell(r1), 1, 0.2)
r2[] <- rnorm(ncell(r2), 1, 0.2)
r3[] <- rnorm(ncell(r3), 1, 0.2)
r4[] <- rnorm(ncell(r4), 1, 0.2)
r5[] <- rnorm(ncell(r5), 1, 0.2)
s <- stack(r1,r2,r3,r4,r5)

r_start[] <- sample(1:2, ncell(r_start),replace = T)
r_end[] <- sample(3:5, ncell(r_end),replace = T)

这篇关于使用其他栅格作为指标对栅格堆栈的每个网格单元求和的功能的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！