[R-sig-Geo] gridded time series analysis

Robert J. Hijmans r.hijmans at gmail.com
Fri Dec 3 20:58:19 CET 2010 

It may have been fixed. Please try this again in raster 1.7-8 (now on
CRAN for windows and linux, mac in a day or so). Robert

On Fri, Dec 3, 2010 at 11:54 AM, steven mosher <moshersteven at gmail.com> wrote:
> Well,
>   I need some more information.
>  before you create the stack do this and report the results
> test1 <- layerNames(rcmraster)
> test2 <- layerNames(emitraster)
> str(test1)
> str(test2)
> Its a bit of a hack but  you can try to copy the layerNames out of the
> invidual stacks.
> Then set them to    ""
> then create one big stack
> then set the layerNames.
> ( I recall having to do that once in the past after a stackApply())
> I know in the past, I've had some similar issues, but I thought that they
>  had been fixed.
> I'll try to take a look at the code  in a bit, but if there is a bug, either
> robert or Jacob would have to fix it.
>
> On Fri, Dec 3, 2010 at 4:28 AM, Advait Godbole <advaitgodbole at gmail.com>
> wrote:
>>
>> With Robert being away until the 15th, can someone else help out with this
>> issue? To recap, I have read in two 507 x  356 x 120 (lon,lat,time) netCDF
>> files into R using the "brick" function. I get the following error when I
>> use "stack" to clump these two multi layer rasters together.
>> > s <- stack(rcmraster,emitraster)
>> Error in function (classes, fdef, mtable)  :
>>   unable to find an inherited method for function "trim", for signature
>> "integer"
>> Robert suggested that it is because of the layerNames being numeric, annd
>> that layerNames() <- layerNames() would solve it but it did not help, even
>> though layerNames() is now of type character.
>> Thanks much!
>> adv
>> On Fri, Dec 3, 2010 at 7:13 AM, Advait Godbole <advaitgodbole at gmail.com>
>> wrote:
>>>
>>> the setMinMax() function worked. Thanks Steve and Robert..
>>> the work around did not work - I still get the trim() error..
>>> adv
>>>
>>> On Thu, Dec 2, 2010 at 1:43 PM, Robert J. Hijmans <r.hijmans at gmail.com>
>>> wrote:
>>>>
>>>> This is the work around:
>>>>
>>>> layerNames(rcmraster) <- layerNames(rcmraster)
>>>>
>>>>
>>>> On Thu, Dec 2, 2010 at 10:43 AM, Robert J. Hijmans <r.hijmans at gmail.com>
>>>> wrote:
>>>> > 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
>>>> > <moshersteven at gmail.com> wrote:
>>>> >>  see the setMinMax()
>>>> >>
>>>> >>  for the error on trim() I'm not sure.
>>>> >>
>>>> >> On Thu, Dec 2, 2010 at 5:16 AM, Advait Godbole
>>>> >> <advaitgodbole at gmail.com>
>>>> >> 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
>>>> >>> <r.hijmans at gmail.com>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
>>>> >>> > <advaitgodbole at gmail.com> 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
>>>> >>> > > <r.hijmans at gmail.com>
>>>> >>> > > 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 <
>>>> >>> > advaitgodbole at gmail.com>
>>>> >>> > >> 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 <
>>>> >>> > r.hijmans at gmail.com>
>>>> >>> > >> > 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
>>>> >>> > >> >> <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
>>>> >>> > >> >> >
>>>> >>> > >> >
>>>> >>> > >> >
>>>> >>> > >> >
>>>> >>> > >> > --
>>>> >>> > >> > advait godbole
>>>> >>> > >> >
>>>> >>> > >
>>>> >>> > >
>>>> >>> > >
>>>> >>> > > --
>>>> >>> > > advait godbole
>>>> >>> > >
>>>> >>> >
>>>> >>>
>>>> >>>
>>>> >>>
>>>> >>> --
>>>> >>> advait godbole
>>>> >>>
>>>> >>>        [[alternative HTML version deleted]]
>>>> >>>
>>>> >>> _______________________________________________
>>>> >>> R-sig-Geo mailing list
>>>> >>> R-sig-Geo at r-project.org
>>>> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>> >>
>>>> >>
>>>> >
>>>
>>>
>>>
>>> --
>>> advait godbole
>>
>>
>>
>> --
>> advait godbole
>
>

More information about the R-sig-Geo mailing list