优化栅格::计算函数 - 函数1 vs 2 - R

Question

我正在计算基于2个栅格（输入ras）和“层”栅格的新栅格（输出ras）。 Stratum栅格值（1到4）指的是偏置和权重数据框中的行。 Strata值'4'用于填充Strata栅格中的任何'NA'，否则该函数将崩溃。以下输入是必需的。

# load library 
library(raster) 

# reproducing the bias and weight data.frames 
bias <- data.frame(
ras_1 = c(56,-7,-30,0), 
ras_2 = c(29,18,-52,0), 
ras_3 = c(44,4,-15,0) 
) 
rownames(bias) <- c("Strat 1","Strat 2","Strat 3","Strat 4") 

weight <- data.frame(
ras_1 = c(0.56,0.66,0.23,0.33), 
ras_2 = c(0.03,0.18,0.5,0.33), 
ras_3 = c(0.41,0.16,0.22,0.34) 
) 
rownames(weight) <- c("Strat 1","Strat 2","Strat 3","Strat 4") 

以下函数（融合）允许我为输入栅格添加“偏差”值。在添加偏差之后，两个经过校正的输入栅格单元值将乘以权重值，具体取决于它们属于哪个层。

输入2栅格值的结果将被累加并使用'calc'返回。

## Create raster data for input 

# create 2 rasters 
r1 <- raster(ncol=10,nrow=10) 
r2 <- raster(ncol=10,nrow=10) 
r1[] <- sample(seq(from = 1, to = 500, by = 1), size = 100, replace = TRUE) 
r2[] <- sample(seq(from = 1, to = 500, by = 1), size = 100, replace = TRUE) 
r2[1:2] <- NA # include NA in input maps for example purpose 

# Create strata raster (4 strata's) 
r3 <- raster(ncol=10,nrow=10) 
r3[] <- sample(seq(from = 1, to = 4, by = 1), size = 100, replace = TRUE) 

Strata.n <- 4 # number of strata values in this example 

fusion <- function(x) { 
    result <- matrix(NA, dim(x)[1], 1) 
    for (n in 1:Strata.n) { 
    ok <- !is.na(x[,3]) & x[,3] == n 
    a <- x[ok,1] + bias[n,1] # add bias to first input raster value    
    b <- x[ok,2] + bias[n,2] # add bias to second input raster value 
    result[ok] <- a * weight[n,1] + b * weight[n,2] # Multiply values by weight 
    } 
    return(result) 
} 

s <- stack(r1,r2,r3) 
Fused.map <- calc(s, fun = fusion, progress = 'text') 

具有上述功能的问题是：

• 它仅适用于2个栅格

• 如果一个光栅具有NA，那么结果将是NA该小区

  is.na(Fused.map@data@values) # check for NA in the fused map 
  

我想要的是：

• 即需要功能的任何输入栅格的数目
• 它可以与NA值工作（忽略栅格NA值）
• 重新调整“权重”，如果一个光栅具有NA值，从而使剩余的权重值加起来1

EDIT

下面的函数WHA我需要，但是比上面的大栅格上的函数慢得多。融合做它在10秒内，fusion2功能下面需要在大栅格8小时...

fusion2 <- function(x) { 
    m <- matrix(x, nrow= 1, ncol=3) # Create matrix per stack of cells 
    n <- m[,3] # get the stratum 
    g <- m[1:(Strata.n-1)] + as.matrix(bias[n,]) # add bias to raster values 
    g[g < 0] <- 0 # set values below 0 to 0 
    w <- weight[n,1:(Strata.n-1)] # get correct strata weight values 
    w[is.na(g)]<- NA # set weight to NA if (g) raster values are NA 
    p <- sum(w, na.rm = T) # calculate sum of weight values 
    pp <- w/p # divide weight values by sum to get the proportion to == 1 
    pp <- as.numeric(pp) 
    result <- as.integer(round(sum(pp*g, na.rm = T))) # return raster value 
    return(result) 
} 

Fused.map <- calc(s, fun = fusion2, progress = 'text') 

任何方式优化fusion2功能的类似方法作为fusion1？

