为从目录中读取光栅、将数据汇总到单个光栅并输出新光栅而创建的循环的索引不起作用

Question

我有几个目录，里面都是每日的气候数据。我需要将每日栅格合并为每周栅格，一些是通过值的总和，一些是通过值的平均值。到目前为止，我已经在目录（其中包含每日光栅文件）中创建了一个文件名向量，并为编写了一个循环来导入前7个光栅，将光栅放入一个堆栈，使用计算将堆栈中光栅的平均值（或总和）写入一个新光栅，并将光栅写入工作目录中的一个文件。然后，我尝试从向量中删除前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:\dailyPrism\10843\prcp\1980_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

下面是如何通过从文件名向量创建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)
}

Answer 2

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

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）。

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

易于修改以解决问题（无需优化代码）：

修改循环：

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))