[Bfast-commits] FW: Running bfast on raster bricks

[White, Joanne]

Hi, Jan.

Thanks for this.

I tried your code and the function does work on the data extracted from the raster brick outside the calc() function. So I'm guessing this means the issue is with calc() function? I have no idea why, so perhaps I should post a question to the raster list? Ultimately I would like to run this over a larger dataset that likely cannot be fully read into memory for processing.

The test image I am using is very small (4 pixels!) and I am happy to share if it would be useful. If you have an ftp site, I could post it there for you. Please let me know.

In terms of output, I would *like* two raster images that have the time and magnitude for each pixel in the image. I know that multi-argument returns are not permitted, so does this mean I have to run two functions (once and if I get it working with calc), one for magnitude and one for time?

Thanks again,

Joanne   

"test.img"
class       : RasterBrick 
dimensions  : 2, 2, 4, 828  (nrow, ncol, ncell, nlayers)
resolution  : 1000, 1000  (x, y)
extent      : 787641.7, 789641.7, 5455728, 5457728  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
values      : T:\ecomonitor\age\prototype\data\fpar\test.img 


-----Original Message-----
From: Jan Verbesselt [mailto:janverbesselt at gmail.com] 
Sent: March-19-12 2:47 PM
To: White, Joanne
Cc: bfast-commits at lists.r-forge.r-project.org
Subject: Re: [Bfast-commits] FW: Running bfast on raster bricks

Hi Joanne,

yes. t() transpose should be in the function - in orde to obtain the data in the right format.
would it be possible to share this test data set e.g. via dropbox so that I can test run your function on it? If it works I could add it as a demo data set to the bfast package :-).
what might help to find the error and debug is to try the function on the data extracted from the raster brick before using it within the
calc() function:

## I have not tried this so might also not work immediately xbfast <- function(data) {
          ndvi <- ts(data, frequency=36, start=c(1985,1))
         result <- bfast(ndvi, h=0.012, season="harmonic", max.iter=1, breaks=2)
        result$Time
 }

bfastfun <- function(x) {
 output <- t(apply(x, 1, xbfast))
 return(output)
}

y = values(data)
bfastfun(y)

Hope this helps,
Jan



On Mon, Mar 19, 2012 at 6:38 PM, White, Joanne <Joanne.White at nrcan-rncan.gc.ca> wrote:
>
> Hello,
>
> I am attempting to run bfast on a raster brick (a time series of AVHRR 
> data).
>
> I have adapted the tutorial code provided in the bfast manual for 
> applying bfastmonitor to raster bricks (code below).
>
> When I try to apply my bfast function to a full image, I get the 
> following
> error:
>
>         Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
>                  cannot use this function
>
> The test of the single pixel works fine and returns the correct value.
>
> I would appreciate any suggestions on how to fix this.
>
> Many thanks,
>
> Joanne
>
>
> ################################################################
> require("raster")
> require("bfast")
> setwd("T:/ecomonitor/age/prototype/data/fpar")
>
> #read the ENVI format file into a raster brick data <- 
> brick("test.img")
>
> #read in the list of dates for the time series date.list <- 
> scan("dates_list.txt", what=character(0))
>
> #assign layer names from date.list
> layerNames(data) <- date.list
>
> #helper function for the calc function xbfast <- function(data, 
> date.list) {
>          ndvi <- ts(data, names=date.list, frequency=36, 
> start=c(1985,1),
> end=c(2007,36))
>         result <- bfast(ndvi, h=0.012, season="harmonic", max.iter=1,
> breaks=2)
>     return(cbind(result$Time))
> }
>
> #apply on a single pixel for testing
> pixel <- as.vector(data[4])
> ndvi <- ts(pixel, names=date.list, frequency=36, start=c(1985,1),
> end=c(2007,36))
> output <- xbfast(ndvi, date.list)
> plot(output)
>
> #apply on the full image                                    #This is 
> the part that throws the error timeofbreak <- calc(data, 
> fun=function(x) { output <- ts(apply(x, 1, xbfast, date.list))     
> #Not sure if this is supposed to be "t" as in the example or "ts". 
> Regardless, the error message is the same.
> return(output)
> })
> plot(timeofbreak)
> ################################################################
>
>
> Note that the test.img is very small (4 x 4 pixels):
>>data
> class       : RasterBrick
> dimensions  : 2, 2, 4, 828  (nrow, ncol, ncell, nlayers) resolution  : 
> 1000, 1000  (x, y) extent      : 787641.7, 789641.7, 5455728, 5457728  
> (xmin, xmax, ymin, ymax) coord. ref. : +proj=utm +zone=16 +datum=WGS84 
> +units=m +no_defs +ellps=WGS84
> +towgs84=0,0,0
> values      : T:\ecomonitor\age\prototype\data\fpar\test.img
> layer names : X1985.01.01, X1985.01.11, X1985.01.21, X1985.02.01, 
> X1985.02.11, X1985.02.21, X1985.03.01,...
>
>
>
>
>
> _______________________________________________
> Bfast-commits mailing list
> Bfast-commits at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/bfast-com
> mits