[R-sig-Geo] gridded time series analysis
Robert J. Hijmans
r.hijmans at gmail.com
Thu Dec 23 01:26:37 CET 2010
Advait,
When asking for help, please provide a self-contained example whenever
possible, as below.
library(raster)
v1 <- matrix(rep(1:12, 16), nrow=16)
v2 <- matrix(rep(1:12, 16), nrow=16, byrow=TRUE)
b <- brick(ncol=4, nrow=4)
b1 <- setValues(b, v1)
b2 <- setValues(b, v2)
s <- stack(b1, b2)
fun <- function(x) {
resp = x[1:12]
expl = x[13:24]
mdl = lm(resp ~ expl)
res = c(mdl$coef[1] , mdl$coef[2])
return(res)
}
# now do:
b2 <- calc(s, fun)
# as an alternative to the above you could do
# but I would not recommend this
datCoef = apply(getValues(s), 1, fun2)
b1 <- setValues(s, t(datCoef))
# I think you were complicating matters with plyr.
# But to get those values in a brick you could do something like
b1 <- setValues(s, as.matrix(datCoef)[,3:4])
# But I am not sure if you would get the values in the correct cells
# (the number matches, but perhaps not the order).
Best, Robert
On Wed, Dec 22, 2010 at 10:39 AM, Advait Godbole
<advaitgodbole at gmail.com> wrote:
> Dear Members,
> I was able to find a solution after some deliberation and posting it here
> for the benefit of forum members.
> To recap, my data was in the form of 2 netCDF files, each of dimension 356 x
> 507 x 120 (lat,lon,time). I wanted to perform regressions between the two
> variables on each grid point over time. I used "brick" to read in the
> netcdfs as raster layers and stack along with the "plyr" package to perform
> these regressions on the grid.
> # create raster object from the netcdf
> rcmraster <- brick("regcm.cumul.hdd.10k.96_05.fixed.latlonvars.nc", varname
> = "hdd")
> emitraster <-
> brick("flip.percap.vulcan.10k.dp.build.381.RES.EIAmy.96_05.this.latest.nc",varname="emitpercap")
> # crop subset of above rasters, test method on subset which is 4 x 4 x 120
> lim <- extent(rcmraster,180,185,300,305)
> rcmraster_crop <- crop(rcmraster,lim)
> emitraster_crop <- crop(emitraster,lim)
> # do some plotting to check
> par(mfrow=c(1,2))
> plot(rcmraster_crop[[1]])
> plot(emitraster_crop[[1]])
> # test regressions for cropped subset
> s_crop <- stack(emitraster_crop,rcmraster_crop)
> # use of plyr (provides a more elegant alternative to lapply)
> # the data frame datCoef contains the slope and intercept for each point in
> the grid (denoted by X1 and X2)
> datCoef = adply( as.array(s_crop) , 1:2 ,
> function(x ) {
> resp = x[1:120]
> expl = x[121:240]
> mdl = lm(resp ~ expl)
> res = c(mdl$coef[1] , mdl$coef[2])
> return(res)
>
> }
> )
> Output of datCoef()
>> datCoef
> X1 X2 (Intercept) expl
> 1 1 1 0.0010523406 2.623927e-05
> 2 2 1 0.0019401572 5.049197e-05
> 3 3 1 0.0171326653 4.425026e-04
> 4 4 1 0.0012244689 3.045755e-05
> 5 1 2 0.0011454348 2.892668e-05
> 6 2 2 0.0013920425 3.460922e-05
> 7 3 2 0.0015973400 3.943349e-05
> 8 4 2 0.0016494511 4.045801e-05
> 9 1 3 0.0018757258 4.704107e-05
> 10 2 3 0.0038066075 9.378512e-05
> 11 3 3 0.0039669120 9.719391e-05
> 12 4 3 0.0038055524 9.246021e-05
> 13 1 4 0.0015700329 3.921758e-05
> 14 2 4 0.0003727142 9.156119e-06
> 15 3 4 0.0005315722 1.308365e-05
> 16 4 4 0.0005872930 1.539454e-05
> Now, what I want to do is for each of these X1,X2, get a "predicted" map
> based on a new X which is an array of dimensions 4 x 4 (same as the rasters.
> in the original dataset, it will be 356 x 507). This array has the following
> values:
>> futurehdda
> [,1] [,2] [,3] [,4]
> [1,] 5707.575 5696.325 5692.960 5683.110
> [2,] 5680.303 5665.092 5660.790 5647.896
> [3,] 5656.163 5641.045 5638.886 5625.231
> [4,] 5634.604 5619.317 5620.226 5601.729
> How can I extract the slope and intercept from datCoef and use the above X
> to get a predicted Y with the same dimensions (4 x 4)?
> Thanks,
> Advait
> On Mon, Dec 6, 2010 at 4:57 PM, Robert J. Hijmans <r.hijmans at gmail.com>
> wrote:
>>
>> Dear Advait,
>>
>> > What I dont understand is:
>> > 1) How can each slope be the same (values of x1)
>>
>> Apparently the values in these cells are all the same (not that you
>> only use the first 8 time steps, not all 120)
>>
>> Here is an example with random values in which the slopes are all
>> different
>>
>> b <- brick(ncol=4, nrow=4, nl=120)
>> b <- setValues(b, matrix(runif(prod(dim(b))), ncol=120))
>> fun = function(x) { if(is.na(x[1])) { NA } else { lm(x[1:60] ~
>> x[61:120])$coefficients[2] } }
>> x <- calc(b, fun)
>> plot(x)
>> as.matrix(x)
>>
>>
>> > 2) The raster layers are of dimensions 4 x 4 x 120
>> > but dim(getValues(rcmraster_crop))
>> > [1] 16 120. Why is that? This must have implications for the way that
>> > the
>> > regressions for the brick are handled wouldnt it?
>> No
>>
>> > Or is it specific to the
>> > way "getValues" handles things?
>>
>> Yes. rows are cells (4*4=16), columns are layers.
>>
>> > Also, does anyone have any pointers of how to accomplish the NA testing
>> > and
>> See an example above, assuming that if the first time step is NA, they
>> all are. You could also do things like
>> if (sum(is.na(x)) > 60)
>> depends on your context.
>>
>>
>> Best, Robert
>>
>> On Mon, Dec 6, 2010 at 8:30 AM, Advait Godbole <advaitgodbole at gmail.com>
>> wrote:
>> > Dear Members,
>> > To recap, I initially had trouble reading in two 3D netCDFs in to R
>> > using
>> > the raster "brick" function. Each brick is a 356 x 507 x 120 (lat, lon,
>> > time) raster. With that resolved (thanks Robert!) I have the following
>> > code
>> > to perform regressions between the two time series (at each grid-cell)
>> > # read the netCDFs
>> > rcmraster <- brick("regcm.cumul.hdd.10k.96_05.fixed.latlonvars.nc",
>> > varname
>> > = "hdd")
>> > emitraster <-
>> >
>> > brick("flip.percap.vulcan.10k.dp.build.381.RES.EIAmy.96_05.this.latest.nc",varname="emitpercap")
>> > # crop subset of above rasters to test
>> > lim <- extent(rcmraster,180,185,300,305) # createss a 4 x 4 x 120 brick
>> > rcmraster_crop <- crop(rcmraster,lim)
>> > emitraster_crop <- crop(emitraster,lim)
>> > #regressions
>> > s_crop <- stack(emitraster_crop,rcmraster_crop)
>> > func <- function(x) { lm(x[1:4] ~ x[5:8])$coefficients[2] }
>> > #regression
>> > for emit* VS rcm*
>> > x1 <- calc(s_crop, func)
>> > While the above regressions run, I have to ultimately do the same for
>> > the
>> > entire domain (356 x 507 cells) and therefore need to check for NAs at
>> > all
>> > time steps for each grid cell for both the bricks, and then have lm run.
>> > To
>> > that end,, I ran the following diagnostics on the cropped dataset:
>> >> dim(rcmraster_crop)
>> > [1] 4 4 120
>> >> dim(emitraster_crop)
>> > [1] 4 4 120
>> >> dim(x1)
>> > [1] 4 4
>> >> getValues(x1)
>> > [1] 3.968207 3.968207 3.968207 3.968207 3.968207 3.968207 3.968207
>> > 3.968207
>> > 3.968207 3.968207 3.968207
>> > [12] 3.968207 3.968207 3.968207 3.968207 3.968207
>> >
>> > What I dont understand is:
>> > 1) How can each slope be the same (values of x1)
>> > 2) The raster layers are of dimensions 4 x 4 x 120
>> > but dim(getValues(rcmraster_crop))
>> > [1] 16 120. Why is that? This must have implications for the way that
>> > the
>> > regressions for the brick are handled wouldnt it? Or is it specific to
>> > the
>> > way "getValues" handles things?
>> > Also, does anyone have any pointers of how to accomplish the NA testing
>> > and
>> > regressions ?
>> > Thanks,
>> > Advait
>> >
>> > On Sun, Dec 5, 2010 at 6:29 PM, Robert J. Hijmans <r.hijmans at gmail.com>
>> > wrote:
>> >>
>> >> Advait,
>> >> It is just an example. In this case the function only checks the first
>> >> time-slice for all cells. That often works because oftentimes either
>> >> all values are NA, or none. Also, I think lm works when not all values
>> >> are NA. If your data is different, you may need a different solution.
>> >> The best thing might be to try it out with a small dataset (you can
>> >> use crop to create one).
>> >> Robert
>> >>
>> >> On Sun, Dec 5, 2010 at 3:14 PM, Advait Godbole
>> >> <advaitgodbole at gmail.com>
>> >> wrote:
>> >> > Earlier in the threade you had mentioned that for two 120 layer
>> >> > bricks,
>> >> > this
>> >> > is what I should do:
>> >> > # regression of values in one brick (or stack) with another (Jacob's
>> >> > suggestion)
>> >> > s <- stack(s1, s2)
>> >> > # s1 and s2 have 120 layers
>> >> > fun <- function(x) { lm(x[1:120] ~ x[121:240])$coefficients[2] }
>> >> > x1 <- calc(s, fun)
>> >> > So if I have to check for NAs for the multi layered raster, I would
>> >> > need
>> >> > to
>> >> > do something like
>> >> > fun <- function (x) { if(is.na(x[1:120])) {NA} else
>> >> > {lm(x[1:120]~x[121:240])$coefficients[2]}}
>> >> > calc(s,fun=fun)
>> >> > That does not seem intuitively right to me although I am not very
>> >> > conversant
>> >> > with how is.na handles 3D variables to go much further. Also, should
>> >> > I
>> >> > not
>> >> > check the other brick for NAs too?
>> >> > Advait
>> >> > On Sun, Dec 5, 2010 at 5:29 PM, Robert J. Hijmans
>> >> > <r.hijmans at gmail.com>
>> >> > wrote:
>> >> >>
>> >> >> >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
>> >> >> >
>> >> >
>> >> >
>> >> >
>> >> > --
>> >> > advait godbole
>> >> >
>> >
>> >
>> >
>> > --
>> > advait godbole
>> >
>
>
>
> --
> advait godbole
>
More information about the R-sig-Geo
mailing list