gridded time series analysis
Dear Advait,
What I dont understand is: 1) How can each slope be the same (values of x1)
Apparently the values in these cells are all the same (not that you
only use the first 8 time steps, not all 120)
Here is an example with random values in which the slopes are all different
b <- brick(ncol=4, nrow=4, nl=120)
b <- setValues(b, matrix(runif(prod(dim(b))), ncol=120))
fun = function(x) { if(is.na(x[1])) { NA } else { lm(x[1:60] ~
x[61:120])$coefficients[2] } }
x <- calc(b, fun)
plot(x)
as.matrix(x)
2) The raster layers are of dimensions 4 x 4 x 120 but?dim(getValues(rcmraster_crop)) [1] ?16 120. Why is that? This must have implications for the way that the regressions for the brick are handled wouldnt it?
No
Or is it specific to the way "getValues" handles things?
Yes. rows are cells (4*4=16), columns are layers.
Also, does anyone have any pointers of how to accomplish the NA testing and
See an example above, assuming that if the first time step is NA, they all are. You could also do things like if (sum(is.na(x)) > 60) depends on your context. Best, Robert
On Mon, Dec 6, 2010 at 8:30 AM, Advait Godbole <advaitgodbole at gmail.com> wrote:
Dear Members,
To recap, I initially had trouble reading in two 3D netCDFs in to R using
the raster "brick" function. Each brick is a 356 x 507 x 120 (lat, lon,
time) raster. With that resolved (thanks Robert!) I have the following code
to perform regressions between the ?two time series (at each grid-cell)
# read the netCDFs
rcmraster <- brick("regcm.cumul.hdd.10k.96_05.fixed.latlonvars.nc", varname
= "hdd")
emitraster <-
brick("flip.percap.vulcan.10k.dp.build.381.RES.EIAmy.96_05.this.latest.nc",varname="emitpercap")
# crop subset of above rasters to test
lim <- extent(rcmraster,180,185,300,305) ?# createss a 4 x 4 x 120 brick
rcmraster_crop <- crop(rcmraster,lim)
emitraster_crop <- crop(emitraster,lim)
#regressions
s_crop <- stack(emitraster_crop,rcmraster_crop)
func <- function(x) { lm(x[1:4] ~ x[5:8])$coefficients[2] } ? #regression
for emit* VS rcm*
x1 <- calc(s_crop, func)
While the above regressions run, I have to ultimately do the same ?for the
?entire domain (356 x 507 cells) and therefore need to check for NAs at all
time steps for each grid cell for both the bricks, and then have lm?run. To
that ?end,, I ran the following diagnostics on the cropped dataset:
dim(rcmraster_crop)
[1] ? 4 ? 4 120
dim(emitraster_crop)
[1] ? 4 ? 4 120
dim(x1)
[1] 4 4
getValues(x1)
?[1] 3.968207 3.968207 3.968207 3.968207 3.968207 3.968207 3.968207 3.968207 3.968207 3.968207 3.968207 [12] 3.968207 3.968207 3.968207 3.968207 3.968207 What I dont understand is: 1) How can each slope be the same (values of x1) 2) The raster layers are of dimensions 4 x 4 x 120 but?dim(getValues(rcmraster_crop)) [1] ?16 120. Why is that? This must have implications for the way that the regressions for the brick are handled wouldnt it? Or is it specific to the way "getValues" handles things? Also, does anyone have any pointers of how to accomplish the NA testing and regressions ? Thanks, Advait On Sun, Dec 5, 2010 at 6:29 PM, Robert J. Hijmans <r.hijmans at gmail.com> wrote:
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:
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
-- advait godbole