Вычисление количества трех последовательных значений выше порога (в растровом стеке)

Мне нужно вычислить количество трех дней подряд, когда значение каждого пикселя в растровом стеке (x) превышает заданный порог (определяемый другим растром y). Я попытался использовать rle для этой цели с calc следующим образом после объединения x и y в новый растр a:

library(raster)    
fn<-function(a) with(rle(a), sum(lengths>=3 & values>a[[nlayers(a)]]))
calc(b,fn)

Однако я получаю сообщение об ошибке:

Ошибка в .calcTest(x[1:5], fun, na.rm, forcefun, forceapply):
не может использовать эту функцию

Воспроизводимый образец:

x1 <- raster(nrows=10, ncols=10)
 x2=x3=x4=x5=x6=x1
 x1[]= runif(ncell(x1))
 x2[]= runif(ncell(x1))
 x3[]= runif(ncell(x1))
 x4[]= runif(ncell(x1))
 x5[]= runif(ncell(x1))
 x6[]= runif(ncell(x1))
 x=stack(x1,x2,x3,x4,x5,x6)
 y=x1
 y[]= runif(ncell(x1))
 a<-stack(x,y)

Может кто-нибудь, пожалуйста, помогите.


r stack raster r-raster run-length-encoding
person rar    schedule 07.11.2016    source источник
comment
Пожалуйста, приведите воспроизводимый пример (вам нужно dput() некоторые примеры данных). Наведите указатель мыши на тег r, чтобы получить дополнительную информацию.   -  person Hack-R    schedule 07.11.2016
comment
@ Hack-R Я соответствующим образом отредактировал вопрос.   -  person rar    schedule 07.11.2016

Ответы (2)


arrow_upward
0
arrow_downward

Вы можете попробовать это:

 x1 <- raster(nrows=10, ncols=10)
 x2=x3=x4=x5=x6=x1
 x1[]= runif(ncell(x1))
 x2[]= runif(ncell(x1))
 x3[]= runif(ncell(x1))
 x4[]= runif(ncell(x1))
 x5[]= runif(ncell(x1))
 x6[]= runif(ncell(x1))
 x=stack(x1,x2,x3,x4,x5,x6)*4
 y=x1
 y[]= runif(ncell(x1))*2
 a<-stack(x,y)


 library(raster)

fn<-function(x) {
  seq <- rle(as.numeric(x)>as.numeric(x)[[nlayers(a)]])
  n = length(seq$lengths > 3 & seq$values == TRUE)
  return(n)
}
calc(a,fn)

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 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 1, 7  (min, max)

(обратите внимание, что я изменил пример набора данных, чтобы получить хорошие последовательности)

ХТХ!

person lbusett    schedule 07.11.2016
comment
Эй, спасибо! Этот код работает очень хорошо. Но когда я использую эту функцию для своих исходных данных, я просто получаю 1 в качестве значения для всех пикселей :( - person rar; 08.11.2016
comment
не знаю .... в тестовом наборе данных я думаю, что код делал то, что должен был делать ... У вас случайно нет NA во входных данных? - person lbusett; 08.11.2016

arrow_upward
0
arrow_downward

Основываясь на ответе Лоренцо, я изменил код, и это именно то, что я искал.

# create data
x1 <- raster(nrows=2, ncols=2)
x2=x3=x4=x5=x6=x1
x1[]= runif(ncell(x1))
x2[]= runif(ncell(x1))
x3[]= runif(ncell(x1))
x4[]= runif(ncell(x1))
x5[]= runif(ncell(x1))
x6[]= runif(ncell(x1))
x=stack(x1,x2,x3,x4,x5,x6)*4
y=x1
y[]= runif(ncell(x1))*2

#Overlay to get greater than values raster (in form of TRUE and FALSE)
a<-overlay(x,y,fun=function(x,y){x>y})
# function for finding total number of 3 consecutive days when condition is TRUE     
fn<-function(a) {
      seq <- rle(a)
      n=sum(seq$lengths[seq$lengths>=3 & seq$values==TRUE])
      return(n)
    }
    calc(a,fn)

Но поскольку мое решение основано на вашем ответе, я принимаю ваш ответ, Лоренцо. Спасибо!

person rar    schedule 08.11.2016