gridded time series analysis
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.hijmans at 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
<jacobvanetten at 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 <moshersteven at gmail.com> wrote: From: steven mosher <moshersteven at gmail.com> Subject: Re: [R-sig-Geo] gridded time series analysis To: "Martin" <martin_brandt at gmx.net> Cc: r-sig-geo at 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_brandt at 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 at 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 at 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 at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo