Переклассифицировать значения в RasterBrick с помощью дополнительного растра (цифровая модель рельефа)

Question

У меня есть RasterBrick , состоящий из ежедневных данных о снежном покрове со значениями 1, 2 и 3 (1 = снег, 2 = нет снега, 3 = затенено облаками).

Пример снежного покрова одного дня:

> 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  

Теперь я хочу интерполировать затененные облаками пиксели (но только в том случае, если в одном слое RasterLayer облачный покров составляет менее 90%, в противном случае исходные значения должны быть сохранены для этого слоя).

Для пространственной интерполяции Я хочу использовать цифровую модель рельефа (та же область исследования и уже в том же разрешении), чтобы извлечь верхний и нижний границы снежной линии для каждого слоя растрового кирпича соответственно . Верхняя снежная линия представляет собой высоту, на которой все безоблачные пиксели классифицируются как снег. Нижняя линия снега определяет высота, ниже которой все безоблачные пиксели также свободны от снега.

> 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)

Для верхних снежных линий Мне нужно минимальное превышение заснеженных пикселей (значение = 1). Теперь все пиксели значения 3 в RasterLayer объекта RasterBrick выше этого минимального уровня должны быть переклассифицированы в значение 1 (предполагается, что оно покрыто снегом).

Для нижней снежной линии , с другой стороны, мне нужно определить максимальную высоту пикселей без снега (значение = 2). Теперь все пиксели значения 3 в RasterLayer объекта RasterBrick выше этого максимального уровня должны быть переклассифицированы в значение 2 (предполагается, что он свободен от снега).

Возможно ли это с помощью R?

Я пытался использовать функцию наложения, но я застрял там.

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

Answer 1

Вот некоторые примеры данных

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)

На графике показаны классы снега (1, 2, 3) со значениями высоты над уровнем моря.

Мы можем использовать маску, но нужно разобраться с пропущенными значениями.

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")

Вы были очень близки, как вы можете

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

Но учтите, что это также перезаписывает высокие ячейки без снега.

Что можно исправить следующим образом:

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

Чтобы выбрать слои с более чем x% ячеек со значением 3 (здесь я использую порог 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