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

Question

我有一个R栅格砖块/堆栈（使用栅格数据包），用于1970年至2015年的45年年度降雨数据。我想计算给定年份的均值，中位数和标准差，例如使用最近5年的2015年，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

要完成此任务，我们可以使用calcraster包中的函数。我们还需要知道如何对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函数示例

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将返回索引。例如，

# 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, 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。