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

Question

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

我尝试使用一个包含2年数据的ENVI文件和2个包含1年数据的ENVI文件。我使用库spatial.tools中的brickstack_to_raster_list()来读取文件，我得到了48层。但是，我想得到两个块(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)

我真的很感谢你对如何做到这一点的建议。非常感谢!

Answer 1

在堆栈中读取后：

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

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

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

Answer 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}      
}
}