更有效地计算栅格数据上的 Thornthwaite 蒸散量

Question

我目前正在对巨大的栅格堆栈运行 Thornthwaite 蒸散量计算。我正在关注另一个堆栈溢出线程中 Robert Hijmans 的解决方案：Thornthwaite evapotranspiration on a raster dataset 。误差公式未矢量化

所述线程中的简短可重现示例如下：

library(SPEI)
library(raster)
library(zoo)

tm = array(20,c(3,4,12*64))
tm = brick(tm)

dates=seq(as.Date("1950-01-01"), as.Date("2013-12-31"), by="month")
tm<- setZ(tm,dates)
names(tm) <- as.yearmon(getZ(tm))

#thornthwaite ET
th <- function(Tave, lat) {
    as.vector(SPEI::thornthwaite(as.vector(Tave), lat))
} 

a <- raster(tm)
lat <- init(a, "y")
#out <- raster::overlay(tm, lat, fun = th)
out <- brick(tm, values=FALSE)

for (i in 1:ncell(tm)) {
    out[i] <- th(tm[i], lat[i])
}

然而，问题是 Thornthwaite PET 的逐单元计算似乎效率相当低，因此我努力寻求使用覆盖或计算的解决方案，我认为这可以提高计算效率。在前面提到的 stackoverflow 线程中，引入了另一个使用覆盖的解决方案 - 但是，我无法重现它。

任何可能提高计算效率的解决方案都将受到欢迎！

Answer 1

我想说一种或另一种方法的效率取决于每个网格是否有更多单元或更多层，基本上分别是您感兴趣的区域/空间分辨率以及您正在研究的周期/时间分辨率的函数。 .. 两种方法可能都是有效的，但这完全取决于您的情况，从我的角度来看。

但是，在您提到的问题中，@mastefan 似乎非常有信心他能够修复 @Robert Hijmans 的解决方案。但这是在2018年回答的，我也无法重现这个问题。所以我想在 CRAN 上尝试 2018 版本的

{SPEI}

（也许我还应该使用旧版本的

{raster}

和

{zoo}

来保持结果一致，但你现在明白了......）：

url <- "https://cran.r-project.org/src/contrib/Archive/SPEI/SPEI_1.7.tar.gz"
install.packages(url, repos = NULL, type = "source")

让我们运行代码（请原谅我没有重复初始化和预处理步骤的所有步骤）：

raster::overlay(tm, lat, fun = Vectorize(th))
#> class      : RasterBrick 
#> dimensions : 3, 4, 12, 768  (nrow, ncol, ncell, nlayers)
#> resolution : 0.25, 0.3333333  (x, y)
#> extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#> crs        : NA 
#> source     : memory
#> names      :  layer.1,  layer.2,  layer.3,  layer.4,  layer.5,  layer.6,  layer.7,  layer.8,  layer.9, layer.10, layer.11, layer.12, layer.13, layer.14, layer.15, ... 
#> min values : 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, ... 
#> max values : 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, ...

太好了，现在可以执行了。看起来

SPEI::thornthwaite()

的行为自 2018 年以来发生了变化。但是让我们看看执行速度，因为你担心效率：

mbm <- microbenchmark::microbenchmark("raster::overlay()" = raster::overlay(tm, lat, fun = Vectorize(th)),
                                      "cell-wise loop" = for (i in 1:ncell(tm)) {
                                        
                                        out[i] <- th(tm[i], lat[i])
                                      }, 
                                      times = 100)
mbm
#> Unit: milliseconds
#>               expr       min        lq     mean    median        uq      max neval
#>  raster::overlay() 7694.4259 7865.7535 8040.483 7995.6493 8126.7056 9043.048   100
#>    cell-wise loop  102.5918  106.5566  116.190  114.4493  121.6464  171.598   100

这看起来一点也不快，恰恰相反。但公平地说，我对

{raster}

不太熟悉，所以也许更熟练的人可以很容易地改进它。我只是想按原样分析这些片段。想必

{terra}

将成为 2024 年的首选方案。

然而，更让我担心的是两种方法的不同结果。也许我们应该更关心获得正确的结果而不是尽可能快。

out1 <- raster::overlay(tm, lat, fun = Vectorize(th))
out1
#> class      : RasterBrick 
#> dimensions : 3, 4, 12, 768  (nrow, ncol, ncell, nlayers)
#> resolution : 0.25, 0.3333333  (x, y)
#> extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#> crs        : NA 
#> source     : memory
#> names      :  layer.1,  layer.2,  layer.3,  layer.4,  layer.5,  layer.6,  layer.7,  layer.8,  layer.9, layer.10, layer.11, layer.12, layer.13, layer.14, layer.15, ... 
#> min values : 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, 125.1804, ... 
#> max values : 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, 125.5248, ...

out2 <- raster::brick(tm, values = FALSE)

for (i in 1:ncell(tm)) {
  
  out2[i] <- th(tm[i], lat[i])
}

out2
#> class      : RasterBrick 
#> dimensions : 3, 4, 12, 768  (nrow, ncol, ncell, nlayers)
#> resolution : 0.25, 0.3333333  (x, y)
#> extent     : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#> crs        : NA 
#> source     : memory
#> names      :  layer.1,  layer.2,  layer.3,  layer.4,  layer.5,  layer.6,  layer.7,  layer.8,  layer.9, layer.10, layer.11, layer.12, layer.13, layer.14, layer.15, ... 
#> min values : 76.06896, 68.79146, 76.29749, 73.89062, 76.37712, 73.92394, 76.38355, 76.36436, 73.87383, 76.21108, 73.64047, 76.04259, 76.06896, 68.79146, 76.29749, ... 
#> max values : 76.27825, 68.91329, 76.32396, 73.97994, 76.56333, 74.14654, 76.59547, 76.49955, 73.89596, 76.30668, 73.82273, 76.27298, 76.27825, 68.91329, 76.32396, ...