干湿月 rasterbrick R 出现频率

Question

我有一个 rasterbrick，a 由 1950-2014 年的月度时间序列组成（见下文）。我试图在 a 上完成的任务是计算每个网格点：

negative/positive number of dry/wet months per total months of dry/wet events*100,

哪里

干值在 −1.49 to −1.00 范围内，湿值在 1.49 to 1.00.

范围内

结果输出是一个单一的 rasterLayer，干 = 负百分比，湿 = 正百分比，这样我就可以用蓝色表示湿，红色表示干。

可以找到示例数据HERE

dd=spei03_df
dd[1:2]<-dd[2:1]#swap lat and lon
a=rasterFromXYZ(dd)
dates=seq(as.Date("1950-01-01"), as.Date("2014-12-31"), by="month")
a=setZ(a,dates)

Answer 1

这里有一个适合您的解决方案。它不会生成正百分比和负百分比，而是生成您的值之一的百分比。百分比的工作方式，您仍然可以规划出您想要的内容。

注意：我没有使用您的数据集，而是创建了一个具有相同值范围的示例数据。这些值将或多或少地随机采样（我加入了偏暖*以获得一种颜色的趋势），所以不要期望有漂亮的图案。

library(raster)

## create sample data

#values
dry <- seq(-1,-1.49,by=-0.01)
wet <- dry * -1

#raster brick

r <- lapply(seq_len(length(1950:2014)*12),function(x){

  r <- raster()
  r[] <- sample(c(wet,dry),ncell(r),T,c(rep(0.3,50),rep(0.7,50)))
  r
})

r <- do.call(brick,r)

names(r) <- seq(as.Date("1950-01-01"), as.Date("2014-12-31"), by="month")

#output

> r
class       : RasterBrick 
dimensions  : 180, 360, 64800, 780  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 
data source : /tmp/Rtmp3BTUH1/raster/r_tmp_2017-08-09_120141_2859_73743.grd 
names       : X1950.01.01, X1950.02.01, X1950.03.01, X1950.04.01, X1950.05.01, X1950.06.01, X1950.07.01, X1950.08.01, X1950.09.01, X1950.10.01, X1950.11.01, X1950.12.01, X1951.01.01, X1951.02.01, X1951.03.01, ... 
min values  :       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49, ... 
max values  :        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49, ... 
  -1.49,   -1.49,    -1.49,    -1.49,    -1.49,    -1.49,    -1.49,    -1.49, ... 
max values  :    1.49,    1.49,    1.49,    1.49,    1.49,    1.49,    1.49,    1.49,    1.49,     1.49,     1.49,     1.49,     1.49,     1.49,     1.49, ...

现在计算频率：

rfreq <- calc(r,fun = function(x){(sum(x > 0) / nlayers(r))*100})

像这样，每个像素的值都在 0 到 100 之间，表示从所有干燥月份到所有潮湿月份的范围，平衡值为 50。因此您可以轻松地将其映射为红色和蓝色。如果您仍然想要百分比，您可以添加带有相应标签的图例：

library(rasterVis)

my.at <- c(0,30,50,70,100)

myColorkey <- list(at=my.at,
                   space = 'bottom',
                   labels=list(
                     labels=c('extreme dry','dry','normal','wet','extreme wet'), ## labels
                     at=c(15,40,50,60,85) ## where to print labels
                   ))

levelplot(rfreq, at=my.at,
          colorkey=myColorkey,scales=list(draw=FALSE),margin=FALSE,
          par.settings=RdBuTheme(),
          main = 'Precipitation 1950 - 2014',
          xlab=NULL, ylab=NULL)

HTH

*

气候变化是真实存在的

干湿月 rasterbrick R 出现频率

Frequency of occurrence of dry and wet months rasterbrick R

r

raster

time-series

frequency