CODEWITHSUNDEEP
rrasterr-raster
Simon
Simon

Reputation: 1111

Compute raster statistics on sequential list of raster layers

Suppose I have a rasterstack consisting of 5 layers:

library(raster)
r <- raster(nrows=10,ncols=10)
r[] <- rnorm(100)
stk <- stack(r,r*2,r+10,r-7, r*6)

I would like to compute a statistic on a sequential set of the layers based on a different list of numbers.

my.seq<-as.numeric(c("2", "1", "2"))

How can I compute a statistic (e.g. mean) for:

Layers 1:2 in s i.e. my.seq[1]

Layer 3 in s i.e. my.seq[2]

Layers 4:5 in s i.e. my.seq[3]

I feel this would be best achieved using a loop but I cannot work out how this would work.

Any help would be gratefully received.

EDIT

This is what I'm trying to do IRL. I thought I would include it here to help give the question greater context.

I'm trying to compute monthly means from daily temperature data. The rasterstack I'm working with has hundreds of layers, one for every day of the month in sequence. So the layers 1:31 are, when taken together, January 1970. In the above example, my.seq is a list containing the number of days in each month. So I'm trying compute the mean for days 1:31 (January) then 32:60 (February) etc.

Upvotes: 1

Views: 1202

Answers (3)

Robert Hijmans
Robert Hijmans

Reputation: 47036

You can use stackApply for that

Example data

library(raster)
r <- raster(nrows=10, ncols=10, vals=rnorm(100))
stk <- stack(r,r*2,r+10,r-7, r*6)
my.seq <- c(2, 1, 2)

Get the indices and use stackApply

i <- rep(1:length(my.seq), my.seq)
x <- stackApply(stk, i, fun='mean')

Upvotes: 2

Hack-R
Hack-R

Reputation: 23216

This function uses my.seq to extract the layers of interest (including potentially repeated layers, as in the example).

You can put them into a new stack or compute/combine stats on the layer-level in the part of the code below that extracts them.

new_stack <- function(x){
  tmp <- stack()
  for(i in x){
    new <- stk[[i]]
    tmp <- stack(tmp, new)  
  }
  return(tmp)
}

tmp2 <- new_stack(my.seq)

Now, I wasn't 100% sure what precisely how you wanted the mean computed, so I did it several ways below:

stack_mean <- function(x){
  tmp         <- stack()
  layer_means <- numeric()

  for(i in x){
    new <- stk[[i]]

    layer_mean <- mean(stk[[i]]@data@values) 
    print(paste("In layer ", i, "the layer mean is ", layer_mean))

    layer_means <- c(layer_means, layer_mean)

    tmp <- stack(tmp, new)  

  }
  print(paste("Mean of means:", mean(layer_means)))
  return(tmp)
}

stack_mean(my.seq)

This computes several means. First, the means of the layers (which I know is not what you want - I'm just showing how for posterity):

[1] "In layer  2 the layer mean is  -0.0981706786020096"
[1] "In layer  1 the layer mean is  -0.0490853393010048"
[1] "In layer  2 the layer mean is  -0.0981706786020096"

Then the mean of the 3 layers together: 

[1] "Mean of means: -0.081808898835008"

Finally, you have the mean of the stack object, which is apparently a thing though I don't know much about this and have no idea what all this really means:

mean(tmp2) # equivalently: mean(stack_means(my.seq))

class       : RasterLayer 
dimensions  : 10, 10, 100  (nrow, ncol, ncell)
resolution  : 36, 18  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 
data source : in memory
names       : layer 
values      : -5.742679, 4.206359  (min, max)

Upvotes: 1

Jeffrey Evans
Jeffrey Evans

Reputation: 2397

I do not understand your example of my.seq in regard to what you explain in wanting different raster subsets. Here is a worked example where I create a list object with your defined subsets. All you need to do is pass the desired numeric index of subset of rasters to a double bracket index of the stack object. The list object has to be subset using a double bracket as well so, it looks a mess s[[idx[[i]]]]. However, if you break it down, for the first raster subset, it is merely: raster::calc(s[[1:2]], mean)

library(raster)
file <- system.file("external/test.grd", package="raster")
s <- stack(file, file, file)
s <- addLayer(s, raster(file)/2, raster(file)*2)

idx <- list( c(1:2), c(3), c(4:5) )
s.mean <- stack()
  for(i in 1:length(idx)) {
    s.mean <- addLayer(s.mean, calc(s[[idx[[i]]]], mean) )
  }

s.mean
plot(s.mean)

In relation to your expanded question, you can use the date class in R to create indices and even more conveniently, the rts package allows for these types of time-series summaries.

Here, let's create a stack with 365 layers representing days of the year. 

library(raster)
f <- raster(nrows=50, ncols=50, xmn=0, xmx=10)
s <- stack()
  for(i in 1:365) {
    x <- f
    x[] <- runif(ncell(x),0,10)
    s <- addLayer(s, x)
}

Now, we can create the corresponding dates vector. 

( d <- seq(as.Date("1970/1/1"), as.Date("1970/12/31"), "days") )

The dates vector can be queried to provide a raster index. Say we want the mean for December, we can use which to produce the index.

( dec.idx <- which( months(d) == "December") )
( dec.mean <- calc(s[[dec.idx]], mean) )

You can easily create a list containing the raster stack indices for each month which would be equivalent to what you describe in your my.seq object.

months.idx <- list()   
  for(m in unique( months(d) ) ) {
    months.idx[[m]] <- which( months(d) == m)  
  }  
months.idx

However, this is built in functionality to perform these type of temporal summaries in the rts package so, we can shortcut any for loop by coercing the stack to an rts object and then using one of the apply function. 

library(rts)
s.date <- rts::rts(s, d)
( s.month <- apply.monthly(s.date, mean) )

Upvotes: 1

Related Questions