r栅格砖中的单元格和值由两个不同的栅格决定，如何加快计算速度？

问题描述 投票：0回答：1

我的工作是处理有日常数据的气候数据文件，所以对于大多数年份来说，365个栅格在一块砖头上，我想对文件中的天数子集进行求和，比如说x天到y天。我想对文件中的天数子集进行求和--比如说从第x天到第y天，这可以通过stackApply来实现。我已经创建了下面的一些代码，生成一些栅格，创建一个砖块，并使用x和y、1和3的特定值应用stackApply。

然而，我需要的是x和y取自两个栅格层。在下面的代码中，它们被称为raster.start和raster.end。在第一组代码下面，我有第二组代码，它可以工作，但速度很慢。

library(raster)
r <- raster(nrows=100, ncols=100)
s <- stack(lapply(1:5, function(i) setValues(r, runif(ncell(r), min = -10*i, max = 10))))
raster.start <- setValues(r, sample(2, ncell(r), replace=TRUE))
raster.end <- raster.start + 3
rasterb <- brick(s)

indices <- format(as.Date(names(rasterb), format = "layer.%d"), format = "%d")
indices <- c(1,1,1,1,1)

datasum.all <- stackApply(rasterb, indices, fun = sum)
datasum.sub1 <- stackApply(rasterb[[c(1:3)]], indices, fun = sum)

我们的想法是通过开始和结束栅格的行和列来步入砖块的子集，并对其进行操作。这是我开发的代码，可以做到这一点。

raster.out <- r
for (i in 1:nrow(r)){
  for (j in 1:ncol(r)){
    start <- raster.start[[1]][i,j] # get the starting day
    end <- raster.end[[1]][i,j] # get the ending day
    raster.out[i,j] <- sum(rasterb[[start:end]][i,j])
  }
}

然而，即使是这个玩具例子，计算时间也很慢。花了大约1.3分钟才完成。我试着用函数替换了一些代码，如下所示，但对完成时间没有影响。任何关于如何加快这个过程的建议都非常感激。

startEnd <- function(raster.start, raster.end, i,j) {
  start <- raster.start[i,j] # get the starting day
  end <- raster.end[i,j] # get the ending day
  return(c(start,end))
}

rasterOutValue <- function(rasterb, i, j, startEnd){
  return(sum(rasterb[[startEnd]][i,j]))
}

for (i in 1:nrow(raster.in1)){
  for (j in 1:ncol(raster.in1)){
    raster.out[i,j] <-rasterOutValue(rasterb, i, j, startEnd(raster.start, raster.end, i,j))
  }
}
r r-raster
1个回答
3
投票

您的示例数据

library(raster)
r <- raster(nrows=100, ncols=100)
set.seed(88)
b <- stack(lapply(1:5, function(i) setValues(r, runif(ncell(r), min = -10*i, max = 10))))
r.start <- setValues(r, sample(2, ncell(r), replace=TRUE))
r.end <- raster.start + 3

首先是你的例子的改进版，可以用，但太慢了。下面的速度大大加快，但还是相当慢。

raster.out <- r
for (i in 1:ncell(r)){
    start <- raster.start[i] # get the starting day
    end <- raster.end[i] # get the ending day
    raster.out[i] <- sum(rasterb[i][start:end])
}

这让我的时间从74秒降到了5秒。但是你永远不应该在单元格上循环，那总是会太慢。相反，你可以这样做（对我来说是0.04秒）。

s <- stack(r.start, r.end, b)
x <- calc(s, fun=function(x) sum(x[(x[1]:x[2])+2]))
#class      : RasterLayer 
#dimensions : 100, 100, 10000  (nrow, ncol, ncell)
#resolution : 3.6, 1.8  (x, y)
#extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#source     : memory
#names      : layer 
#values     : -129.5758, 30.31813  (min, max)

这似乎是正确的

a <- s[1]
a
#     layer.1.1 layer.2.1 layer.1.2 layer.2.2  layer.3   layer.4   layer.5
#[1,]         1         4 -1.789974  2.640807 4.431439 -23.09203 -5.688119    

