R 中的 clusterApply 是否会影响每层栅格砖的平均值计算？

Question

我有来自 Daymet 的每日最低和最高温度的光栅砖。地理范围是大陆 U.S。每块砖有 365 层，一层代表一年中的每一天。我想从一对分钟计算平均每日温度。和最大。按年砖。 R raster 包中的 overlay() 函数似乎可以满足我的要求。例如：

library(raster)
tmin <- brick("tmin_2018.nc")
tmax <- brick("tmax_2018.nc")
tmean <- overlay(tmin, tmax, fun=mean)

问题是，这需要很长时间和大量 RAM。内存比我多。我可以访问高性能计算集群。我在使用 OpenMPI 实用程序并可以提供并行处理的集群上拥有此过程的代码。我的测试用例只使用了 5 层。我从集群中请求 6 个节点来完成这项工作。

library(iterators)
library(foreach)
library(snow)
library(parallel)
library(Rmpi)
library(compiler)
library(raster)
library(rgdal)
library(ncdf4)
print(mpi.universe.size())

nras <- brick("infil/ras_in/daymet_v3_tmax_2009_na.nc4")
xras <- brick("infil/ras_in/daymet_v3_tmin_2009_na.nc4")

cl <- makeCluster( mpi.universe.size() - 1, type='MPI' )

clusterExport(cl, "nras")
clusterExport(cl, "xras")
clusterExport(cl, c('overlay', 'writeRaster', 'brick'))

easy_fn <- function(i) {
        overlay(nras[[i]], xras[[i]], fun = mean)
        }

out <- clusterApply(cl, 1:5, fun=easy_fn)
outbrickb <- brick(out)
writeRaster(outbrickb, "outfil/rasout/testbrick.tif", wopt=list(datatype="FLT4S", filetype="GTiff"), overwrite=T)

这也行得通。但是价值观不同。当我计算一层（天）的平均值时，使用上面的第一个代码块，我得到的值介于 -55.255 和 31.08 之间。在使用第二个代码块对 HPC 集群进行运行少量层测试后，我得到的值介于 -50 和 44.5 之间。如果事情是按小数位计算的，我不会担心，但这些值有很大的不同。在 clusterApply() 或我的小 easy_fn() 工作方式与仅调用 overlay() 会导致不同值的方式中，我是否遗漏了什么？

Answer 1

使用 terra 可能比 raster

更简单、更快捷

示例数据

library(terra)
tmin <- rast(system.file("ex/logo.tif", package="terra")) / 10
tmax <- tmin + 10

解决方案

tavg1 <- mean(tmin, tmax)

请注意，这在 terra 中有效，因为它计算平行平均值（按层），而 raster 将 return 所有层的平均值。要用 terra 得到它，你会做

m <- mean(c(tmin, tmax))

raster 中 overlay 的等价物是 terra::lapp，在更复杂的情况下可能需要它。

tavg2 <- lapp(sds(tmin, tmax), function(x,y) (x + y)/2)

all(values(tavg1) == values(tavg2))
#[1] TRUE

现在 raster

示例数据

library(raster)
tmin <- brick(tmin)
tmax <- brick(tmax)

使用代数方法可能比 overlay

更快

tavg3 <- (tmin + tmax) / 2
tavg4 <- overlay(tmin, tmax, fun = mean)

用你的函数

easy_fn <- function(i) {
    overlay(tmin[[i]], tmax[[i]], fun = mean)
}

x <- lapply(1:nlayers(tmin), easy_fn)
tavg5 <- stack(x)

现在并行

library(parallel)
cl <- makeCluster(getOption("cl.cores", 2))
clusterExport(cl, c('overlay', 'tmin', 'tmax'))
out <- clusterApply(cl, 1:nlayers(tmin), fun=easy_fn)
stopCluster(cl)
tavg6 <- stack(out)

all(values(tavg1) == values(tavg3))
#[1] TRUE
all(values(tavg1) == values(tavg4))
#[1] TRUE
all(values(tavg1) == values(tavg5))
#[1] TRUE
all(values(tavg1) == values(tavg6))
#[1] TRUE

所有结果都一样。

现在有了更大的数据集，并为不同的方法计时

fmx <- "daymet_v4_daily_hi_tmax_2015.nc"
fmn <- "daymet_v4_daily_hi_tmin_2015.nc"
library(terra)
tmin <- rast(fmn)
tmax <- rast(fmx)
system.time(tavg1 <- mean(tmin, tmax))
#   user  system elapsed 
#   1.88    0.75    2.62 

system.time(tavg2 <- lapp(sds(tmin, tmax), function(x,y) (x + y)/2))
#   user  system elapsed 
#   2.34    1.40    3.76

这些显然是最快的方法（见下文）。您可以将文件名参数添加到 lapp，但是对于 mean，您还必须使用 writeRaster

和 raster

library(raster)
tmin <- brick(fmn)
tmax <- brick(fmx)
system.time(tavg3 <- (tmin + tmax) / 2)
#   user  system elapsed 
#  13.94    4.42   18.47

在这种情况下 overlay 非常慢

system.time(tavg4 <- overlay(tmin, tmax, fun = mean))
#  user  system elapsed 
# 478.06    5.78  488.16 

system.time(x <- lapply(1:nlayers(tmin), easy_fn))
#   user  system elapsed 
# 474.03    5.14  482.42 
tavg5 <- stack(x)

并行化确实更快；但还不足以让它值得追求

library(parallel)
cl <- makeCluster(getOption("cl.cores", 2))
clusterExport(cl, c('overlay', 'tmin', 'tmax'))
system.time(out <- clusterApply(cl, 1:nlayers(tmin), fun=easy_fn))
#   user  system elapsed 
#   0.53    0.78  184.82 
stopCluster(cl)
tavg6 <- stack(out)

我看不出顺序方法和并行方法之间有什么有意义的区别

d <- tavg6 - tavg4
d
#class      : RasterBrick 
#dimensions : 584, 284, 165856, 365  (nrow, ncol, ncell, nlayers)
#resolution : 1000, 1000  (x, y)
#extent     : -5802750, -5518750, -622500, -38500  (xmin, xmax, ymin, ymax)
#crs        : +proj=lcc +lat_0=42.5 +lon_0=-100 +lat_1=25 +lat_2=60 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs 
#source     : r_tmp_2021-11-20_103534_54328_41274.grd 
#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 : -9.536743e-07, -9.536743e-07, -9.536743e-07, -9.536743e-07, -9.536743e-07, -9.536743e-07, -9.536743e-07, -9.536743e-07, -9.536743e-07, -9.536743e-07, -9.536743e-07, -9.536743e-07, -9.536743e-07, -9.536743e-07, -9.536743e-07, ... 
#max values :  9.536743e-07,  9.536743e-07,  9.536743e-07,  9.536743e-07,  9.536743e-07,  9.536743e-07,  9.536743e-07,  9.536743e-07,  9.536743e-07,  9.536743e-07,  9.536743e-07,  9.536743e-07,  9.536743e-07,  9.536743e-07,  9.536743e-07, ...

R 中的 clusterApply 是否会影响每层栅格砖的平均值计算？

Does clusterApply in R affect calculation of average for each layer of raster brick?

parallel-processing

r

raster