gridded time series analysis
On Wed, Dec 22, 2010 at 12:39 PM, Advait Godbole
<advaitgodbole at gmail.com> wrote:
Dear Members,
I was able to find a solution after some deliberation and posting it here
for the benefit of forum members.
To recap, my data was in the form of 2 netCDF files, each of dimension 356 x
507 x 120 (lat,lon,time). I wanted to perform regressions between the two
variables on each grid point *over time*. I used "brick" to read in the
netcdfs as raster layers and stack along with the "*plyr"* package to
perform these regressions on the grid.
# create raster object from the netcdf
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, test method on subset which is 4 x 4 x 120
lim <- extent(rcmraster,180,185,300,305)
rcmraster_crop <- crop(rcmraster,lim)
emitraster_crop <- crop(emitraster,lim)
# do some plotting to check
?par(mfrow=c(1,2))
?plot(rcmraster_crop[[1]])
?plot(emitraster_crop[[1]])
# test regressions for cropped subset
s_crop <- stack(emitraster_crop,rcmraster_crop)
# use of plyr (provides a more elegant alternative to lapply)
# the data frame datCoef contains the slope and intercept for each point in
the grid (denoted by X1 and X2)
datCoef = adply( as.array(s_crop) , 1:2 ,
? ? ? ? ? ? ?function(x ) {
? ? ? ? ? ? ? ? ? ? ? ? ? ?resp = x[1:120]
? ? ? ? ? ? ? ? ? ? ? ? ? ?expl = x[121:240]
? ? ? ? ? ? ? ? ? ? ? ? ? ?mdl = lm(resp ~ expl)
? ? ? ? ? ? ? ? ? ? ? ? ? ?res = c(mdl$coef[1] , mdl$coef[2])
? ? ? ? ? ? ? ? ? ? ? ? ? ?return(res)
? ? ? ? ? ? ? ? ? ? ? ? ?}
? ? ? ? ? ? ? ? ? ? ? ? ?)
I'd do this in two steps:
models <- alply( as.array(s_crop) , 1:2, function(x) {
resp = x[1:120]
expl = x[121:240]
lm(resp ~ expl)
}
Then you can access coefficients with
ldply(models, coef)
and predictions with
laply(models, predict)
but I think you might need to manipulate the output of predict to get
the matrix shape that you are expecting.
Hadley
Assistant Professor / Dobelman Family Junior Chair Department of Statistics / Rice University http://had.co.nz/