计算堆栈栅格中不同时间步长的均值、中值和标准差

Question

我在 R 中有一个栅格 brick/stack（使用栅格包）从 1970 年到 2015 年 45 年的年降雨量数据。我想计算给定年份的平均值、中值和标准偏差，例如2015年使用最近5年、10年、15年、20年、30年。我想从 2000 年到 2015 年执行此操作，其中每年使用堆叠数据重复此过程并保存给定年份的新派生栅格。这是示例代码。任何帮助是极大的赞赏。

raster <- raster(ncol=10, nrow=10)
raster_brick <- brick( sapply(1:45, function(i) setValues(r, rnorm(ncell(r), i, 3))))
plot(raster_brick)
str(raster_brick)

Answer 1

要完成此任务，我们可以使用栅格包中的 calc 函数。我们还需要知道如何对 RasterBrick 对象进行子集化。

资料准备

library(raster)
set.seed(123)
r <- raster(ncol=10, nrow=10)
r_brick <- brick(sapply(1:45, function(i) setValues(r, rnorm(ncell(r), i, 3))))

计算函数示例

calc 可以对RasterBrick 对象中的所有图层应用一个函数。最终结果是栅格图层。

# Calculate mean
r_mean <- calc(r_brick, mean)
# Calculate median
r_median <- calc(r_brick, median)
# Calculate sd
r_sd <- calc(r_brick, sd)

注意 r_mean、r_median 和 r_sd 都是 RasterLayer。

RasterBrick 子集示例

我们可以使用索引对图层进行子集化。例如，

r_sub <- r_brick[[1:3]]

r_sub是r_brick

的前三层

进行分析的函数

知道calc和子集的技巧，我们可以设计一个函数来进行分析。

首先要做的是创建一个向量作为年份和索引的参考。

# Create the index
ind <- 1:45
names(ind) <- 1971:2015

调用年份编号 ind 将 return 索引。例如，

# Get the index of 2015
ind[as.character(2015)]
#2015 
#  45

现在设计函数，它有五个参数

end_year: 分析结束年份

n_year：最后n年的结束年份

FUN：一个函数，例如mean、median和sd

index: 年份指数 (ind)

ras_brick：RasterBrick 与

一起工作

# Define the function

raster_stat <- function(end_year, n_year, FUN, index, ras_brick){
  
  # Subset the raster
  index_temp <- index[as.character((end_year - n_year + 1):end_year)]
  ras_brick_temp <- ras_brick[[index_temp]]
  
  # Calculate the statistics
  ras_result <- calc(ras_brick_temp, FUN)
  
  # Set the name
  names(ras_result) <- paste("Y", end_year, n_year, substitute(FUN), sep = "_")
  
  return(ras_result)
}

现在我们可以测试功能了。

raster_stat(2015, 5, FUN = sd, index = ind, ras_brick = r_brick)
#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       : Y_2015_5_sd 
#values      : 12.05333, 14.61298  (min, max)

请注意 raster_stat 函数的结果名称为 Y_2015_5_sd。这有助于确定应用了哪些 end_year、n_year 和 FUN。

应用函数

我们可以使用 for 循环将 raster_stat 应用于所有 end_year 和 n_year。这是计算 mean.

的示例

# Set the range of end_year and n_year
end_year_vec <- 2000:2015
n_year_vec <- c(5, 10 , 15, 20, 30)

# Create an empty list to store result
r_mean_list <- list()

for (i in end_year_vec){
  for(j in n_year_vec){
      result_temp <- raster_stat(end_year = i, n_year = j, FUN = mean, 
                                 index = ind, ras_brick = r_brick)
      # Add the raster layer to the result_list
      r_mean_list[[names(result_temp)]] <- result_temp
  }
}

所有结果都以唯一的名称存储在 r_mean_list 中。我们可以对 median 和 sd.

使用相同的方法

计算堆栈栅格中不同时间步长的均值、中值和标准差

Calculating mean, median and standard deviation in stack raster for different time steps

stack

r

raster

mean

资料准备

计算函数示例

RasterBrick 子集示例

进行分析的函数

应用函数