[R-sig-Geo] gridded time series analysis
Robert J. Hijmans
r.hijmans at gmail.com
Tue Nov 30 23:18:27 CET 2010
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
>
More information about the R-sig-Geo
mailing list