[R-sig-Geo] gridded time series analysis

[Hadley Wickham]

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/