gridded time series analysis

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

nc$var[[zvar]]$dim[[2]]$len

xx <- nc$var[[zvar]]$dim[[1]]$vals
# should be a vector of 507 longitudes but:
dim(xx)

# works anyway
xrange <- c(min(xx), max(xx))
xrange

yy <- nc$var[[zvar]]$dim[[2]]$vals
# should be a vector of 356 latitudes but:
dim(yy)

# this fails
yrange <- c(min(yy), max(yy))
yrange

# and that is the only reason why raster failed.

#OK
yy[,,1][1:10,350:355]

#bad values here (and note the NaN):
yy[,,120][1:10,350:355]

# this would have worked
yrange <- c(min(yy[,,1]), max(yy[,,1]))
yrange

close.ncdf(nc)

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.

??? [[alternative HTML version deleted]]

--
advait godbole

--
advait godbole

--
advait godbole