Search code examples
rrasterr-raster

Looping through a raster brick

I know there are a large number of questions on here about looping through raster bricks, but none of them provide quite the answer/advice I'm looking for.

I have a large (17.2GB, 7901 layers) netcdf file that I've imported into R as a RasterBrick. 

> KK10Brick
class       : RasterBrick 
dimensions  : 2160, 4320, 9331200, 7901  (nrow, ncol, ncell, nlayers)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : D:\LandUse\KK10.nc 
names       : X8000, X7999, X7998, X7997, X7996, X7995, X7994, X7993, X7992, X7991, X7990, X7989, X7988, X7987, X7986, ... 
z-value     : 100, 8000 (min, max)
varname     : land_use

Each layer in the file represents 1 year and I need to create a temporal moving average of each pixel in the brick. Even though the variable seems categorical (land_use) it is actually a % cover.

I want to create a 30 year moving average, with a 10 year sliding window. e.g. the first window would produce a raster of the average values from layers 1:30, the next window produces another raster of the average values from layers 11:40...7871:7901.

I was thinking a for loop would probably be the best way to accomplish this, but I'm not sure if I'm headed down the right path here e.g.

for (i in 1:7901){
subsetLayers <- code to subset relevant layers
out <- stackApply(KK10Brick, indices = subsetLayers, fun = "mean", na.rm = TRUE, filename = paste("./Output/", "meanLU_window_", i, ".tif", sep = ""))
rm(out)}

Where I'm stuck is writing the code to produce the sequences for subsetLayers. Any help would be hugely appreciated.

EDIT.

library(raster)
exBrick <- brick(nrow = 180, ncol = 360, nl = 100)
values(exBrick) <- runif(ncell(exBrick) * nlayers(exBrick))
crs(exBrick) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" 
exBrick
Solution

  • This should work on your example data. I'm not sure how well it will scale to your very large netcdf data in terms of speed and RAM usage - please let me know if it works on large data.

    starts = seq(1, nlayers(exBrick)-30, 10)
nout = length(starts)
out = brick(nrow = 180, ncol = 360, nl = nout)
values(out) = NA
crs(out) = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" 

for (i in 1:nout) {
  start = starts[i]
  out[[i]] = mean(exBrick[[start:(start+30)]]) 
}

    In case RAM usage is the limiting factor from allocating a large brick to store the results, we can save RAM at the cost of some speed by saving each result layer to disk, one raster at a time:

    for (i in starts) {
  out = mean(exBrick[[i:(i+30)]]) 
  writeRaster(out, filename=paste0("out",i,".grd"), overwrite=TRUE)
}