沿时间维度组合多个 spatRaster 文件时出错

Question

我正在处理来自 PRISM 的温度数据。到目前为止，我已经能够下载个别日期的数据，然后裁剪为形状文件，并相应地提取信息。 然而，我希望沿着时间维度组合多个 spatRaster 文件，所以我最终得到一个数据帧。 在此处链接我的两个 .rds 文件

到目前为止，我已经尝试了两种方法。 首先，我尝试将两天的数据合并到一个 .rds 文件（stars 对象）中，然后转换为 spatRaster：

file1<- readRDS("file1.rds")
file2<- readRDS("file2.rds")
combinedfiles<- c(file1, file2, along = 3)

合并的文件具有以下描述/属性：

stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
                                  Min. 1st Qu. Median    Mean 3rd Qu.  Max.  NA's
PRISM_tmin_stable_4kmD2_202...  -1.447    8.65 10.344 11.7053  16.168 20.97 60401
dimension(s):
        from   to offset    delta refsys x/y
x          1 1405   -125  0.04167  NAD83 [x]
y          1  621  49.94 -0.04167  NAD83 [y]
new_dim    1    2     NA       NA     NA    

但是，我似乎无法成功地将其转换为 spatRaster：

temprast<-stack(combinedfiles)
Error in stack.default(combinedfiles) : 
  at least one vector element is required

我还尝试将每一天转换为 spatRaster，然后组合 - 这也不起作用：

temprast<-rast(file1)
temprast2 <-rast(file2)

raster_brick <- brick(temprast, temprast2)
Error in .local(x, ...) : 
  unused argument (new("SpatRaster", pnt = new("Rcpp_SpatRaster", .xData = <environment>)))

最终目标是能够将此光栅文件裁剪为形状文件。如果我首先按时间顺序组合 .rds 文件，然后再处理形状文件裁剪等，那么会容易得多。但如果有必要，我想我的最后一个选择是将每个 .rds 文件单独裁剪为形状文件，提取数据，然后合并（可能使用bind_rows）？但这似乎效率很低。

Answer 1

PRISM 文件的结构均相同：621 行和 1405 列的原始二进制 32 位数字数据（BIL 格式）。这是 GIS 软件通常使用的格式（很久以前），

terra

和

stars

包都可以轻松处理该格式。然而，这两个包都会产生大量开销，并且对文件名中包含的日期信息没有良好的支持。由于您遇到了错误，因此您最好使用

readBin()

使用基本的基本 R 方法，然后自己设置一些附加信息。合并所有文件后，您可以在

terra

或

stars

中导入组合数据，并进行裁剪和进一步分析。

作为一个函数，合并目录中的所有数据文件将如下所示：

mergePRISM <- function(dir) {
  # Need package abind
  requireNamespace("abind")
    
  # The files to merge
  lf <- list.files(dir, "\\.bil$", full.names = TRUE)
  
  # The dates of each of the files
  dates <- sapply(strsplit(basename(lf), "_"), function(p) p[5])
  
  # Read the data, then rearrange to go from row-major (BIL) to column-major (R)
  data <- lapply(lf, function(fn) {
    d <- readBin(fn, what = "numeric", n = 621 * 1405, size = 4)
    dim(d) <- c(1405, 621)
    aperm(d, c(2, 1))
  })
  
  # Merge the data
  data <- abind(data, along = 3)
  
  # Set the dimension names
  lon <- seq(from = -125, to = -125 + 1404 * 0.0416666666667, by = 0.0416666666667)
  lat <- seq(from = 49.9166666666664, to = 49.9166666666664 - 620 * 0.0416666666667, by = -0.0416666666667)
  dimnames(data) <- list(lat, lon, dates)
  data
}

请注意，如果您的数据出现乱码，您可能必须设置

endian

的

readBin()

参数（在我的 Macbook 上不需要）。

您可以对一个月甚至一年的数据执行此操作，但由于 RAM 问题，我不会使时间序列更长。

您可以使用

stars

将结果转换为

st_as_stars()

对象，然后进行进一步处理，例如裁剪为 shapefile。

Answer 2

您正在使用

stack

想要使用

raster::stack

您正在使用

raster::brick

，您应该使用“terra”中的函数，即

c()

一个明智的方法可能是

ff <- list.files(pattern=".bil$")
library(terra)
x <- terra::rast(ff)