gridded time series analysis

classic Classic list List threaded Threaded
42 messages Options
123
Reply | Threaded
Open this post in threaded view
|

gridded time series analysis

Advait Godbole
Dear Members,

I have two variables on a 507 x 356 grid and each grid cell has a time
series of length 120.  Is it possible to perform a time series analysis on
the entire grid all at the same time, with the ultimate goal of doing linear
regressions for each grid point for the 2 variables? Are there functions
available or will I have to utilize for loops and repeat it for each grid
cell? I have been trying out the package "ts" but not sure if it can handle
3d variables. Any pointers will be greatly appreciated.

Thanks,

--
advait godbole

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

Martin
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..

Martin Brandt

Department of Geography and
Regional Research
UZA II, Althanstr. 14
1090 Wien, AUSTRIA
http://geooekologie.univie.ac.at
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

steven mosher
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 <[hidden email]> 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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

Jacob van Etten
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 <[hidden email]> wrote:

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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo



     
        [[alternative HTML version deleted]]


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

Advait Godbole
Thanks everyone for these leads..However, I am running into problems with
data handling. I used ncdf to extract the two variables (507 x 356 x 120 -
lon,lat,time) and converted them into array objects. The suggestions
provided all utilize the raster package, but using "brick" on the file
returns this error: "rcmraster <- brick("regcm.cumul.hdd.10k.96_05.new.nc")
Error in if (e@xmin > -400 & e@xmax < 400 & e@ymin > -90.1 & e@ymax <  :
missing value where TRUE/FALSE needed"

I have a fair number of NAs (all missing values in the original variables).
Is there a way to set a flag for that? Or any other workaround?

Thanks much,

advait

On Fri, Nov 26, 2010 at 6:25 PM, Robert J. Hijmans <[hidden email]>wrote:

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



--
advait godbole

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

Robert Hijmans
Advait,

That suggests that you have missing values in the lon or lat
dimensions. Is that true? That would be awkward, but perhaps something
else is going on. Can you send me your file? If not can you send the
results of this:

filename <- "regcm.cumul.hdd.10k.96_05.new.nc"
nc <- open.ncdf(filename)
length(nc$dim[[1]]$vals)
length(nc$dim[[2]]$vals)
nc$dim[[1]]$vals
nc$dim[[2]]$vals
close.ncdf(nc)

Thanks, Robert

On Sat, Nov 27, 2010 at 8:36 AM, Advait Godbole <[hidden email]> wrote:

> Thanks everyone for these leads..However, I am running into problems with
> data handling. I used ncdf to extract the two variables (507 x 356 x 120 -
> lon,lat,time) and converted them into array objects. The suggestions
> provided all utilize the raster package, but using "brick" on the file
> returns this error: "rcmraster <- brick("regcm.cumul.hdd.10k.96_05.new.nc")
> Error in if (e@xmin > -400 & e@xmax < 400 & e@ymin > -90.1 & e@ymax <  :
> missing value where TRUE/FALSE needed"
>
> I have a fair number of NAs (all missing values in the original variables).
> Is there a way to set a flag for that? Or any other workaround?
> Thanks much,
> advait
> On Fri, Nov 26, 2010 at 6:25 PM, Robert J. Hijmans <[hidden email]>
> wrote:
>>
>> 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 <[hidden email]>
>> 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
>> > <[hidden email]> 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 <[hidden email]> wrote:
>> >>
>> >> From: steven mosher <[hidden email]>
>> >> Subject: Re: [R-sig-Geo] gridded time series analysis
>> >> To: "Martin" <[hidden email]>
>> >> Cc: [hidden email]
>> >> 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 <[hidden email]> 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
>> >>> [hidden email]
>> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >>>
>> >>
>> >>     [[alternative HTML version deleted]]
>> >>
>> >> _______________________________________________
>> >> R-sig-Geo mailing list
>> >> [hidden email]
>> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >>
>> >>
>> >>
>> >>
>> >>        [[alternative HTML version deleted]]
>> >>
>> >>
>> >> _______________________________________________
>> >> R-sig-Geo mailing list
>> >> [hidden email]
>> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >>
>> >>
>> >
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
> --
> advait godbole
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

Martin
In reply to this post by Robert Hijmans
thanks Robert, this works fine.
but could it be that the functions have a problem with NA-values?
your Regression functions return this error:


> x3 <- calc(tt, fun)
error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
  0 (non-NA) cases

when i use rasters with NA values...it works without problems with rasters without NAs..

cheers,
Martin



Martin Brandt

Department of Geography and
Regional Research
UZA II, Althanstr. 14
1090 Wien, AUSTRIA
http://geooekologie.univie.ac.at
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

Robert Hijmans
Martin,

I think the problem is whith the lm function that returns an error
rather than NA:

> a <- rep(NA,10)
> lm(a~a)
Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
  0 (non-NA) cases

#rather than something like this:
> fun=function(x) { lm(x[1]~x[2])$coefficients[2] }
> fun(c(NA,NA))
Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
  0 (non-NA) cases

#you could catch this in your function and do something like:

> fun=function(x) { if (is.na(x[1])){ NA } else { lm(x[1]~x[2])$coefficients[2] }}
> fun(c(NA,NA))
[1] NA


and then:

calc(s, fun=fun)

lm and relatives have an na.action argument but I do not think that
can come to the rescue here.

Roberrt

On Sun, Nov 28, 2010 at 4:33 AM, Martin <[hidden email]> wrote:

>
> thanks Robert, this works fine.
> but could it be that the functions have a problem with NA-values?
> your Regression functions return this error:
>
>
>> x3 <- calc(tt, fun)
> error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
>  0 (non-NA) cases
>
> when i use rasters with NA values...it works without problems with rasters
> without NAs..
>
> cheers,
> Martin
>
>
>
>
> --
> View this message in context: http://r-sig-geo.2731867.n2.nabble.com/gridded-time-series-analysis-tp5775651p5781554.html
> Sent from the R-sig-geo mailing list archive at Nabble.com.
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

Martin
that works Robert, thanks a lot!
Martin Brandt

Department of Geography and
Regional Research
UZA II, Althanstr. 14
1090 Wien, AUSTRIA
http://geooekologie.univie.ac.at
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

Martin
another thing which might be interesting:

I know I can call the intercept with "coefficients[1]" and the slope with "coefficients[2]". But is it also possible to get a gridded result with the other lm regression results (normaly available via "summary()") , like R², p-value, or residuals?

Martin

Martin Brandt

Department of Geography and
Regional Research
UZA II, Althanstr. 14
1090 Wien, AUSTRIA
http://geooekologie.univie.ac.at
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

Advait Godbole
In reply to this post by Robert Hijmans
sorry to disrupt the flow of the conversation, but after doing what Robert
suggested, this is what was returned:

> length(frcm$dim[[1]]$vals)
[1] 21659040

> length(frcm$dim[[2]]$vals)
[1] 21659040

the other two functions return the latitude and longitude arrays not
populated with NAs/missing values. I havent reproduced them here because of
their number.

Regards,

advait

On Sat, Nov 27, 2010 at 3:48 PM, Robert J. Hijmans <[hidden email]>wrote:

> Advait,
>
> That suggests that you have missing values in the lon or lat
> dimensions. Is that true? That would be awkward, but perhaps something
> else is going on. Can you send me your file? If not can you send the
> results of this:
>
> filename <- "regcm.cumul.hdd.10k.96_05.new.nc"
> nc <- open.ncdf(filename)
> length(nc$dim[[1]]$vals)
> length(nc$dim[[2]]$vals)
> nc$dim[[1]]$vals
> nc$dim[[2]]$vals
> close.ncdf(nc)
>
> Thanks, Robert
>
> On Sat, Nov 27, 2010 at 8:36 AM, Advait Godbole <[hidden email]>
> wrote:
> > Thanks everyone for these leads..However, I am running into problems with
> > data handling. I used ncdf to extract the two variables (507 x 356 x 120
> -
> > lon,lat,time) and converted them into array objects. The suggestions
> > provided all utilize the raster package, but using "brick" on the file
> > returns this error: "rcmraster <- brick("
> regcm.cumul.hdd.10k.96_05.new.nc")
> > Error in if (e@xmin > -400 & e@xmax < 400 & e@ymin > -90.1 & e@ymax <  :
> > missing value where TRUE/FALSE needed"
> >
> > I have a fair number of NAs (all missing values in the original
> variables).
> > Is there a way to set a flag for that? Or any other workaround?
> > Thanks much,
> > advait
> > On Fri, Nov 26, 2010 at 6:25 PM, Robert J. Hijmans <[hidden email]>
> > wrote:
> >>
> >> 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 <[hidden email]
> >
> >> 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
> >> > <[hidden email]> 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 <[hidden email]> wrote:
> >> >>
> >> >> From: steven mosher <[hidden email]>
> >> >> Subject: Re: [R-sig-Geo] gridded time series analysis
> >> >> To: "Martin" <[hidden email]>
> >> >> Cc: [hidden email]
> >> >> 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 <[hidden email]>
> 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
> >> >>> [hidden email]
> >> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >> >>>
> >> >>
> >> >>     [[alternative HTML version deleted]]
> >> >>
> >> >> _______________________________________________
> >> >> R-sig-Geo mailing list
> >> >> [hidden email]
> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >> >>
> >> >>
> >> >>
> >> >>
> >> >>        [[alternative HTML version deleted]]
> >> >>
> >> >>
> >> >> _______________________________________________
> >> >> R-sig-Geo mailing list
> >> >> [hidden email]
> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >> >>
> >> >>
> >> >
> >>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> [hidden email]
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> >
> >
> > --
> > advait godbole
> >
>



--
advait godbole

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

Robert Hijmans
Advait,
I do not think I can help you further unless you send me your file;
but this suggests that it is either not following the CF standards, or
something is wrong with the file, or that there is something I have
not understood about these files.
You said you have 507 x 356 x 120 values. I would expect
length(frcm$dim[[1]]$vals) == 507 and length(frcm$dim[[2]]$vals) ==
356
Robert

On Mon, Nov 29, 2010 at 6:11 AM, Advait Godbole <[hidden email]> wrote:

> sorry to disrupt the flow of the conversation, but after doing what Robert
> suggested, this is what was returned:
>> length(frcm$dim[[1]]$vals)
> [1] 21659040
>> length(frcm$dim[[2]]$vals)
> [1] 21659040
> the other two functions return the latitude and longitude arrays not
> populated with NAs/missing values. I havent reproduced them here because of
> their number.
> Regards,
> advait
> On Sat, Nov 27, 2010 at 3:48 PM, Robert J. Hijmans <[hidden email]>
> wrote:
>>
>> Advait,
>>
>> That suggests that you have missing values in the lon or lat
>> dimensions. Is that true? That would be awkward, but perhaps something
>> else is going on. Can you send me your file? If not can you send the
>> results of this:
>>
>> filename <- "regcm.cumul.hdd.10k.96_05.new.nc"
>> nc <- open.ncdf(filename)
>> length(nc$dim[[1]]$vals)
>> length(nc$dim[[2]]$vals)
>> nc$dim[[1]]$vals
>> nc$dim[[2]]$vals
>> close.ncdf(nc)
>>
>> Thanks, Robert
>>
>> On Sat, Nov 27, 2010 at 8:36 AM, Advait Godbole <[hidden email]>
>> wrote:
>> > Thanks everyone for these leads..However, I am running into problems
>> > with
>> > data handling. I used ncdf to extract the two variables (507 x 356 x 120
>> > -
>> > lon,lat,time) and converted them into array objects. The suggestions
>> > provided all utilize the raster package, but using "brick" on the file
>> > returns this error: "rcmraster <-
>> > brick("regcm.cumul.hdd.10k.96_05.new.nc")
>> > Error in if (e@xmin > -400 & e@xmax < 400 & e@ymin > -90.1 & e@ymax <  :
>> > missing value where TRUE/FALSE needed"
>> >
>> > I have a fair number of NAs (all missing values in the original
>> > variables).
>> > Is there a way to set a flag for that? Or any other workaround?
>> > Thanks much,
>> > advait
>> > On Fri, Nov 26, 2010 at 6:25 PM, Robert J. Hijmans <[hidden email]>
>> > wrote:
>> >>
>> >> 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
>> >> <[hidden email]>
>> >> 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
>> >> > <[hidden email]> 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 <[hidden email]> wrote:
>> >> >>
>> >> >> From: steven mosher <[hidden email]>
>> >> >> Subject: Re: [R-sig-Geo] gridded time series analysis
>> >> >> To: "Martin" <[hidden email]>
>> >> >> Cc: [hidden email]
>> >> >> 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 <[hidden email]>
>> >> >> 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
>> >> >>> [hidden email]
>> >> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >> >>>
>> >> >>
>> >> >>     [[alternative HTML version deleted]]
>> >> >>
>> >> >> _______________________________________________
>> >> >> R-sig-Geo mailing list
>> >> >> [hidden email]
>> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >> >>
>> >> >>
>> >> >>
>> >> >>
>> >> >>        [[alternative HTML version deleted]]
>> >> >>
>> >> >>
>> >> >> _______________________________________________
>> >> >> R-sig-Geo mailing list
>> >> >> [hidden email]
>> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >> >>
>> >> >>
>> >> >
>> >>
>> >> _______________________________________________
>> >> R-sig-Geo mailing list
>> >> [hidden email]
>> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >
>> >
>> >
>> > --
>> > advait godbole
>> >
>
>
>
> --
> advait godbole
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

Robert Hijmans
In reply to this post by Martin
You can return anything you like, as long as it is numeric (or a
vector of numerics). R^2 should come out like this:

 fun=function(x) { if (is.na(x[1])){ NA } else { m <-
lm(x[1]~x[2])$coefficients[2]; summary(m)$r.squared  }}


slope and R^2

 fun=function(x) { if (is.na(x[1])){ NA } else { m <-
lm(x[1]~x[2])$coefficients[2]; c(m$coefficients[2],
summary(m)$r.squared  }}

Robert

On Mon, Nov 29, 2010 at 1:58 AM, Martin <[hidden email]> wrote:

>
> another thing which might be interesting:
>
> I know I can call the intercept with "coefficients[1]" and the slope with
> "coefficients[2]". But is it also possible to get a gridded result with the
> other lm regression results (normaly available via "summary()") , like R²,
> p-value, or residuals?
>
> Martin
>
>
> --
> View this message in context: http://r-sig-geo.2731867.n2.nabble.com/gridded-time-series-analysis-tp5775651p5783798.html
> Sent from the R-sig-geo mailing list archive at Nabble.com.
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

Advait Godbole
In reply to this post by Robert Hijmans
Dear Robert,

here is the link to the file:
http://dl.dropbox.com/u/4030944/regcm.cumul.hdd.10k.96_05.new.nc. It is a
large file ~400 MB. Let me know if dropbox is not convenient and you would
rather do ftp -  but our machines are down for maintenance and I may not be
immediately able to get access.

I was going over the NCL code that processes the files but did not find
anything glaringly amiss, so I will have to look at it in greater detail
since it is clear that the number of entries in the variables is not right
as you've pointed out. There is a chance that I have not written out the
co-ordinates properly.

Thanks,

Advait

On Mon, Nov 29, 2010 at 1:14 PM, Robert J. Hijmans <[hidden email]>wrote:

> Advait,
> I do not think I can help you further unless you send me your file;
> but this suggests that it is either not following the CF standards, or
> something is wrong with the file, or that there is something I have
> not understood about these files.
> You said you have 507 x 356 x 120 values. I would expect
> length(frcm$dim[[1]]$vals) == 507 and length(frcm$dim[[2]]$vals) ==
> 356
> Robert
>
> On Mon, Nov 29, 2010 at 6:11 AM, Advait Godbole <[hidden email]>
> wrote:
> > sorry to disrupt the flow of the conversation, but after doing what
> Robert
> > suggested, this is what was returned:
> >> length(frcm$dim[[1]]$vals)
> > [1] 21659040
> >> length(frcm$dim[[2]]$vals)
> > [1] 21659040
> > the other two functions return the latitude and longitude arrays not
> > populated with NAs/missing values. I havent reproduced them here because
> of
> > their number.
> > Regards,
> > advait
> > On Sat, Nov 27, 2010 at 3:48 PM, Robert J. Hijmans <[hidden email]>
> > wrote:
> >>
> >> Advait,
> >>
> >> That suggests that you have missing values in the lon or lat
> >> dimensions. Is that true? That would be awkward, but perhaps something
> >> else is going on. Can you send me your file? If not can you send the
> >> results of this:
> >>
> >> filename <- "regcm.cumul.hdd.10k.96_05.new.nc"
> >> nc <- open.ncdf(filename)
> >> length(nc$dim[[1]]$vals)
> >> length(nc$dim[[2]]$vals)
> >> nc$dim[[1]]$vals
> >> nc$dim[[2]]$vals
> >> close.ncdf(nc)
> >>
> >> Thanks, Robert
> >>
> >> On Sat, Nov 27, 2010 at 8:36 AM, Advait Godbole <
> [hidden email]>
> >> wrote:
> >> > Thanks everyone for these leads..However, I am running into problems
> >> > with
> >> > data handling. I used ncdf to extract the two variables (507 x 356 x
> 120
> >> > -
> >> > lon,lat,time) and converted them into array objects. The suggestions
> >> > provided all utilize the raster package, but using "brick" on the file
> >> > returns this error: "rcmraster <-
> >> > brick("regcm.cumul.hdd.10k.96_05.new.nc")
> >> > Error in if (e@xmin > -400 & e@xmax < 400 & e@ymin > -90.1 & e@ymax <
>  :
> >> > missing value where TRUE/FALSE needed"
> >> >
> >> > I have a fair number of NAs (all missing values in the original
> >> > variables).
> >> > Is there a way to set a flag for that? Or any other workaround?
> >> > Thanks much,
> >> > advait
> >> > On Fri, Nov 26, 2010 at 6:25 PM, Robert J. Hijmans <
> [hidden email]>
> >> > wrote:
> >> >>
> >> >> 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
> >> >> <[hidden email]>
> >> >> 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
> >> >> > <[hidden email]> 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 <[hidden email]>
> wrote:
> >> >> >>
> >> >> >> From: steven mosher <[hidden email]>
> >> >> >> Subject: Re: [R-sig-Geo] gridded time series analysis
> >> >> >> To: "Martin" <[hidden email]>
> >> >> >> Cc: [hidden email]
> >> >> >> 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 <[hidden email]>
> >> >> >> 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
> >> >> >>> [hidden email]
> >> >> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >> >> >>>
> >> >> >>
> >> >> >>     [[alternative HTML version deleted]]
> >> >> >>
> >> >> >> _______________________________________________
> >> >> >> R-sig-Geo mailing list
> >> >> >> [hidden email]
> >> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >> >> >>
> >> >> >>
> >> >> >>
> >> >> >>
> >> >> >>        [[alternative HTML version deleted]]
> >> >> >>
> >> >> >>
> >> >> >> _______________________________________________
> >> >> >> R-sig-Geo mailing list
> >> >> >> [hidden email]
> >> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >> >> >>
> >> >> >>
> >> >> >
> >> >>
> >> >> _______________________________________________
> >> >> R-sig-Geo mailing list
> >> >> [hidden email]
> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >> >
> >> >
> >> >
> >> > --
> >> > advait godbole
> >> >
> >
> >
> >
> > --
> > advait godbole
> >
>



--
advait godbole

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

Robert Hijmans
Advait,

The is an error in the file, I think, see below:

> library(ncdf)
> filename = "c:/downloads/regcm.cumul.hdd.10k.96_05.new.nc"
> nc <- open.ncdf(filename)
>
> zvar <- 'hdd'
> nc$var[[zvar]]$dim[[1]]$len
[1] 507
> nc$var[[zvar]]$dim[[2]]$len
[1] 356
>
> xx <- nc$var[[zvar]]$dim[[1]]$vals
> # should be a vector of 507 longitudes but:
> dim(xx)
[1] 507 356 120
>
> # works anyway
> xrange <- c(min(xx), max(xx))
> xrange
[1] -137.25705  -62.18677
>
> yy <- nc$var[[zvar]]$dim[[2]]$vals
> # should be a vector of 356 latitudes but:
> dim(yy)
[1] 507 356 120
>
> # this fails
> yrange <- c(min(yy), max(yy))
> yrange
[1] NaN NaN
> # and that is the only reason why raster failed.
>
> #OK
> yy[,,1][1:10,350:355]
          [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
 [1,] 51.49148 51.57122 51.65092 51.73055 51.81014 51.88967
 [2,] 51.52888 51.60866 51.68839 51.76807 51.84769 51.92727
 [3,] 51.56616 51.64598 51.72575 51.80547 51.88513 51.96474
 [4,] 51.60332 51.68318 51.76299 51.84275 51.92245 52.00210
 [5,] 51.64036 51.72026 51.80011 51.87991 51.95965 52.03934
 [6,] 51.67728 51.75722 51.83711 51.91695 51.99673 52.07646
 [7,] 51.71409 51.79407 51.87399 51.95387 52.03369 52.11345
 [8,] 51.75077 51.83079 51.91076 51.99067 52.07053 52.15033
 [9,] 51.78734 51.86739 51.94740 52.02735 52.10725 52.18709
[10,] 51.82378 51.90388 51.98392 52.06391 52.14385 52.22373
>
> #bad values here (and note the NaN):
> yy[,,120][1:10,350:355]
                [,1]           [,2]           [,3]           [,4]
     [,5]           [,6]
 [1,] -1.313847e-229 -4.895880e-272 -7.553765e-266   1.293814e+58
9.551441e-293  -3.291580e-36
 [2,] -5.321424e+272 -3.076997e-223  -1.515780e-33 -9.895836e+171
3.115686e-149  -4.970539e+92
 [3,]  5.259277e-203   5.980335e-56  2.893449e-212  2.974200e+274
3.186547e-218  2.503689e+296
 [4,] -3.977129e+225 -3.696719e-283   1.284612e+42   1.887923e+28
-2.322520e+268  5.005820e-108
 [5,]  -6.569744e+88   1.824506e-74   7.663356e+39  1.608759e-201
-1.433761e+110 -8.804183e-211
 [6,] -4.237822e-162            NaN  7.405904e-275  5.859877e-258
1.339820e+110  1.631993e+229
 [7,]   5.603325e-64   4.188764e+37 -2.188887e+243 -4.657907e+280
1.795268e-185  1.166057e-175
 [8,] -2.112804e-229  4.536625e+137  2.910322e-278  4.407645e-274
8.886205e-239 -1.913860e+238
 [9,]  -9.126350e+44 -1.866422e+177 -2.578975e+256  3.275379e+208
-2.657457e+75  -1.258259e-90
[10,]  6.698325e+216 -2.958355e+137 -1.058748e-178 -1.776776e+215
-6.368807e+78 -5.337653e+299
>
> # this would have worked
> yrange <- c(min(yy[,,1]), max(yy[,,1]))
> yrange
[1] 22.17641 57.31006
>
> close.ncdf(nc)


Best, Robert


On Tue, Nov 30, 2010 at 10:37 AM, Advait Godbole
<[hidden email]> wrote:

> Dear Robert,
> here is the link to the file:
> http://dl.dropbox.com/u/4030944/regcm.cumul.hdd.10k.96_05.new.nc. It is a
> large file ~400 MB. Let me know if dropbox is not convenient and you would
> rather do ftp -  but our machines are down for maintenance and I may not be
> immediately able to get access.
> I was going over the NCL code that processes the files but did not find
> anything glaringly amiss, so I will have to look at it in greater detail
> since it is clear that the number of entries in the variables is not right
> as you've pointed out. There is a chance that I have not written out the
> co-ordinates properly.
> Thanks,
> Advait
> On Mon, Nov 29, 2010 at 1:14 PM, Robert J. Hijmans <[hidden email]>
> wrote:
>>
>> Advait,
>> I do not think I can help you further unless you send me your file;
>> but this suggests that it is either not following the CF standards, or
>> something is wrong with the file, or that there is something I have
>> not understood about these files.
>> You said you have 507 x 356 x 120 values. I would expect
>> length(frcm$dim[[1]]$vals) == 507 and length(frcm$dim[[2]]$vals) ==
>> 356
>> Robert
>>
>> On Mon, Nov 29, 2010 at 6:11 AM, Advait Godbole <[hidden email]>
>> wrote:
>> > sorry to disrupt the flow of the conversation, but after doing what
>> > Robert
>> > suggested, this is what was returned:
>> >> length(frcm$dim[[1]]$vals)
>> > [1] 21659040
>> >> length(frcm$dim[[2]]$vals)
>> > [1] 21659040
>> > the other two functions return the latitude and longitude arrays not
>> > populated with NAs/missing values. I havent reproduced them here because
>> > of
>> > their number.
>> > Regards,
>> > advait
>> > On Sat, Nov 27, 2010 at 3:48 PM, Robert J. Hijmans <[hidden email]>
>> > wrote:
>> >>
>> >> Advait,
>> >>
>> >> That suggests that you have missing values in the lon or lat
>> >> dimensions. Is that true? That would be awkward, but perhaps something
>> >> else is going on. Can you send me your file? If not can you send the
>> >> results of this:
>> >>
>> >> filename <- "regcm.cumul.hdd.10k.96_05.new.nc"
>> >> nc <- open.ncdf(filename)
>> >> length(nc$dim[[1]]$vals)
>> >> length(nc$dim[[2]]$vals)
>> >> nc$dim[[1]]$vals
>> >> nc$dim[[2]]$vals
>> >> close.ncdf(nc)
>> >>
>> >> Thanks, Robert
>> >>
>> >> On Sat, Nov 27, 2010 at 8:36 AM, Advait Godbole
>> >> <[hidden email]>
>> >> wrote:
>> >> > Thanks everyone for these leads..However, I am running into problems
>> >> > with
>> >> > data handling. I used ncdf to extract the two variables (507 x 356 x
>> >> > 120
>> >> > -
>> >> > lon,lat,time) and converted them into array objects. The suggestions
>> >> > provided all utilize the raster package, but using "brick" on the
>> >> > file
>> >> > returns this error: "rcmraster <-
>> >> > brick("regcm.cumul.hdd.10k.96_05.new.nc")
>> >> > Error in if (e@xmin > -400 & e@xmax < 400 & e@ymin > -90.1 & e@ymax <
>> >> >  :
>> >> > missing value where TRUE/FALSE needed"
>> >> >
>> >> > I have a fair number of NAs (all missing values in the original
>> >> > variables).
>> >> > Is there a way to set a flag for that? Or any other workaround?
>> >> > Thanks much,
>> >> > advait
>> >> > On Fri, Nov 26, 2010 at 6:25 PM, Robert J. Hijmans
>> >> > <[hidden email]>
>> >> > wrote:
>> >> >>
>> >> >> 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
>> >> >> <[hidden email]>
>> >> >> 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
>> >> >> > <[hidden email]> 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 <[hidden email]>
>> >> >> >> wrote:
>> >> >> >>
>> >> >> >> From: steven mosher <[hidden email]>
>> >> >> >> Subject: Re: [R-sig-Geo] gridded time series analysis
>> >> >> >> To: "Martin" <[hidden email]>
>> >> >> >> Cc: [hidden email]
>> >> >> >> 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 <[hidden email]>
>> >> >> >> 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
>> >> >> >>> [hidden email]
>> >> >> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >> >> >>>
>> >> >> >>
>> >> >> >>     [[alternative HTML version deleted]]
>> >> >> >>
>> >> >> >> _______________________________________________
>> >> >> >> R-sig-Geo mailing list
>> >> >> >> [hidden email]
>> >> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >> >> >>
>> >> >> >>
>> >> >> >>
>> >> >> >>
>> >> >> >>        [[alternative HTML version deleted]]
>> >> >> >>
>> >> >> >>
>> >> >> >> _______________________________________________
>> >> >> >> R-sig-Geo mailing list
>> >> >> >> [hidden email]
>> >> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >> >> >>
>> >> >> >>
>> >> >> >
>> >> >>
>> >> >> _______________________________________________
>> >> >> R-sig-Geo mailing list
>> >> >> [hidden email]
>> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >> >
>> >> >
>> >> >
>> >> > --
>> >> > advait godbole
>> >> >
>> >
>> >
>> >
>> > --
>> > advait godbole
>> >
>
>
>
> --
> advait godbole
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

Advait Godbole
That worked perfectly. Thanks for pointing this out. I changed the way my
variables were being written to the netCDF and now "brick" works. However, I
am getting an error while using the stack function, like so:

"Error in function (classes, fdef, mtable)  :
  unable to find an inherited method for function "trim", for signature
"integer"

The two rasters show the following:
> rcmraster
class       : RasterBrick
filename    : regcm.cumul.hdd.10k.96_05.fixed.latlonvars.nc
nlayers     : 120
nrow        : 356
ncol        : 507
ncell       : 180492
projection  : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
min value   : ? ? ? ? ? ? ? ? ? ?
max value   : ? ? ? ? ? ? ? ? ? ?
xmin        : -137.3312
xmax        : -62.11259
ymin        : 22.12693
ymax        : 57.35954
xres        : 0.1483602
yres        : 0.09896801

> emitraster
class       : RasterBrick
filename    :
flip.percap.vulcan.10k.dp.build.381.RES.EIAmy.96_05.this.latest.nc
nlayers     : 120
nrow        : 356
ncol        : 507
ncell       : 180492
projection  : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
min value   : ? ? ? ? ? ? ? ? ? ?
max value   : ? ? ? ? ? ? ? ? ? ?
xmin        : -137.3312
xmax        : -62.11259
ymin        : 22.12693
ymax        : 57.35954
xres        : 0.1483602
yres        : 0.09896801

A getValues call shows a lot of NAs. Now, my original file does have a large
number of missing values so that is why this should happen, but I suspect
something is going wrong since the min and max values are not being
displayed for either of these rasters. Is that the source of the problem in
the "stack" function?

Thanks,

Advait

On Tue, Nov 30, 2010 at 5:18 PM, Robert J. Hijmans <[hidden email]>wrote:

> Advait,
>
> The is an error in the file, I think, see below:
>
> > library(ncdf)
> > filename = "c:/downloads/regcm.cumul.hdd.10k.96_05.new.nc"
> > nc <- open.ncdf(filename)
> >
> > zvar <- 'hdd'
> > nc$var[[zvar]]$dim[[1]]$len
> [1] 507
> > nc$var[[zvar]]$dim[[2]]$len
> [1] 356
> >
> > xx <- nc$var[[zvar]]$dim[[1]]$vals
> > # should be a vector of 507 longitudes but:
> > dim(xx)
> [1] 507 356 120
> >
> > # works anyway
> > xrange <- c(min(xx), max(xx))
> > xrange
> [1] -137.25705  -62.18677
> >
> > yy <- nc$var[[zvar]]$dim[[2]]$vals
> > # should be a vector of 356 latitudes but:
> > dim(yy)
> [1] 507 356 120
> >
> > # this fails
> > yrange <- c(min(yy), max(yy))
> > yrange
> [1] NaN NaN
> > # and that is the only reason why raster failed.
> >
> > #OK
> > yy[,,1][1:10,350:355]
>          [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
>  [1,] 51.49148 51.57122 51.65092 51.73055 51.81014 51.88967
>  [2,] 51.52888 51.60866 51.68839 51.76807 51.84769 51.92727
>  [3,] 51.56616 51.64598 51.72575 51.80547 51.88513 51.96474
>  [4,] 51.60332 51.68318 51.76299 51.84275 51.92245 52.00210
>  [5,] 51.64036 51.72026 51.80011 51.87991 51.95965 52.03934
>  [6,] 51.67728 51.75722 51.83711 51.91695 51.99673 52.07646
>  [7,] 51.71409 51.79407 51.87399 51.95387 52.03369 52.11345
>  [8,] 51.75077 51.83079 51.91076 51.99067 52.07053 52.15033
>  [9,] 51.78734 51.86739 51.94740 52.02735 52.10725 52.18709
> [10,] 51.82378 51.90388 51.98392 52.06391 52.14385 52.22373
> >
> > #bad values here (and note the NaN):
> > yy[,,120][1:10,350:355]
>                [,1]           [,2]           [,3]           [,4]
>     [,5]           [,6]
>  [1,] -1.313847e-229 -4.895880e-272 -7.553765e-266   1.293814e+58
> 9.551441e-293  -3.291580e-36
>  [2,] -5.321424e+272 -3.076997e-223  -1.515780e-33 -9.895836e+171
> 3.115686e-149  -4.970539e+92
>  [3,]  5.259277e-203   5.980335e-56  2.893449e-212  2.974200e+274
> 3.186547e-218  2.503689e+296
>  [4,] -3.977129e+225 -3.696719e-283   1.284612e+42   1.887923e+28
> -2.322520e+268  5.005820e-108
>  [5,]  -6.569744e+88   1.824506e-74   7.663356e+39  1.608759e-201
> -1.433761e+110 -8.804183e-211
>  [6,] -4.237822e-162            NaN  7.405904e-275  5.859877e-258
> 1.339820e+110  1.631993e+229
>  [7,]   5.603325e-64   4.188764e+37 -2.188887e+243 -4.657907e+280
> 1.795268e-185  1.166057e-175
>  [8,] -2.112804e-229  4.536625e+137  2.910322e-278  4.407645e-274
> 8.886205e-239 -1.913860e+238
>  [9,]  -9.126350e+44 -1.866422e+177 -2.578975e+256  3.275379e+208
> -2.657457e+75  -1.258259e-90
> [10,]  6.698325e+216 -2.958355e+137 -1.058748e-178 -1.776776e+215
> -6.368807e+78 -5.337653e+299
> >
> > # this would have worked
> > yrange <- c(min(yy[,,1]), max(yy[,,1]))
> > yrange
> [1] 22.17641 57.31006
> >
> > close.ncdf(nc)
>
>
> Best, Robert
>
>
> On Tue, Nov 30, 2010 at 10:37 AM, Advait Godbole
> <[hidden email]> wrote:
> > Dear Robert,
> > here is the link to the file:
> > http://dl.dropbox.com/u/4030944/regcm.cumul.hdd.10k.96_05.new.nc. It is
> a
> > large file ~400 MB. Let me know if dropbox is not convenient and you
> would
> > rather do ftp -  but our machines are down for maintenance and I may not
> be
> > immediately able to get access.
> > I was going over the NCL code that processes the files but did not find
> > anything glaringly amiss, so I will have to look at it in greater detail
> > since it is clear that the number of entries in the variables is not
> right
> > as you've pointed out. There is a chance that I have not written out the
> > co-ordinates properly.
> > Thanks,
> > Advait
> > On Mon, Nov 29, 2010 at 1:14 PM, Robert J. Hijmans <[hidden email]>
> > wrote:
> >>
> >> Advait,
> >> I do not think I can help you further unless you send me your file;
> >> but this suggests that it is either not following the CF standards, or
> >> something is wrong with the file, or that there is something I have
> >> not understood about these files.
> >> You said you have 507 x 356 x 120 values. I would expect
> >> length(frcm$dim[[1]]$vals) == 507 and length(frcm$dim[[2]]$vals) ==
> >> 356
> >> Robert
> >>
> >> On Mon, Nov 29, 2010 at 6:11 AM, Advait Godbole <
> [hidden email]>
> >> wrote:
> >> > sorry to disrupt the flow of the conversation, but after doing what
> >> > Robert
> >> > suggested, this is what was returned:
> >> >> length(frcm$dim[[1]]$vals)
> >> > [1] 21659040
> >> >> length(frcm$dim[[2]]$vals)
> >> > [1] 21659040
> >> > the other two functions return the latitude and longitude arrays not
> >> > populated with NAs/missing values. I havent reproduced them here
> because
> >> > of
> >> > their number.
> >> > Regards,
> >> > advait
> >> > On Sat, Nov 27, 2010 at 3:48 PM, Robert J. Hijmans <
> [hidden email]>
> >> > wrote:
> >> >>
> >> >> Advait,
> >> >>
> >> >> That suggests that you have missing values in the lon or lat
> >> >> dimensions. Is that true? That would be awkward, but perhaps
> something
> >> >> else is going on. Can you send me your file? If not can you send the
> >> >> results of this:
> >> >>
> >> >> filename <- "regcm.cumul.hdd.10k.96_05.new.nc"
> >> >> nc <- open.ncdf(filename)
> >> >> length(nc$dim[[1]]$vals)
> >> >> length(nc$dim[[2]]$vals)
> >> >> nc$dim[[1]]$vals
> >> >> nc$dim[[2]]$vals
> >> >> close.ncdf(nc)
> >> >>
> >> >> Thanks, Robert
> >> >>
> >> >> On Sat, Nov 27, 2010 at 8:36 AM, Advait Godbole
> >> >> <[hidden email]>
> >> >> wrote:
> >> >> > Thanks everyone for these leads..However, I am running into
> problems
> >> >> > with
> >> >> > data handling. I used ncdf to extract the two variables (507 x 356
> x
> >> >> > 120
> >> >> > -
> >> >> > lon,lat,time) and converted them into array objects. The
> suggestions
> >> >> > provided all utilize the raster package, but using "brick" on the
> >> >> > file
> >> >> > returns this error: "rcmraster <-
> >> >> > brick("regcm.cumul.hdd.10k.96_05.new.nc")
> >> >> > Error in if (e@xmin > -400 & e@xmax < 400 & e@ymin > -90.1 &
> e@ymax <
> >> >> >  :
> >> >> > missing value where TRUE/FALSE needed"
> >> >> >
> >> >> > I have a fair number of NAs (all missing values in the original
> >> >> > variables).
> >> >> > Is there a way to set a flag for that? Or any other workaround?
> >> >> > Thanks much,
> >> >> > advait
> >> >> > On Fri, Nov 26, 2010 at 6:25 PM, Robert J. Hijmans
> >> >> > <[hidden email]>
> >> >> > wrote:
> >> >> >>
> >> >> >> 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
> >> >> >> <[hidden email]>
> >> >> >> 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
> >> >> >> > <[hidden email]> 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 <[hidden email]>
> >> >> >> >> wrote:
> >> >> >> >>
> >> >> >> >> From: steven mosher <[hidden email]>
> >> >> >> >> Subject: Re: [R-sig-Geo] gridded time series analysis
> >> >> >> >> To: "Martin" <[hidden email]>
> >> >> >> >> Cc: [hidden email]
> >> >> >> >> 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 <[hidden email]
> >
> >> >> >> >> 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
> >> >> >> >>> [hidden email]
> >> >> >> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >> >> >> >>>
> >> >> >> >>
> >> >> >> >>     [[alternative HTML version deleted]]
> >> >> >> >>
> >> >> >> >> _______________________________________________
> >> >> >> >> R-sig-Geo mailing list
> >> >> >> >> [hidden email]
> >> >> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >> >> >> >>
> >> >> >> >>
> >> >> >> >>
> >> >> >> >>
> >> >> >> >>        [[alternative HTML version deleted]]
> >> >> >> >>
> >> >> >> >>
> >> >> >> >> _______________________________________________
> >> >> >> >> R-sig-Geo mailing list
> >> >> >> >> [hidden email]
> >> >> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >> >> >> >>
> >> >> >> >>
> >> >> >> >
> >> >> >>
> >> >> >> _______________________________________________
> >> >> >> R-sig-Geo mailing list
> >> >> >> [hidden email]
> >> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >> >> >
> >> >> >
> >> >> >
> >> >> > --
> >> >> > advait godbole
> >> >> >
> >> >
> >> >
> >> >
> >> > --
> >> > advait godbole
> >> >
> >
> >
> >
> > --
> > advait godbole
> >
>



--
advait godbole

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

steven mosher
 see the setMinMax()

 for the error on trim() I'm not sure.

On Thu, Dec 2, 2010 at 5:16 AM, Advait Godbole <[hidden email]>wrote:

> That worked perfectly. Thanks for pointing this out. I changed the way my
> variables were being written to the netCDF and now "brick" works. However,
> I
> am getting an error while using the stack function, like so:
>
> "Error in function (classes, fdef, mtable)  :
>  unable to find an inherited method for function "trim", for signature
> "integer"
>
> The two rasters show the following:
> > rcmraster
> class       : RasterBrick
> filename    : regcm.cumul.hdd.10k.96_05.fixed.latlonvars.nc
> nlayers     : 120
> nrow        : 356
> ncol        : 507
> ncell       : 180492
> projection  : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> min value   : ? ? ? ? ? ? ? ? ? ?
> max value   : ? ? ? ? ? ? ? ? ? ?
> xmin        : -137.3312
> xmax        : -62.11259
> ymin        : 22.12693
> ymax        : 57.35954
> xres        : 0.1483602
> yres        : 0.09896801
>
> > emitraster
> class       : RasterBrick
> filename    :
> flip.percap.vulcan.10k.dp.build.381.RES.EIAmy.96_05.this.latest.nc
> nlayers     : 120
> nrow        : 356
> ncol        : 507
> ncell       : 180492
> projection  : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> min value   : ? ? ? ? ? ? ? ? ? ?
> max value   : ? ? ? ? ? ? ? ? ? ?
> xmin        : -137.3312
> xmax        : -62.11259
> ymin        : 22.12693
> ymax        : 57.35954
> xres        : 0.1483602
> yres        : 0.09896801
>
> A getValues call shows a lot of NAs. Now, my original file does have a
> large
> number of missing values so that is why this should happen, but I suspect
> something is going wrong since the min and max values are not being
> displayed for either of these rasters. Is that the source of the problem in
> the "stack" function?
>
> Thanks,
>
> Advait
>
> On Tue, Nov 30, 2010 at 5:18 PM, Robert J. Hijmans <[hidden email]
> >wrote:
>
> > Advait,
> >
> > The is an error in the file, I think, see below:
> >
> > > library(ncdf)
> > > filename = "c:/downloads/regcm.cumul.hdd.10k.96_05.new.nc"
> > > nc <- open.ncdf(filename)
> > >
> > > zvar <- 'hdd'
> > > nc$var[[zvar]]$dim[[1]]$len
> > [1] 507
> > > nc$var[[zvar]]$dim[[2]]$len
> > [1] 356
> > >
> > > xx <- nc$var[[zvar]]$dim[[1]]$vals
> > > # should be a vector of 507 longitudes but:
> > > dim(xx)
> > [1] 507 356 120
> > >
> > > # works anyway
> > > xrange <- c(min(xx), max(xx))
> > > xrange
> > [1] -137.25705  -62.18677
> > >
> > > yy <- nc$var[[zvar]]$dim[[2]]$vals
> > > # should be a vector of 356 latitudes but:
> > > dim(yy)
> > [1] 507 356 120
> > >
> > > # this fails
> > > yrange <- c(min(yy), max(yy))
> > > yrange
> > [1] NaN NaN
> > > # and that is the only reason why raster failed.
> > >
> > > #OK
> > > yy[,,1][1:10,350:355]
> >          [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
> >  [1,] 51.49148 51.57122 51.65092 51.73055 51.81014 51.88967
> >  [2,] 51.52888 51.60866 51.68839 51.76807 51.84769 51.92727
> >  [3,] 51.56616 51.64598 51.72575 51.80547 51.88513 51.96474
> >  [4,] 51.60332 51.68318 51.76299 51.84275 51.92245 52.00210
> >  [5,] 51.64036 51.72026 51.80011 51.87991 51.95965 52.03934
> >  [6,] 51.67728 51.75722 51.83711 51.91695 51.99673 52.07646
> >  [7,] 51.71409 51.79407 51.87399 51.95387 52.03369 52.11345
> >  [8,] 51.75077 51.83079 51.91076 51.99067 52.07053 52.15033
> >  [9,] 51.78734 51.86739 51.94740 52.02735 52.10725 52.18709
> > [10,] 51.82378 51.90388 51.98392 52.06391 52.14385 52.22373
> > >
> > > #bad values here (and note the NaN):
> > > yy[,,120][1:10,350:355]
> >                [,1]           [,2]           [,3]           [,4]
> >     [,5]           [,6]
> >  [1,] -1.313847e-229 -4.895880e-272 -7.553765e-266   1.293814e+58
> > 9.551441e-293  -3.291580e-36
> >  [2,] -5.321424e+272 -3.076997e-223  -1.515780e-33 -9.895836e+171
> > 3.115686e-149  -4.970539e+92
> >  [3,]  5.259277e-203   5.980335e-56  2.893449e-212  2.974200e+274
> > 3.186547e-218  2.503689e+296
> >  [4,] -3.977129e+225 -3.696719e-283   1.284612e+42   1.887923e+28
> > -2.322520e+268  5.005820e-108
> >  [5,]  -6.569744e+88   1.824506e-74   7.663356e+39  1.608759e-201
> > -1.433761e+110 -8.804183e-211
> >  [6,] -4.237822e-162            NaN  7.405904e-275  5.859877e-258
> > 1.339820e+110  1.631993e+229
> >  [7,]   5.603325e-64   4.188764e+37 -2.188887e+243 -4.657907e+280
> > 1.795268e-185  1.166057e-175
> >  [8,] -2.112804e-229  4.536625e+137  2.910322e-278  4.407645e-274
> > 8.886205e-239 -1.913860e+238
> >  [9,]  -9.126350e+44 -1.866422e+177 -2.578975e+256  3.275379e+208
> > -2.657457e+75  -1.258259e-90
> > [10,]  6.698325e+216 -2.958355e+137 -1.058748e-178 -1.776776e+215
> > -6.368807e+78 -5.337653e+299
> > >
> > > # this would have worked
> > > yrange <- c(min(yy[,,1]), max(yy[,,1]))
> > > yrange
> > [1] 22.17641 57.31006
> > >
> > > close.ncdf(nc)
> >
> >
> > Best, Robert
> >
> >
> > On Tue, Nov 30, 2010 at 10:37 AM, Advait Godbole
> > <[hidden email]> wrote:
> > > Dear Robert,
> > > here is the link to the file:
> > > http://dl.dropbox.com/u/4030944/regcm.cumul.hdd.10k.96_05.new.nc. It
> is
> > a
> > > large file ~400 MB. Let me know if dropbox is not convenient and you
> > would
> > > rather do ftp -  but our machines are down for maintenance and I may
> not
> > be
> > > immediately able to get access.
> > > I was going over the NCL code that processes the files but did not find
> > > anything glaringly amiss, so I will have to look at it in greater
> detail
> > > since it is clear that the number of entries in the variables is not
> > right
> > > as you've pointed out. There is a chance that I have not written out
> the
> > > co-ordinates properly.
> > > Thanks,
> > > Advait
> > > On Mon, Nov 29, 2010 at 1:14 PM, Robert J. Hijmans <
> [hidden email]>
> > > wrote:
> > >>
> > >> Advait,
> > >> I do not think I can help you further unless you send me your file;
> > >> but this suggests that it is either not following the CF standards, or
> > >> something is wrong with the file, or that there is something I have
> > >> not understood about these files.
> > >> You said you have 507 x 356 x 120 values. I would expect
> > >> length(frcm$dim[[1]]$vals) == 507 and length(frcm$dim[[2]]$vals) ==
> > >> 356
> > >> Robert
> > >>
> > >> On Mon, Nov 29, 2010 at 6:11 AM, Advait Godbole <
> > [hidden email]>
> > >> wrote:
> > >> > sorry to disrupt the flow of the conversation, but after doing what
> > >> > Robert
> > >> > suggested, this is what was returned:
> > >> >> length(frcm$dim[[1]]$vals)
> > >> > [1] 21659040
> > >> >> length(frcm$dim[[2]]$vals)
> > >> > [1] 21659040
> > >> > the other two functions return the latitude and longitude arrays not
> > >> > populated with NAs/missing values. I havent reproduced them here
> > because
> > >> > of
> > >> > their number.
> > >> > Regards,
> > >> > advait
> > >> > On Sat, Nov 27, 2010 at 3:48 PM, Robert J. Hijmans <
> > [hidden email]>
> > >> > wrote:
> > >> >>
> > >> >> Advait,
> > >> >>
> > >> >> That suggests that you have missing values in the lon or lat
> > >> >> dimensions. Is that true? That would be awkward, but perhaps
> > something
> > >> >> else is going on. Can you send me your file? If not can you send
> the
> > >> >> results of this:
> > >> >>
> > >> >> filename <- "regcm.cumul.hdd.10k.96_05.new.nc"
> > >> >> nc <- open.ncdf(filename)
> > >> >> length(nc$dim[[1]]$vals)
> > >> >> length(nc$dim[[2]]$vals)
> > >> >> nc$dim[[1]]$vals
> > >> >> nc$dim[[2]]$vals
> > >> >> close.ncdf(nc)
> > >> >>
> > >> >> Thanks, Robert
> > >> >>
> > >> >> On Sat, Nov 27, 2010 at 8:36 AM, Advait Godbole
> > >> >> <[hidden email]>
> > >> >> wrote:
> > >> >> > Thanks everyone for these leads..However, I am running into
> > problems
> > >> >> > with
> > >> >> > data handling. I used ncdf to extract the two variables (507 x
> 356
> > x
> > >> >> > 120
> > >> >> > -
> > >> >> > lon,lat,time) and converted them into array objects. The
> > suggestions
> > >> >> > provided all utilize the raster package, but using "brick" on the
> > >> >> > file
> > >> >> > returns this error: "rcmraster <-
> > >> >> > brick("regcm.cumul.hdd.10k.96_05.new.nc")
> > >> >> > Error in if (e@xmin > -400 & e@xmax < 400 & e@ymin > -90.1 &
> > e@ymax <
> > >> >> >  :
> > >> >> > missing value where TRUE/FALSE needed"
> > >> >> >
> > >> >> > I have a fair number of NAs (all missing values in the original
> > >> >> > variables).
> > >> >> > Is there a way to set a flag for that? Or any other workaround?
> > >> >> > Thanks much,
> > >> >> > advait
> > >> >> > On Fri, Nov 26, 2010 at 6:25 PM, Robert J. Hijmans
> > >> >> > <[hidden email]>
> > >> >> > wrote:
> > >> >> >>
> > >> >> >> 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
> > >> >> >> <[hidden email]>
> > >> >> >> 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
> > >> >> >> > <[hidden email]> 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 <[hidden email]>
> > >> >> >> >> wrote:
> > >> >> >> >>
> > >> >> >> >> From: steven mosher <[hidden email]>
> > >> >> >> >> Subject: Re: [R-sig-Geo] gridded time series analysis
> > >> >> >> >> To: "Martin" <[hidden email]>
> > >> >> >> >> Cc: [hidden email]
> > >> >> >> >> 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 <
> [hidden email]
> > >
> > >> >> >> >> 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
> > >> >> >> >>> [hidden email]
> > >> >> >> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > >> >> >> >>>
> > >> >> >> >>
> > >> >> >> >>     [[alternative HTML version deleted]]
> > >> >> >> >>
> > >> >> >> >> _______________________________________________
> > >> >> >> >> R-sig-Geo mailing list
> > >> >> >> >> [hidden email]
> > >> >> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > >> >> >> >>
> > >> >> >> >>
> > >> >> >> >>
> > >> >> >> >>
> > >> >> >> >>        [[alternative HTML version deleted]]
> > >> >> >> >>
> > >> >> >> >>
> > >> >> >> >> _______________________________________________
> > >> >> >> >> R-sig-Geo mailing list
> > >> >> >> >> [hidden email]
> > >> >> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > >> >> >> >>
> > >> >> >> >>
> > >> >> >> >
> > >> >> >>
> > >> >> >> _______________________________________________
> > >> >> >> R-sig-Geo mailing list
> > >> >> >> [hidden email]
> > >> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > >> >> >
> > >> >> >
> > >> >> >
> > >> >> > --
> > >> >> > advait godbole
> > >> >> >
> > >> >
> > >> >
> > >> >
> > >> > --
> > >> > advait godbole
> > >> >
> > >
> > >
> > >
> > > --
> > > advait godbole
> > >
> >
>
>
>
> --
> advait godbole
>
>        [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: gridded time series analysis

Robert Hijmans
The question marks indicate that the min and max values are not known
to the object (because the file type does not store them, and it is
too 'expensive' to determine them on the fly). setMinMax if you really
need to know.
The error on trim seems to be a problem with layerNames being numeric.
I'll fix that. This may be a work-around






On Thu, Dec 2, 2010 at 10:23 AM, steven mosher <[hidden email]> wrote:

>  see the setMinMax()
>
>  for the error on trim() I'm not sure.
>
> On Thu, Dec 2, 2010 at 5:16 AM, Advait Godbole <[hidden email]>
> wrote:
>>
>> That worked perfectly. Thanks for pointing this out. I changed the way my
>> variables were being written to the netCDF and now "brick" works. However,
>> I
>> am getting an error while using the stack function, like so:
>>
>> "Error in function (classes, fdef, mtable)  :
>>  unable to find an inherited method for function "trim", for signature
>> "integer"
>>
>> The two rasters show the following:
>> > rcmraster
>> class       : RasterBrick
>> filename    : regcm.cumul.hdd.10k.96_05.fixed.latlonvars.nc
>> nlayers     : 120
>> nrow        : 356
>> ncol        : 507
>> ncell       : 180492
>> projection  : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
>> min value   : ? ? ? ? ? ? ? ? ? ?
>> max value   : ? ? ? ? ? ? ? ? ? ?
>> xmin        : -137.3312
>> xmax        : -62.11259
>> ymin        : 22.12693
>> ymax        : 57.35954
>> xres        : 0.1483602
>> yres        : 0.09896801
>>
>> > emitraster
>> class       : RasterBrick
>> filename    :
>> flip.percap.vulcan.10k.dp.build.381.RES.EIAmy.96_05.this.latest.nc
>> nlayers     : 120
>> nrow        : 356
>> ncol        : 507
>> ncell       : 180492
>> projection  : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
>> min value   : ? ? ? ? ? ? ? ? ? ?
>> max value   : ? ? ? ? ? ? ? ? ? ?
>> xmin        : -137.3312
>> xmax        : -62.11259
>> ymin        : 22.12693
>> ymax        : 57.35954
>> xres        : 0.1483602
>> yres        : 0.09896801
>>
>> A getValues call shows a lot of NAs. Now, my original file does have a
>> large
>> number of missing values so that is why this should happen, but I suspect
>> something is going wrong since the min and max values are not being
>> displayed for either of these rasters. Is that the source of the problem
>> in
>> the "stack" function?
>>
>> Thanks,
>>
>> Advait
>>
>> On Tue, Nov 30, 2010 at 5:18 PM, Robert J. Hijmans
>> <[hidden email]>wrote:
>>
>> > Advait,
>> >
>> > The is an error in the file, I think, see below:
>> >
>> > > library(ncdf)
>> > > filename = "c:/downloads/regcm.cumul.hdd.10k.96_05.new.nc"
>> > > nc <- open.ncdf(filename)
>> > >
>> > > zvar <- 'hdd'
>> > > nc$var[[zvar]]$dim[[1]]$len
>> > [1] 507
>> > > nc$var[[zvar]]$dim[[2]]$len
>> > [1] 356
>> > >
>> > > xx <- nc$var[[zvar]]$dim[[1]]$vals
>> > > # should be a vector of 507 longitudes but:
>> > > dim(xx)
>> > [1] 507 356 120
>> > >
>> > > # works anyway
>> > > xrange <- c(min(xx), max(xx))
>> > > xrange
>> > [1] -137.25705  -62.18677
>> > >
>> > > yy <- nc$var[[zvar]]$dim[[2]]$vals
>> > > # should be a vector of 356 latitudes but:
>> > > dim(yy)
>> > [1] 507 356 120
>> > >
>> > > # this fails
>> > > yrange <- c(min(yy), max(yy))
>> > > yrange
>> > [1] NaN NaN
>> > > # and that is the only reason why raster failed.
>> > >
>> > > #OK
>> > > yy[,,1][1:10,350:355]
>> >          [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
>> >  [1,] 51.49148 51.57122 51.65092 51.73055 51.81014 51.88967
>> >  [2,] 51.52888 51.60866 51.68839 51.76807 51.84769 51.92727
>> >  [3,] 51.56616 51.64598 51.72575 51.80547 51.88513 51.96474
>> >  [4,] 51.60332 51.68318 51.76299 51.84275 51.92245 52.00210
>> >  [5,] 51.64036 51.72026 51.80011 51.87991 51.95965 52.03934
>> >  [6,] 51.67728 51.75722 51.83711 51.91695 51.99673 52.07646
>> >  [7,] 51.71409 51.79407 51.87399 51.95387 52.03369 52.11345
>> >  [8,] 51.75077 51.83079 51.91076 51.99067 52.07053 52.15033
>> >  [9,] 51.78734 51.86739 51.94740 52.02735 52.10725 52.18709
>> > [10,] 51.82378 51.90388 51.98392 52.06391 52.14385 52.22373
>> > >
>> > > #bad values here (and note the NaN):
>> > > yy[,,120][1:10,350:355]
>> >                [,1]           [,2]           [,3]           [,4]
>> >     [,5]           [,6]
>> >  [1,] -1.313847e-229 -4.895880e-272 -7.553765e-266   1.293814e+58
>> > 9.551441e-293  -3.291580e-36
>> >  [2,] -5.321424e+272 -3.076997e-223  -1.515780e-33 -9.895836e+171
>> > 3.115686e-149  -4.970539e+92
>> >  [3,]  5.259277e-203   5.980335e-56  2.893449e-212  2.974200e+274
>> > 3.186547e-218  2.503689e+296
>> >  [4,] -3.977129e+225 -3.696719e-283   1.284612e+42   1.887923e+28
>> > -2.322520e+268  5.005820e-108
>> >  [5,]  -6.569744e+88   1.824506e-74   7.663356e+39  1.608759e-201
>> > -1.433761e+110 -8.804183e-211
>> >  [6,] -4.237822e-162            NaN  7.405904e-275  5.859877e-258
>> > 1.339820e+110  1.631993e+229
>> >  [7,]   5.603325e-64   4.188764e+37 -2.188887e+243 -4.657907e+280
>> > 1.795268e-185  1.166057e-175
>> >  [8,] -2.112804e-229  4.536625e+137  2.910322e-278  4.407645e-274
>> > 8.886205e-239 -1.913860e+238
>> >  [9,]  -9.126350e+44 -1.866422e+177 -2.578975e+256  3.275379e+208
>> > -2.657457e+75  -1.258259e-90
>> > [10,]  6.698325e+216 -2.958355e+137 -1.058748e-178 -1.776776e+215
>> > -6.368807e+78 -5.337653e+299
>> > >
>> > > # this would have worked
>> > > yrange <- c(min(yy[,,1]), max(yy[,,1]))
>> > > yrange
>> > [1] 22.17641 57.31006
>> > >
>> > > close.ncdf(nc)
>> >
>> >
>> > Best, Robert
>> >
>> >
>> > On Tue, Nov 30, 2010 at 10:37 AM, Advait Godbole
>> > <[hidden email]> wrote:
>> > > Dear Robert,
>> > > here is the link to the file:
>> > > http://dl.dropbox.com/u/4030944/regcm.cumul.hdd.10k.96_05.new.nc. It
>> > > is
>> > a
>> > > large file ~400 MB. Let me know if dropbox is not convenient and you
>> > would
>> > > rather do ftp -  but our machines are down for maintenance and I may
>> > > not
>> > be
>> > > immediately able to get access.
>> > > I was going over the NCL code that processes the files but did not
>> > > find
>> > > anything glaringly amiss, so I will have to look at it in greater
>> > > detail
>> > > since it is clear that the number of entries in the variables is not
>> > right
>> > > as you've pointed out. There is a chance that I have not written out
>> > > the
>> > > co-ordinates properly.
>> > > Thanks,
>> > > Advait
>> > > On Mon, Nov 29, 2010 at 1:14 PM, Robert J. Hijmans
>> > > <[hidden email]>
>> > > wrote:
>> > >>
>> > >> Advait,
>> > >> I do not think I can help you further unless you send me your file;
>> > >> but this suggests that it is either not following the CF standards,
>> > >> or
>> > >> something is wrong with the file, or that there is something I have
>> > >> not understood about these files.
>> > >> You said you have 507 x 356 x 120 values. I would expect
>> > >> length(frcm$dim[[1]]$vals) == 507 and length(frcm$dim[[2]]$vals) ==
>> > >> 356
>> > >> Robert
>> > >>
>> > >> On Mon, Nov 29, 2010 at 6:11 AM, Advait Godbole <
>> > [hidden email]>
>> > >> wrote:
>> > >> > sorry to disrupt the flow of the conversation, but after doing what
>> > >> > Robert
>> > >> > suggested, this is what was returned:
>> > >> >> length(frcm$dim[[1]]$vals)
>> > >> > [1] 21659040
>> > >> >> length(frcm$dim[[2]]$vals)
>> > >> > [1] 21659040
>> > >> > the other two functions return the latitude and longitude arrays
>> > >> > not
>> > >> > populated with NAs/missing values. I havent reproduced them here
>> > because
>> > >> > of
>> > >> > their number.
>> > >> > Regards,
>> > >> > advait
>> > >> > On Sat, Nov 27, 2010 at 3:48 PM, Robert J. Hijmans <
>> > [hidden email]>
>> > >> > wrote:
>> > >> >>
>> > >> >> Advait,
>> > >> >>
>> > >> >> That suggests that you have missing values in the lon or lat
>> > >> >> dimensions. Is that true? That would be awkward, but perhaps
>> > something
>> > >> >> else is going on. Can you send me your file? If not can you send
>> > >> >> the
>> > >> >> results of this:
>> > >> >>
>> > >> >> filename <- "regcm.cumul.hdd.10k.96_05.new.nc"
>> > >> >> nc <- open.ncdf(filename)
>> > >> >> length(nc$dim[[1]]$vals)
>> > >> >> length(nc$dim[[2]]$vals)
>> > >> >> nc$dim[[1]]$vals
>> > >> >> nc$dim[[2]]$vals
>> > >> >> close.ncdf(nc)
>> > >> >>
>> > >> >> Thanks, Robert
>> > >> >>
>> > >> >> On Sat, Nov 27, 2010 at 8:36 AM, Advait Godbole
>> > >> >> <[hidden email]>
>> > >> >> wrote:
>> > >> >> > Thanks everyone for these leads..However, I am running into
>> > problems
>> > >> >> > with
>> > >> >> > data handling. I used ncdf to extract the two variables (507 x
>> > >> >> > 356
>> > x
>> > >> >> > 120
>> > >> >> > -
>> > >> >> > lon,lat,time) and converted them into array objects. The
>> > suggestions
>> > >> >> > provided all utilize the raster package, but using "brick" on
>> > >> >> > the
>> > >> >> > file
>> > >> >> > returns this error: "rcmraster <-
>> > >> >> > brick("regcm.cumul.hdd.10k.96_05.new.nc")
>> > >> >> > Error in if (e@xmin > -400 & e@xmax < 400 & e@ymin > -90.1 &
>> > e@ymax <
>> > >> >> >  :
>> > >> >> > missing value where TRUE/FALSE needed"
>> > >> >> >
>> > >> >> > I have a fair number of NAs (all missing values in the original
>> > >> >> > variables).
>> > >> >> > Is there a way to set a flag for that? Or any other workaround?
>> > >> >> > Thanks much,
>> > >> >> > advait
>> > >> >> > On Fri, Nov 26, 2010 at 6:25 PM, Robert J. Hijmans
>> > >> >> > <[hidden email]>
>> > >> >> > wrote:
>> > >> >> >>
>> > >> >> >> 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
>> > >> >> >> <[hidden email]>
>> > >> >> >> 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
>> > >> >> >> > <[hidden email]> 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 <[hidden email]>
>> > >> >> >> >> wrote:
>> > >> >> >> >>
>> > >> >> >> >> From: steven mosher <[hidden email]>
>> > >> >> >> >> Subject: Re: [R-sig-Geo] gridded time series analysis
>> > >> >> >> >> To: "Martin" <[hidden email]>
>> > >> >> >> >> Cc: [hidden email]
>> > >> >> >> >> 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
>> > >> >> >> >> <[hidden email]
>> > >
>> > >> >> >> >> 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
>> > >> >> >> >>> [hidden email]
>> > >> >> >> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> > >> >> >> >>>
>> > >> >> >> >>
>> > >> >> >> >>     [[alternative HTML version deleted]]
>> > >> >> >> >>
>> > >> >> >> >> _______________________________________________
>> > >> >> >> >> R-sig-Geo mailing list
>> > >> >> >> >> [hidden email]
>> > >> >> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> > >> >> >> >>
>> > >> >> >> >>
>> > >> >> >> >>
>> > >> >> >> >>
>> > >> >> >> >>        [[alternative HTML version deleted]]
>> > >> >> >> >>
>> > >> >> >> >>
>> > >> >> >> >> _______________________________________________
>> > >> >> >> >> R-sig-Geo mailing list
>> > >> >> >> >> [hidden email]
>> > >> >> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> > >> >> >> >>
>> > >> >> >> >>
>> > >> >> >> >
>> > >> >> >>
>> > >> >> >> _______________________________________________
>> > >> >> >> R-sig-Geo mailing list
>> > >> >> >> [hidden email]
>> > >> >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> > >> >> >
>> > >> >> >
>> > >> >> >
>> > >> >> > --
>> > >> >> > advait godbole
>> > >> >> >
>> > >> >
>> > >> >
>> > >> >
>> > >> > --
>> > >> > advait godbole
>> > >> >
>> > >
>> > >
>> > >
>> > > --
>> > > advait godbole
>> > >
>> >
>>
>>
>>
>> --
>> advait godbole
>>
>>        [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
123