Ограничения с помощью RasterStack по математике и условным выражениям

Question

У меня есть большой RasterStack (114 изображений геотифов), который я успешно геообработал (замаскировал и т. Д. c.) В R, но у меня возникают трудности с его применением простого условного оператора к каждому растру (ко всем слоям растра) имеют одинаковую степень, разрешение и совместно регистрируются). Я хочу установить для всех значений пикселей, которые составляют менее 95% от максимального значения каждого растра, значение NA. Например, если максимальный уровень слоя был 85, то значения пикселей <80,75 = NA. Вот мой код: </p>

#Get max value from each raster layer
r_max <- maxValue(rstack)

#Set all values < 95% of max to NA
rstack[rstack < (r_max * 0.95)] = NA

Когда я запускаю этот код во всем стеке растров, я получаю «Ошибка в значении [j,]: неверное количество измерений». Однако, если я запускаю его на меньшем наборе (14 или около того), он работает точно так же, как и должен. Поскольку я успешно выполнил ряд аналогичных операций (другие условные операторы, маскирование и т. Д. c.) На всем стеке без ошибок, я не уверен, почему он выдает эту ошибку сейчас. Любые идеи?

Я прошу прощения, если это обсуждалось ранее, но я не смог найти такой пост. Если оно существует, пожалуйста, укажите мне в этом направлении.

Спасибо

Answer 1

Пожалуйста, всегда включайте минимальный воспроизводимый пример, такой как:

library(raster)
s <- stack(system.file("external/rlogo.grd", package="raster")) 
cutoff <- maxValue(s) * .95
cutoff
#[1] 242.25 242.25 242.25

Теперь вы можете сделать

s[s < cutoff] = NA
s
#class      : RasterBrick 
#dimensions : 77, 101, 7777, 3  (nrow, ncol, ncell, nlayers)
#resolution : 1, 1  (x, y)
#extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
#crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#source     : memory
#names      : red, green, blue 
#min values : 243,   243,  243 
#max values : 255,   255,  255 

Но есть ошибка, когда RasterStack велик (и должен записать в файл) --- и это то, что вы наткнулись. Мы можем подражать этой ситуации с rasterOptions(todisk=TRUE)

 rasterOptions(todisk=TRUE)
 s[s < cutoff] = NA
 #Error in value[j, ] : incorrect number of dimensions

Я попытаюсь это исправить. Вот обходной путь

s <- stack(system.file("external/rlogo.grd", package="raster")) 
cutoff <- maxValue(s) * .95
x <- sapply(1:nlayers(s), function(i) reclassify(s[[i]], cbind(-Inf, cutoff[i], NA)))
x <- stack(x)