计算栅格堆栈中高于特定阈值的连续天数的最大长度 [英] Calculate maximum length of consecutive days above a certain threshold in a raster stack

问题描述

我想计算超过阈值 t 的最大连续天数，给定栅格堆栈 s，如下所示:

I would like to calculate the maximum length of consecutive days above a threshold t given a raster stack s as shown below:

library(raster)

set.seed(112)

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))
s=stack(x1,x2,x3,x4,x5,x6)*56

这是我目前的功能.

fun <- function(x,t){
  y <- rle((x > t)*1)
  z <- y$lengths[y$values==1]
  return(max(z,0))
}

我还按照 cluster {raster} 函数中的建议设置了一个参数 q 用于导出

I have also set a parameter q for export as advised in the cluster {raster} function

q <- 0

我希望将栅格图层作为输出，但会弹出以下错误.

I expect a raster layer as an output but instead the error below pops up.

[1] "cannot use this function"
attr(,"class")
[1] "snow-try-error" "try-error"     
Error in clusterR(s, calc, args = list(fun = fun), export = "q") : 
  cluster error

可能是什么问题?

推荐答案

首先，如果您在示例数据中使用随机值，请同时设置随机种子以使其可重现.

First, if you use random values in your example data, please also set the random seed so it's reproducible.

library(raster)

set.seed(42)

x1 <- raster(nrows=10, ncols=10)

s <- do.call(stack,lapply(1:6,function(x) setValues(x1,runif(ncell(x1)))*56))

至于您的问题，您唯一需要的是一个简单的函数，可以将其传递给 calc 以获得所需的结果:

As to your question, the only thing you need is a simple function that can be passed into calc to obtain the desired results:

cd <- function(x,t){

  y <- rle((x > t)*1)

  z <- y$lengths[y$values==1]

  return(max(z,0))

}

该函数使用 rle 或运行长度编码来计算向量中连续运行的次数.在这种情况下，我正在寻找连续 1 的最大数量，它来自将 TRUE 值(值高于阈值 t)与 1 相乘.

This function uses rle, or run length encoding, to calculate the number of consecutive runs in a vector. In this case I'm looking for the maximum number of consecutive 1s, which come from multiplying TRUE values (value is above threshold t) with 1.

最后，您希望返回值 1 的最大运行，0 是在没有发生的情况下的回退(旁注:1 表示单个非连续发生).

In the end you want to return the maximum run of a value 1, with 0 being a fallback in case there's no occurrence (sidenote: 1 indicates a single, non-consecutive occurence).

最后，cd 可以传递给 calc，在这种情况下使用 40 的阈值:

Finally, cd can be passed into calc, in this case using a threshold of 40:

plot(calc(s,function(x) cd(x,40)))

这篇关于计算栅格堆栈中高于特定阈值的连续天数的最大长度的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！