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

[Jan Verbesselt]

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