index not working on for loop created to read rasters from a directory, summarize data into single raster, 并输出新的栅格

Question

我有几个充满栅格的目录，这些栅格是每日气候数据。我需要将每日栅格合并为每周栅格，一些按值的总和，一些按值的平均值。到目前为止，我已经在目录中创建了一个文件名向量（其中包含每日栅格文件）并编写了一个 for 循环来导入前 7 个栅格，将栅格放入 stack , 使用 calc 将 stack 中栅格的平均值（或总和）转换为单个新栅格，并将栅格写入工作目录中的文件。然后，我通过尝试从向量中删除前 7 个名称并对向量中剩余的前 7 个文件名重复循环来处理新的文件名向量。我遇到的问题是第一个索引文件名没有从向量中删除。这是代码。

#file names for rasters are in a column of data frame
fname <- (repDf$fname)

#get rid of first 5 values to start on Sunday
fname <- fname[-c(1:5)]

#look at beginning of fname vector
head(fname)
[1] "1980_6.asc"  "1980_7.asc"  "1980_8.asc"  "1980_9.asc"  "1980_10.asc" "1980_11.asc"

for (i in seq_along(fname)){
  f1 <- fname[[i]]
  f2 <- fname[[i+1]]  
  f3 <- fname[[i+2]]
  f4 <- fname[[i+3]]
  f5 <- fname[[i+4]]
  f6 <- fname[[i+5]]
  f7 <- fname[[i+6]]
  r1 <- raster(f1)
  r2 <- raster(f2)
  r3 <- raster(f3)
  r4 <- raster(f4)
  r5 <- raster(f5)
  r6 <- raster(f6)
  r7 <- raster(f7)
  s <- stack(r1,r2,r3,r4,r5,r6,r7)
  r <- calc(s, fun=sum)
  r <- r * 0.0393701
  r <- round(r, 2)
  writeRaster(x=r, filename=paste0("week_", i, sep=""), format="ascii", overwrite=TRUE)
  remove <- c(f1,f2,f3,f4,f5,f6,f7)
  fname <- fname[! fname %in% remove]
}

#Example RasterLayer (after loop has run twice)
print(r1)
class       : RasterLayer 
dimensions  : 227, 199, 45173  (nrow, ncol, ncell)
resolution  : 994.9749, 994.9749  (x, y)
extent      : 367500, 565500, -1325500, -1099641  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : G:\dailyPrism843\prcp80_30.asc 
names       : X1980_30 

#Example indexed file names (after loop has run twice)
f1
[1] "1980_30.asc"

#Example output raster (after loop has run twice)
print(r)
class       : RasterLayer 
dimensions  : 227, 199, 45173  (nrow, ncol, ncell)
resolution  : 994.9749, 994.9749  (x, y)
extent      : 367500, 565500, -1325500, -1099641  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : in memory
names       : layer 
values      : 1.02, 3.54  (min, max)
The problem can be seen here;

head(fname)
[1] "1980_13.asc" "1980_21.asc" "1980_29.asc" "1980_30.asc" "1980_31.asc" "1980_32.asc"

出于某种原因，将用作一周第一天的文件名保留在 fname 向量中，并且进入每周计算的文件不代表他们需要的天数。任何帮助是极大的赞赏。我会尝试提供一些示例文件以供使用，但这些文件非常大。

Answer 1

在您的代码中添加一些调试代码，您可以识别发生了什么（我已禁用处理部分，因为我没有文件或计算逻辑）：

fname <- c(paste0("1980_", 6:40, ".asc"))
fname[1]
seq_along(fname) # vector from 1 to length of vector!!!

for (i in seq_along(fname)){
  print(i)
  print(paste("Size of fname:", length(fname)))
  print(head(fname))
  print(fname[i])
  f1 <- fname[[i]]
  f2 <- fname[[i+1]]  
  f3 <- fname[[i+2]]
  f4 <- fname[[i+3]]
  f5 <- fname[[i+4]]
  f6 <- fname[[i+5]]
  f7 <- fname[[i+6]]
#   r1 <- raster(f1)
#   r2 <- raster(f2)
#   r3 <- raster(f3)
#   r4 <- raster(f4)
#   r5 <- raster(f5)
#   r6 <- raster(f6)
#   r7 <- raster(f7)
#  s <- stack(r1,r2,r3,r4,r5,r6,r7)
#  r <- calc(s, fun=sum)
#  r <- r * 0.0393701
#  r <- round(r, 2)
#  writeRaster(x=r, filename=paste0("week_", i, sep=""), format="ascii", overwrite=TRUE)
  remove <- c(f1,f2,f3,f4,f5,f6,f7)
  fname <- fname[! fname %in% remove]
}

这将导致：

[1] 1
[1] "Size of fname: 35"
[1] "1980_6.asc"  "1980_7.asc"  "1980_8.asc"  "1980_9.asc"  "1980_10.asc" "1980_11.asc"
[1] 2
[1] "Size of fname: 28"
[1] "1980_13.asc" "1980_14.asc" "1980_15.asc" "1980_16.asc" "1980_17.asc" "1980_18.asc"
[1] 3
[1] "Size of fname: 21"
[1] "1980_13.asc" "1980_21.asc" "1980_22.asc" "1980_23.asc" "1980_24.asc" "1980_25.asc"
[1] 4
[1] "Size of fname: 14"
[1] "1980_13.asc" "1980_21.asc" "1980_29.asc" "1980_30.asc" "1980_31.asc" "1980_32.asc"
[1] 5
[1] "Size of fname: 7"
[1] "1980_13.asc" "1980_21.asc" "1980_29.asc" "1980_37.asc" "1980_38.asc" "1980_39.asc"
Error in fname[[i + 3]] : subscript out of bounds

原因是您遍历了预定义数量的 fname 向量项（在我的示例中 "seq_along" = 35）。

因此您删除了已处理的项目，但 i 并未重置为 1，而是在每次循环时递增。

而且，还有很多改进代码的空间（例如，不需要删除元素，只需循环向量元素；如何处理最后一个循环中错误的向量大小；为什么要使用 double方括号访问 fnames...)

轻松修改解决问题（无需优化代码）：

修改循环：

for (i in seq(1,length(fname), by=7)) {

删除两行：

remove <- c(f1,f2,f3,f4,f5,f6,f7)
fname <- fname[! fname %in% remove]

并改为添加调试输出：

print( paste("processing files", f1, "to", f7))

Answer 2

这是通过从文件名向量创建 RasterStack 来使代码更简洁的方法

fname <- repDf$fname[-c(1:5)]

for (i in seq(1,length(fname), by=7)){
    s <- stack(fname[i:(i+6)])
    r <- sum(s) * 0.0393701
    r <- round(r, 2)
    writeRaster(r, filename=paste0("week_", i), format="ascii", overwrite=TRUE)
}

index not working on for loop created to read rasters from a directory, summarize data into single raster, 并输出新的栅格

index not working on for loop created to read rasters from a directory, summarize data into single raster, and output new raster

for-loop

r

raster