Re: [R-sig-Geo] 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 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
> 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 wrote:
>>
>> From: steven mosher
>> Subject: Re: [R-sig-Geo] gridded time series analysis
>> To: "Martin"
>> Cc: [email protected]
>> 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 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
>>> [email protected]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>> [[alternative HTML version deleted]]
>>
>> ___
>> R-sig-Geo mailing list
>> [email protected]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>>
>> ___
>> R-sig-Geo mailing list
>> [email protected]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>
___
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] gridded time series analysis
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
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 wrote:
>
> From: steven mosher
> Subject: Re: [R-sig-Geo] gridded time series analysis
> To: "Martin"
> Cc: [email protected]
> 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 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
>> [email protected]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> [email protected]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
>
> [[alternative HTML version deleted]]
>
>
> ___
> R-sig-Geo mailing list
> [email protected]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
___
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] gridded time series analysis
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 wrote: From: steven mosher Subject: Re: [R-sig-Geo] gridded time series analysis To: "Martin" Cc: [email protected] 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 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 > [email protected] > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > [[alternative HTML version deleted]] ___ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] gridded time series analysis
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 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
> [email protected]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] gridded time series analysis
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
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
