[R-sig-Geo] gridded time series analysis
Robert J. Hijmans
r.hijmans at gmail.com
Sun Dec 5 23:29:06 CET 2010
>If so, should I use a for loop instead?
No, that is the point of the calc function, that you do not do that.
Internally, calc uses apply as in
apply(values, 1, fun)
On Sat, Dec 4, 2010 at 1:50 PM, Advait Godbole <advaitgodbole at gmail.com> wrote:
> Roger, that fix worked, its there in the latest version. Thanks much.
> I am running into the problem that Martin mentioned earlier in the thread:
> Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
> 0 (non-NA) cases
>
> For which you suggested this:
> #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)
>
> What I am confused is that my stack "s" that I will pass via calc is of
> dimensions 356 x 507 x 240, and therefore is.na will act on the entire array
> and therefore still fail at those points where it is not NA. Is that right?
> If so, should I use a for loop instead?
> Thanks,
> Advait
> On Fri, Dec 3, 2010 at 2:58 PM, Robert J. Hijmans <r.hijmans at gmail.com>
> wrote:
>>
>> 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
>> >
>> >
>
>
>
> --
> advait godbole
>
More information about the R-sig-Geo
mailing list