There are some difference in the behavior of 'calc' between functions,
because of attempts to accommodate different functions & intentions.
But in 'raster' 1.7-4 (available from R-Forge in ~ 24 hrs; and from
CRAN soon), the below will work:
library(raster)
#creating some data
r <- raster(nrow=10, ncol=10)
s1 <- s2<- list()
for (i in 1:12) {
s1[i] <- setValues(r, rnorm(ncell(r), i, 3) )
s2[i] <- setValues(r, rnorm(ncell(r), i, 3) )
}
s1 <- stack(s1)
s2 <- stack(s2)
# regression of values in one brick (or stack) with another (Jacob's suggestion)
s <- stack(s1, s2)
# s1 and s2 have 12 layers
fun <- function(x) { lm(x[1:12] ~ x[13:24])$coefficients[2] }
x1 <- calc(s, fun)
# regression of values in one brick (or stack) with 'time'
time <- 1:nlayers(s)
fun <- function(x) { lm(x ~ time)$coefficients[2] }
x2 <- calc(s, fun)
# get multiple layers, e.g. the slope _and_ intercept
fun <- function(x) { lm(x ~ time)$coefficients }
x3 <- calc(s, fun)
If this does not work in your version, you can use apply( ) as in what
I send earlier.
Robert
On Fri, Nov 26, 2010 at 2:56 PM, Robert J. Hijmans <r.hijm...@gmail.com> wrote:
> It seems that 'calc' does not like this (any more; another thing to
> fix) . If your rasters are not very large you can use apply, which
> makes it much faster:
>
> library(raster)
> #creating some data
> r <- raster(nrow=10, ncol=10)
> s <- list()
> for (i in 1:25) { s[i] <- setValues(r, rnorm(ncell(r), i, 3) ) }
> s <- stack(s)
>
> # regression
> time <- 1:nlayers(s)
> fun <- function(x) { lm(x ~ time)$coefficients[2] }
> r[] <- apply(getValues(s), 1, fun)
>
> Robert
>
>
>
> On Fri, Nov 26, 2010 at 2:51 PM, Jacob van Etten
> <jacobvanet...@yahoo.com> wrote:
>> you could try this approach (use calc whenever you can):
>>
>> (supposing your bricks have 12 layers)
>>
>> br3 <- stack(brick1, brick2)
>> lmS <- function(x) lm(x[1:12] ~ x[13:24)$coefficients[2]
>> r <- calc(br3, lmS)
>>
>> Jacob.
>>
>> --- On Fri, 26/11/10, steven mosher <mosherste...@gmail.com> wrote:
>>
>> From: steven mosher <mosherste...@gmail.com>
>> Subject: Re: [R-sig-Geo] gridded time series analysis
>> To: "Martin" <martin_bra...@gmx.net>
>> Cc: r-sig-geo@stat.math.ethz.ch
>> Date: Friday, 26 November, 2010, 23:33
>>
>> that's cool, I'm also interested in a similar problem but just with one
>> brick ending up with a slope raster as the output. It may be possible with
>> stackApply(). have a look. or maybe robert will chime in
>>
>>
>>
>> On Fri, Nov 26, 2010 at 1:35 PM, Martin <martin_bra...@gmx.net> wrote:
>>
>>>
>>> this is what I did to perform a regression between two bricks (each brick
>>> represent a time series):
>>>
>>> r <- raster(brick1)
>>> for (i in 1:ncell(r)) {
>>> r[i] = lm(as.ts(cellValues(brick1, i)) ~ as.ts(cellValues(brick2,
>>> i)))$coefficients[2]
>>> }
>>>
>>> The result will be a slope raster, but it really takes a lot of time, so
>>> maybe there is a better solution..
>>>
>>>
>>> --
>>> View this message in context:
>>> http://r-sig-geo.2731867.n2.nabble.com/gridded-time-series-analysis-tp5775651p5778472.html
>>> Sent from the R-sig-geo mailing list archive at Nabble.com.
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo@stat.math.ethz.ch
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo@stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo@stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
