[R-sig-Geo] gridded time series analysis

[Robert J. 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 <advaitgodbole at gmail.com> 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 at xmin > -400 & e at xmax < 400 & e at ymin > -90.1 & e at 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 <r.hijmans at gmail.com>
> 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 <r.hijmans at gmail.com>
>> wrote:
>> > It seems that 'calc' does not like this (any more; another thing to
>> > fix) . If your rasters are not very large you can use apply, which
>> > makes it much faster:
>> >
>> > library(raster)
>> > #creating some data
>> > r <- raster(nrow=10, ncol=10)
>> > s <- list()
>> > for (i in 1:25) { s[i] <- setValues(r, rnorm(ncell(r), i, 3) ) }
>> > s <- stack(s)
>> >
>> > # regression
>> > time <- 1:nlayers(s)
>> > fun <- function(x) { lm(x ~ time)$coefficients[2] }
>> > r[] <- apply(getValues(s), 1, fun)
>> >
>> > Robert
>> >
>> >
>> >
>> > On Fri, Nov 26, 2010 at 2:51 PM, Jacob van Etten
>> > <jacobvanetten at yahoo.com> wrote:
>> >> you could try this approach (use calc whenever you can):
>> >>
>> >> (supposing your bricks have 12 layers)
>> >>
>> >> br3 <- stack(brick1, brick2)
>> >> lmS <- function(x) lm(x[1:12] ~ x[13:24)$coefficients[2]
>> >> r <- calc(br3, lmS)
>> >>
>> >> Jacob.
>> >>
>> >> --- On Fri, 26/11/10, steven mosher <moshersteven at gmail.com> wrote:
>> >>
>> >> From: steven mosher <moshersteven at gmail.com>
>> >> Subject: Re: [R-sig-Geo] gridded time series analysis
>> >> To: "Martin" <martin_brandt at gmx.net>
>> >> Cc: r-sig-geo at stat.math.ethz.ch
>> >> Date: Friday, 26 November, 2010, 23:33
>> >>
>> >> that's cool, I'm also interested in a similar problem but just with one
>> >> brick ending up with a slope raster as the output. It may be possible
>> >> with
>> >> stackApply(). have a look. or maybe robert will chime in
>> >>
>> >>
>> >>
>> >> On Fri, Nov 26, 2010 at 1:35 PM, Martin <martin_brandt at gmx.net> wrote:
>> >>
>> >>>
>> >>> this is what I did to perform a regression between two bricks (each
>> >>> brick
>> >>> represent a time series):
>> >>>
>> >>> r <- raster(brick1)
>> >>> for (i in 1:ncell(r)) {
>> >>> r[i] = lm(as.ts(cellValues(brick1, i)) ~ as.ts(cellValues(brick2,
>> >>> i)))$coefficients[2]
>> >>> }
>> >>>
>> >>> The result will be a slope raster, but it really takes a lot of time,
>> >>> so
>> >>> maybe there is a better solution..
>> >>>
>> >>>
>> >>> --
>> >>> View this message in context:
>> >>>
>> >>> http://r-sig-geo.2731867.n2.nabble.com/gridded-time-series-analysis-tp5775651p5778472.html
>> >>> Sent from the R-sig-geo mailing list archive at Nabble.com.
>> >>>
>> >>> _______________________________________________
>> >>> R-sig-Geo mailing list
>> >>> R-sig-Geo at stat.math.ethz.ch
>> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >>>
>> >>
>> >>     [[alternative HTML version deleted]]
>> >>
>> >> _______________________________________________
>> >> R-sig-Geo mailing list
>> >> R-sig-Geo at stat.math.ethz.ch
>> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >>
>> >>
>> >>
>> >>
>> >>        [[alternative HTML version deleted]]
>> >>
>> >>
>> >> _______________________________________________
>> >> R-sig-Geo mailing list
>> >> R-sig-Geo at stat.math.ethz.ch
>> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >>
>> >>
>> >
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
> --
> advait godbole
>