[R-sig-Geo] gridded time series analysis
Robert J. Hijmans
r.hijmans at gmail.com
Thu Dec 2 19:43:50 CET 2010
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
>>
>>
>
More information about the R-sig-Geo
mailing list