[Bfast-commits] FW: Running bfast on raster bricks
White, Joanne
Joanne.White at NRCan-RNCan.gc.ca
Thu Mar 22 19:13:32 CET 2012
Hi, Jan.
The code worked, with only one minor problem:
The last line:
plot(magntime)
Gives the following error:
Error in as.double(y) :
cannot coerce type 'S4' to vector of type 'double'
I removed the plot function and added the following to output an ENVI format file:
writeRaster(magntime, filename="output_bfast", format="ENVI", overwrite=TRUE)
Many thanks for your help with this.
Joanne
-----Original Message-----
From: Jan Verbesselt [mailto:janverbesselt at gmail.com]
Sent: March-22-12 9:58 AM
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,
Hope this works. Please try the code below out. Multi-arguments are possible now :-) so you can save any number of bfast output away to a file. In the following example the time and magnitude of the biggest break are returned in two rasters. It all takes quite some time so a multi-core computer, or a server might help. Let me know if anything about the bfast manual can be improved or if you have questions.
cheers,
Jan
## SCRIPT FOR Joanne
require("raster")
require("bfast")
## read the ENVI format file into a raster brick temp <- tempdir() download.file("http://dl.dropbox.com/u/8150174/test.img",
paste(temp,"/test.img",sep=""))
download.file("http://dl.dropbox.com/u/8150174/test.hdr",
paste(temp,"/test.hdr",sep=""))
data <- brick(paste(temp,"/test.img",sep=""))
## helper function for the calc function xbfast <- function(data) {
ndvi <- ts(data, frequency=36, start=c(1985,1))
result <- bfast(ndvi, season="harmonic", max.iter=1, breaks=2)
return(cbind(result$Magnitude,result$Time)) ## save magnitude and time of the biggest break }
## apply on a single pixel for testing
pixel <- as.vector(data[4])
ndvi <- ts(pixel, frequency=36, start=c(1985,1))
plot(ndvi)
result <- bfast(ndvi, season="harmonic", max.iter=1, breaks=2)
plot(result)
result$Magnitude # of the biggest break
result$Time # time of the biggest break
## apply on a single pixel using the xbfast function output <- xbfast(pixel) output
## optimise function that takes into account the percentage of NA's within a time series bfastfun <- function(y) {
percNA <- apply(y, 1, FUN=function(x) (sum(is.na(x))/length(x)) ) ## checks the percentage of NA's with the time series
i <- (percNA<0.2)
res <- matrix(NA, length(i), 2)
if (sum(i) > 0) {
res[i,] <- t(apply(y[i,], 1, xbfast))
}
res
}
#apply on the full satellite image time series magntime <- calc(data, fun=bfastfun)
layerNames(magntime) <- c("Magnitude", "Time of biggest Break")
plot(magntime)
## two rasters are returned: i.e. the magnitude and time of the biggest break ## WARNING: this is requires some time :-).
On Tue, Mar 20, 2012 at 5:35 PM, White, Joanne <Joanne.White at nrcan-rncan.gc.ca> wrote:
> 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-co
>> m
>> mits
More information about the Bfast-commits
mailing list