[R-sig-Geo] gridded time series analysis
Robert J. Hijmans
r.hijmans at gmail.com
Mon Dec 6 00:29:45 CET 2010
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
>
More information about the R-sig-Geo
mailing list