[Bfast-commits] running bfast on a 3D image

[Cheney Shreve]

Hi,

 I have a question regarding running bfast on 3D images in R (for
example, a MODIS time series). Using caTools, I’ve found a method for
easily reading and writing ENVI files [I’m using 16 bit int., bil].
Based on your previous notes to Brian, I’ve learned more about how to
extract various components of bfast. I’ve figured out a brute force
method (very likely an ugly, inefficient one, because I’m not used to
object oriented programming) for turning a pixel into a time series
object and running bfast. I’m having trouble storing the output in a
2D matrix. What I’d like to have is an output of the timing of the
major disturbance per pixel. So, if my input array is [1200,1200,400],
I’d have an output of [1200,1200] with time of the first disturbance.
My three primary questions are: (1) is there a better way to run bfast
over all the pixels in an image (avoiding for loops), (2) how to fill
the output array with the results, and (3) I get the following
warning: In sqrt(fr) : NaNs produced when I run bfast for some pixels
(there are no NaN values in the actual time series). I've searched the
R help for multidimensional matrices, but haven't found much
applicable information.

What I have so far:

require(bfast)
require(caTools)
require(tseries)

loop.test<-function(infile){
    data=read.ENVI(infile)
#infile is just a test (subset) array of time series data  [20,20,383]	

	#find the dimensions of input image
	x<-ncol(data)
	y<-nrow(data)
	z<-dim(data)[3]
	
     #set up an output array to hold results
    output.matrix<-matrix(NA,x,y)
    for (i in 1:x){
    for (j in 1:y) {
	#isolate the time series for a single pixel
	pixel<-data[i,j,1:z]
 	gpp<-ts(pixel,start=c(2002,23),end=c(2010,38),frequency=46)
    fit<-bfast(gpp,h=0.15,max.iter=1,breaks=NULL)
    niter <- length(fit$output)
    out <- fit$output[[niter]]

     if (out$Vt.bp[1]>0) {
       (tnrbps <- time(gpp)[(out$Vt.bp)])
        majorbps<-tnrbps[1]
      } else {
      	
      	majorbps=0     	
      	}

     #I know this is wrong, but don’t know what my sin is.
    output.matrix[i,j]<-as.vector(majorbps)
   # I also tried just output.matrix[i,j]=majorbps but got an error
about dimensions being incorrect. It doesn't make sense to me that
  # it would need to be a vector, but for some reason using, as.vector
doesn't throw an error?

     }
}
write.ENVI(output.matrix, “output.bil”, interleave=c(“bil”))
}


Thanks in advance for your help and patience.

Kind regards,

Cheney