R：将栅格图层作为一个对象读取以在每个 ENVI 文件中应用函数答案

【问题标题】：R: Read raster layers as one object to apply a function in each ENVI fileR：将栅格图层作为一个对象读取以在每个 ENVI 文件中应用函数
【发布时间】：2017-06-02 18:35:17
【问题描述】：

我有一个 ENVI 文件 (82_83_test.envi)，其中包含从 1982 年到 1983 年的双周栅格图层。即每年 24 层，总共 48 层。我想创建一个 for 循环来应用一个函数来每年执行时间序列分析，即 R 将在一个像素中运行 24 层，并使用函数“fun”为该年的所有像素计算 5 个参数。最终，我希望每年有 5 个地块（5 个参数），所以两年总共有 10 个地块。

我尝试在每个文件中使用 1 个具有 2 年数据的 ENVI 文件和 2 个具有 1 年数据的 ENVI 文件。我使用了库 spatial.tools 中的brickstack_to_raster_list() 来读取文件，我得到了 48 层。但是，我想获得 2 个块（1982 年和 1983 年），每个块由 24 层组成，以便我可以运行方程式。

也许像brickstack_to_raster_list() 之类的东西然后将第1 层到第24 层合并为一个，然后将第25 到第48 层合并为一个？

new <- stack("82_83_test.envi")
new1<- brickstack_to_raster_list(new)

new1 返回 48 个栅格图层。例如，

new1
[[1]]
class       : RasterLayer 
band        : 1  (of  48  bands)
dimensions  : 151, 101, 15251  (nrow, ncol, ncell)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -105.0833, -96.66667, 56.66667, 69.25  (xmin, xmax, ymin, 
ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 
+towgs84=0,0,0 
data source : C:\*\82_83_test.envi 
names       : Band.1 
values      : -32768, 5038  (min, max)

另一种方法是将多个年度 ENVI 文件连接到一个列表中。

new <- stack("1982_test.envi")
new1<- stack(new,new)
new2<- brickstack_to_raster_list(new1)

上述两种方法产生相同的结果，尽管我不确定它的效率。因为设置完之后，我会生成 1982 年到 2015 年的数据，所以效率很重要。

下面是我想在 for 循环中应用的函数。

# A is an unknown that will be the number of components in the list. 
for (i in length(A)) {
new1[new1<=-1000]<-0
Data_value<-new1/10000
# assign 0 to pixel value that is less than -1000 and divide by 10000 in order to use the equation
DOY<-(1:nlayers(new1)*15)
# so that the unit will be in days instead of the number of weeks.

fun<- function(x) { if (all(is.na(x[1]))) { return(rep(NA,5)) } else { 
fitForThisData  <-nls(x~ a+((b/(1+ exp(-c*(DOY-e))))- (g/(1+ exp(-d*(DOY-
f))))), alg="port",start=list(a=0.1,b=1,g=1,c=0.04,d=0.04,e=112,f=218),
lower=list(a=0,b=0.3,g=0.3,c=-1,d=-1,e=20,f=100),
upper=list(a=0.4,b=2,g=2,c=1,d=1,e=230,f=365),
control=nls.control(maxiter=2000, tol = 1e-15, minFactor = 1/1024, 
warnOnly=TRUE))
SOS<-(coef(fitForThisData)[6] -(4.562/(2*coef(fitForThisData)[4])))
EOS<-(coef(fitForThisData)[7] -(4.562/(2*coef(fitForThisData)[5])))
LOS<-(EOS-SOS)
SPUDOY<-(1.317*((-1/coef(fitForThisData)[4])+ coef(fitForThisData)[6]))
P_TAmplitude<-(SPUDOY-SOS)
return (c(SOS,EOS,LOS,SPUDOY,P_TAmplitude))
}
}
}

equation<-calc(Data_value,fun,forceapply=TRUE)
plot(equation)

非常感谢您对如何执行此操作的建议。非常感谢你！

【问题讨论】：

标签： r spatial raster r-raster

【解决方案1】：

阅读你的堆栈后：

library(raster)    
new <- stack("82_83_test.envi")

只需使用基本索引将堆栈拆分为年度子堆栈：

year1 <- new[[1:24]]
year2 <- new[[25:48]]

【讨论】：

索引听起来是个不错的方法。做一个嵌套的for循环怎么样？或者还有其他更有效的方法吗？
嗨@maRtin，我尝试了索引，但我不能让 writeRaster 为不同年份的我编写两个单独的文件。 year1
@Candice 抱歉，但很难复制和理解问题所在。一个最小的可重现示例：stackoverflow.com/questions/5963269/… 将非常有帮助。也许您可以消除有效的代码并构建一个关于无效部分的简短示例以及输出的外观。如果您遵循这些准则，我相信您的问题会吸引更多答案。
示例应包括以下内容：-) 重现错误所需的最小数据集 -) 重现错误所需的最小可运行代码，可在给定数据集上运行。

【解决方案2】：

更新：我能够循环该函数，但我的猜测是在将其与真值进行比较后，计算是在所有栅格图层上完成的。但是，由于两个文件具有相同的摘要，所以会生成两个相同内容但文件名不同的文件。

new <- stack("82_83_test.envi")
new[new<=-1000]<-0
Data_value<-new/10000
nlayers <- nlayers(new)
nyears <- nlayers(new)/24
DOY<-((1:nlayers(new))/nyears)*15
dummy<- FALSE

for (i in 1:nyears) {
  for (j in (1+24*(i-1)):(24*i)) { 
    fun<-function (x)
    equation<-calc(Data_value,fun,forceapply=TRUE)
    date<- 1981+i
    writeRaster(equation,filename=paste("Output",date,".envi",sep=""),
    format="ENVI",overwrite=T)
    if (j == nlayers){
    dummy<-TRUE
    break
    if (dummy) {break}      
}
}

【讨论】：