Dear Robert.

Thank you very much.

The code run perfectly well for a rasterstack that has a dimension of 1200
rows x 1300 columns and 600 layers. It only took about an hour to process.

# fun working fine with real data
fun <- function(x) {
if (all(is.na(x))) {
return(cbind(NA,NA))
}
m <- lm(x ~ time)
s  <- summary(m)
r2 <- s$r.squared
p <- s$coefficients[2,4]
cbind(r2, p)
}


Just wondering if I would run the lm() for a certain time period e.g from
layer 100 to 300, is it correct to use subset as below? The code seems
running but I am not sure wheter it subsetted the rows or layers:

fun <- function(x) {
if (all(is.na(x))) {
return(cbind(NA,NA))
}
m <- lm(x ~ time, subset=100:300)
s  <- summary(m)
r2 <- s$r.squared
p <- s$coefficients[2,4]
cbind(r2, p)
}

Thanks in advance.




On Sun, Nov 24, 2013 at 3:19 AM, Robert J. Hijmans <r.hijmans@gmail.com>wrote:

> Eddie, here is how you can do that (presented as a general worked
> example for writing functions for calc) Robert
>
> #1 create test data
> library(raster)
> r <- raster(nrow=10, ncol=10)
> set.seed(0)
> s <- stack(lapply(1:12, function(i) setValues(r, rnorm(ncell(r), i, 3) )))
> time <- 1:nlayers(s)
>
> #2 create an appropriate function that returns the values of interest
>
> fun <- function(x) {
>   m <- lm(x ~ time)
>   s  <- summary(m)
>   r2 <- s$r.squared
>   p <- s$coefficients[2,4]
>   cbind(r2, p)
> }
>
>
> #3 test the function for a single cell (vector)
> test <- as.vector(s[1])
> fun(test)
> #            r2           p
> #[1,] 0.453078 0.01643562
> # OK
>
> #4 also test with some NA values
> test[3] <- NA
> fun(test)
> #            r2          p
> #[1,] 0.4025646 0.03600484
> #OK
>
> #5 also test with only NA values
> test[] <- NA
> fun(test)
> #Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
> #  0 (non-NA) cases
>
> #6 not good, fix function
> fun <- function(x) {
> if (all(is.na(x))) {
> return(cbind(NA,NA))
> }
> m <- lm(x ~ time)
> s  <- summary(m)
> r2 <- s$r.squared
> p <- s$coefficients[2,4]
> cbind(r2, p)
> }
>
> #7 test again
> test <- as.vector(s[1])
> fun(test)
> #            r2           p
> #[1,] 0.453078 0.01643562
> test[3] <- NA
> fun(test)
> #            r2          p
> #[1,] 0.4025646 0.03600484
> test[] <- NA
> fun(test)
> #     [,1] [,2]
> #[1,]   NA   NA
> # all good now
>
> #8 test within apply (although this may not be strictly required)
>
> t(apply(s[1:5], 1, fun))
> #          [,1]         [,2]
> #[1,] 0.4530780 0.0164356218
> #[2,] 0.6958559 0.0007421297
> #[3,] 0.6469321 0.0016099809
> #[4,] 0.3184170 0.0559753606
> #[5,] 0.4496749 0.0170002701
> # looks good
>
> #9 test with Raster* data
>  s[1:5] <- NA
>  r <- calc(s, fun)
>  plot(r)
>
> #10 all good, now try with the real data
>
> On Sat, Nov 23, 2013 at 8:12 AM, Eddie Smith <eddieatr@gmail.com> wrote:
> > Dear list,
> >
> > 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)
> >
> > stest <- stack(s1, s2)
> > fun <- function(x) { lm(x ~ time)$coefficients }
> > x <- calc(stest, fun)
> >
> > >From the raster document, x contains the intercept and slope of the
> model.
> >
> > Anybody could help me in producing other details of the model e.g
> standard
> > error, r-squared, p-value etc?
> >
> > Thanks in advance.
> >
> >         [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]

