Substitute pixels values of all layers in a raster stack

Question

I would need help with raster stack. I have a 365-layers raster stack (corresponding to an hydrological year) and I want to substitute the pixel with value 3 with the value of the last previous layer with different value. For example, for a specific pixel I have the first 10 layers that give values: 1,1,2,3,3,3,3,2,1,1. The layers where the value is 3 should be substituted with, in this case, 2, which is the last value, before the "3-window", different from 3. This should be done for all layers (except, obviously, the first one) for each pixel. Any idea? I just convert my raster in a matrix and I work with matrices but this requires a lot of time. Let's suppose that vel_3D_snow_P4 is my matrix (retrieved from raster) where I have 365 rows (days) and more than 400000 columns (pixels), I writed this code:

vel_3D_snow_P4_back=matrix(nrow=nrow(vel_3D_snow_P4),ncol=ncol(vel_3D_snow_P4))
for (i in 1:ncol(vel_3D_snow_P4)){
  y <- which(vel_3D_snow_P4[,i]==3)
  startIndex <- y[!(y-1) %in% y]
  stopIndex <- y[!(y+1) %in% y]

  matrix=matrix(nrow=length(startIndex),ncol=2)
  matrix[,1]=startIndex
  matrix[,2]=stopIndex 

  new_vector_back=vel_3D_snow_P4[,i]
  for (j in 1:nrow(matrix)){
  if (matrix[j,1]==1) next
  new_vector_back[matrix[j,1]:matrix[j,2]]=new_vector_back[matrix[j,1]-1]
  }
vel_3D_snow_P4_back=cbind(vel_3D_snow_P4_back,new_vector_back)
print(c("fine",toString(i)))
}

But, as you can imagine, with numerous pixels it is impossible! This is the reason why I was asking for a solution/idea by maintaining raster format (maybe using calc function?)

Thanks in advance.

Show us what do you have right now. We cannot create a solution for you from scratch as that is actually a job . — M.K

Robert Hijmans Robert Hijmans · Accepted Answer · 2019-02-09T20:27:31

Typically, the first step with problems like this is in R, is to write a function that operates on a vector. Or search for an existing one. My first google query pointed me to zoo::na.locf

library(zoo)
x <- c(1,1,2,3,3,3,3,2,1,1)
x[x==3] <- NA
na.locf(x)
# [1] 1 1 2 2 2 2 2 2 1 1

Then create example raster data

library(raster)
r <- raster(ncol=10, nrow=10)
s <- stack(lapply(c(1,2,3,3,4), function(i) setValues(r, i)))

And combine the two. You can do

A)

x <- reclassify(s, cbind(3, NA))
z <- calc(x, fun=na.locf)

or

B)

f <- function(v) {
    v[v==3] <- NA
    na.locf(v)
}   
zz <- calc(s, f)

Substitute pixels values of all layers in a raster stack

1 Answers