Reclassify values in a RasterBrick by the use of an additional Raster (Digital elevation model)

I have a RasterBrick consisting of daily snow cover data with the values 1, 2 and 3 (1= snow, 2= no snow, 3= cloud-obscured).

Example of snow cover of one day:

> snowcover
class       : Large RasterBrick 
dimensions  : 26, 26, 2938  (nrow, ncol, nlayers)
resolution  : 231, 232  (x, y)
extent      : 718990, 724996, 5154964, 5160996  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84       
              +towgs84=0,0,0  

Now I wish to interpolate the cloud-obscured pixels (but only where there are less than 90 % cloud cover in a single RasterLayer, otherwise the original values should be retained for this Layers).

For spatial interpolation I want to use a digital elevation model (same study area and already in same resolution) to extract upper and lower snowline boundaries for each Layer of the RasterBrick respectively. The upper snow line represents the elevation where all cloud-free pixels are classified as snow. The lower snowline identifies the altitude below which all cloud-free pixels are also snow-free.

> dem
class       : RasterLayer 
resolution  : 231, 232  (x, y)
extent      : 718990.2, 724996.2, 5154964, 5160996 (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84       
              +towgs84=0,0,0 
values      : 1503, 2135  (min, max)

For the upper snowlines I need the minimum elevation of the snow-covered pixels (value = 1). Now all pixels of value 3 in a RasterLayer of the RasterBrick above this minimum elevation, should be reclassified as value 1 (assumed to be snow-covered).

For the lower snowline on the other hand I need to identify the maximum elevation of the no-snow-pixels (value = 2). Now all pixels of value 3 in a RasterLayer of the RasterBrick above this maximum elevation should be reclassified as value 2 (assumed to be snow-free).

Is this possible using R?

I tried to make use of the overlay function, but I got stuck there.

# For the upper snowline:
overlay <- overlay(snowcover, dem, fun=function(x,y){ x[y>=minValue(y[x == 1])] <- 1; x})

Here is some example data

library(raster)
dem <- raster(ncol=8, nrow=7, xmn=720145, xmx=721993, ymn=5158211, ymx=5159835, crs='+proj=utm +zone=32 +datum=WGS84')
values(dem) <- ncell(dem):1
snow <- setValues(dem, c(1, 1, rep(1:3, each=18)))
snow[,c(2,5)] <- NA
snow[3] <- 3


plot(snow)
lines(as(dem, 'SpatialPolygons'))
text(dem)

The plot shows the snow classes (1, 2, 3) with the elevation values on top.

We can use mask, but need to deal with the missing values.

msnow <- reclassify(snow, cbind(NA, 0))
# mask to get only the snow elevations
x <- mask(dem, msnow, maskvalue=1, inverse=TRUE)

# minimum elevation of the snow-covered cells
minsnow <- minValue(x)
minsnow 
#[1] 37

# snow elevation = 1
snowy <- reclassify(dem, rbind(c(-Inf, minsnow, NA), c(minsnow, Inf, 1)))
newsnow <- cover(snow, snowy)

s <- stack(dem, snow, newsnow)
names(s) <- c("elevation", "old_snow", "new_snow")

You were very close, as you can do

 r <- overlay(dem, snow, fun=function(e, s){ s[e >= minsnow] <- 1; s})

But note that that also overwrites high cells with no snow.

Which could be fixed like this:

r <- overlay(dem, snow, fun=function(e, s){ s[e >= minsnow & is.na(s)] <- 1; s})

To select layers with more than x% cells with value 3 (here I use a threshold of 34%):

threshold = .34
s <- stack(snow, snow+1, snow+2)
f <- freq(snow)
f 
#     value count
#[1,]     1    14
#[2,]     2    13
#[3,]     3    15
#[4,]    NA    14

nas <- f[is.na(f[,1]), 2]

ss <- subs(s, data.frame(from=3, to=1, subsWithNA=TRUE))
cs <- cellStats(ss, sum)
csf <- cs / (ncell(snow) - nas)
csf
#  layer.1   layer.2   layer.3 
#0.3571429 0.3095238 0.3333333 

i <- which(csf < threshold)
use <- s[[i]]
#use
class       : RasterStack 
dimensions  : 7, 8, 56, 2  (nrow, ncol, ncell, nlayers)
resolution  : 231, 232  (x, y)
extent      : 720145, 721993, 5158211, 5159835  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=32 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
names       : layer.2, layer.3 
min values  :       2,       3 
max values  :       4,       5