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