> sessionInfo() 
R version 3.3.2 (2016-10-31) 
Platform: x86_64-w64-mingw32/x64 (64-bit) 
Running under: Windows >= 8 x64 (build 9200)  

谢谢你的时间！

Answer 1

似乎有很多不必要的格式转换正在进行，使用最简单的数据结构是最快的。 calc参数是一个数字向量，因此您可以在任何地方使用数字向量。而且，舍入和转换为整数是多余的。

fusion3 <- function(x) { 
    n <- x[3] # get the stratum 
    g <- x[1:(Strata.n-1)] + as.numeric(bias[n,]) # add bias to raster values 
    g[g < 0] <- 0 # set values below 0 to 0 
    w <- as.numeric(weight[n,1:(Strata.n-1)]) # get correct strata weight values 
    w[is.na(g)]<- NA # set weight to NA if (g) raster values are NA 
    p <- sum(w, na.rm = T) # calculate sum of weight values 
    pp <- w/p # divide weight values by sum to get the proportion to == 1 
    result <- as.integer(sum(pp*g, na.rm = T)+0.5) # return raster value 
    return(result) 
} 

在一个100×100光栅，你原来的功能需要：

system.time(Fused.map <- calc(s, fun = fusion, progress = 'text')) 
    user system elapsed 
    0.015 0.000 0.015 
system.time(Fused.map <- calc(s, fun = fusion2, progress = 'text')) 
    user system elapsed 
    8.270 0.078 8.312 

修改后的功能已经是快5倍：从数据帧

system.time(Fused.map3 <- calc(s, fun = fusion3, progress = 'text')) 
    user system elapsed 
    1.970 0.026 1.987 

接下来，预先计算矩阵所以你不需要为每个像素做到这一点：

bias_matrix = as.matrix(bias) 
weight_matrix = as.matrix(weight) 

fusion3 <- function(x) { 
    n <- x[3] # get the stratum 
    g <- x[1:(Strata.n-1)] + bias_matrix[n,] # add bias to raster values 
    g[g < 0] <- 0 # set values below 0 to 0 
    w <- weight_matrix[n,1:(Strata.n-1)] # get correct strata weight values 
    w[is.na(g)]<- NA # set weight to NA if (g) raster values are NA 
    p <- sum(w, na.rm = T) # calculate sum of weight values 
    pp <- w/p # divide weight values by sum to get the proportion to == 1 
    result <- as.integer(sum(pp*g, na.rm = T)+0.5) # return raster value 
    return(result) 
} 

我们得到：

system.time(Fused.map3 <- calc(s, fun = fusion3, progress = 'text')) 
    user system elapsed 
0.312 0.008 0.318 

最后，还预先计算1:(Strata.n-1)：

bias_matrix = as.matrix(bias) 
weight_matrix = as.matrix(weight) 
Strata.minus1 = 1:(Strata.n-1) 

fusion3 <- function(x) { 
    n <- x[3] # get the stratum 
    g <- x[Strata.minus1] + bias_matrix[n,] # add bias to raster values 
    g[g < 0] <- 0 # set values below 0 to 0 
    w <- weight_matrix[n,Strata.minus1] # get correct strata weight values 
    w[is.na(g)]<- NA # set weight to NA if (g) raster values are NA 
    p <- sum(w, na.rm = T) # calculate sum of weight values 
    pp <- w/p # divide weight values by sum to get the proportion to == 1 
    result <- as.integer(sum(pp*g, na.rm = T)+0.5) # return raster value 
    return(result) 
} 

我们得到：

system.time(Fused.map3 <- calc(s, fun = fusion3, progress = 'text')) 
    user system elapsed 
    0.252 0.011 0.262 

这不太0.015，但你还必须考虑到您的原始函数不输出整数，也不会将值设置为低于0，也不会将比例总和设为1，也不会像您提到的那样与NAs一起。

你要知道，这个功能仍然只只有两个栅格的作品，因为你硬编码层数为3层，您应该改用raster::overlay有两个参数，地层栅格和图层本身（或使用calc层数光栅作为层1，但这不是calc的设计目的）。