fun <- function(x) sum(x[(x[1]:x[2])+2])
fun(a)
#[1] -17.80976
x[1]
#[1] -17.80976

calc 对光栅对象来说，就像 apply 是对矩阵。(这就是为什么它被称为 app 在 terra.

开始的地方是先写一个函数来实现你对一个向量的要求。

x <- 1:10
test1 <- function(start, end, values) {
    mean(values[start:end]) 
}
test1(2, 5, x)
test1(5, 8, x)

calc 只接受一个参数，所以像这样的函数是 

test2 <- function(values) {
    # the +2 to skip the first two elements in the computation
    start <- values[1] + 2
    end <- values[2] + 2
    mean(values[start:end]) 
}

test2(c(2, 5, x))
test2(c(5, 8, x))

还有一个更简洁的版本 

test3 <- function(v) {
    mean(v[ (v[1]:v[2])+2 ] ) 
}
 test3(c(2, 5, x))
 #[1] 3.5
 test3(c(5, 8, x))
 #[1] 6.5

第二项补充（提醒大家一定要用NA值检查！）。test3 当其中一个指数（开始和结束）为 NA (如果是其他的也可以 NA)

test3(c(NA, 5, x))
#Error in v[1]:v[2] : NA/NaN argument

所以我们需要一个函数来捕获这些

test4 <- function(v) {
    if (any(is.na(v[1:2]))) {
        NA
    } else {
        mean(v[ (v[1]:v[2])+2 ] ) 
    }
}

test4(c(NA, 5, x))
#[1] NA
test4(c(1, 5, x))
#[1] 3

通常 "开始 "和 "结束 "都会是 NA 同时，所以一个更简单的版本，应该也是可以的。

test5 <- function(v) {
    if (is.na(v[1])) {
        NA
    } else {
        mean(v[ (v[1]:v[2])+2 ] ) 
    }
}

这种方法与 calc 可能会很慢，因为它把一个RasterBrick变成了一个有365+2层的RasterStack。这将大大降低数据的读取速度。所以你可以尝试用这种方法 overlay 取而代之（这里用 sum 再次)

f <- function(i, v) {
    j <- !is.na(i[,1])
    r <- rep(NA, nrow(i))
    x <- cbind(i[j,,drop=FALSE], v[j,,drop=FALSE])
    r[j] <- apply(x, 1, function(y) sum(y[ (y[1]:y[2])+2 ] )) 
    r
}
cal <-stack(r.start, r.end)
x <- overlay(cal, b, fun= f, recycle=FALSE)
x
#class      : RasterLayer 
# ...
#values     : -129.5758, 30.31813  (min, max)

你可以通过用RcppC++编写来加快算法的速度。

library(Rcpp)
cppFunction('std::vector<double> gtemp(NumericMatrix cal, NumericMatrix wth) {
    std::vector<double> out(cal.nrow(), NAN);
    for (int i=0; i<cal.nrow(); i++) {
      if (!std::isnan(cal(i,0))){
         NumericVector v = wth(i,_);
         size_t start = cal(i,0)-1;
         size_t end = cal(i,1);
         out[i] = std::accumulate(v.begin()+start, v.begin()+end, 0.0);
      }  
    }
    return out;
}')

x <- overlay(cal, b, fun=gtemp, recycle=FALSE)

以下是您如何使用 terra (版本>= 0.6-14)和。rapp (range-apply)方法。

示例数据 

library(terra)
d <- rast(nrows=100, ncols=100, nl=5)
rstart <- rast(d, nlyr=1)
nc <- ncell(d) 
set.seed(88)
values(d) <- t(sapply(1:5, function(i) runif(nc, min = -10*i, max = 10)))
values(rstart) <- sample(2, nc, replace=TRUE)
rend <- rstart + 3

解决办法 

idx <- c(rstart, rend)
z <- rapp(d, idx, "sum")
z  
#class       : SpatRaster 
#dimensions  : 100, 100, 1  (nrow, ncol, nlyr)
#resolution  : 3.6, 1.8  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#data source : memory 
#names       :      lyr1 
#min values  : -184.6918 
#max values  :  34.93876
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.