如何在R中过滤Hadisst栅格数据

Question

我正在尝试研究一定范围内的热带环隆活动的海面温度(SST(相关性。我使用的数据来自HADSSTR软件包的get_anual_ssts()函数来自Hadley Center(以NETCDF格式(。

get_annual_ssts <- function(hadsst_raster, years = 1969:2011) {
    mean_rasts <-
        apply(matrix(years), 1, function(x) {
            yearIDx <- which(chron::years(hadsst_raster@z$Date) == x)
            subset_x <- raster::subset(hadsst_raster, yearIDx)
            means <- raster::calc(subset_x, mean, na.rm = TRUE)
            names(means) <- as.character(x)
            return(means)
        })
    mean_brick <- raster::brick(mean_rasts)
    mean_brick <- raster::setZ(mean_brick, as.Date(paste0(years, '-01-01')), 'Date')
    return(mean_brick)
}

我需要的是要有一个附加参数，该参数使我可以通过数月的飓风活动过滤，而不是计算全年平均SST。

例如，对于西南太平洋，我应该能够称呼get_annual_ssts(hadsst_raster, 12:04, 1966:2007)，是12月至4月的飓风活动的月份。设置包括两个不同年份的几个月的范围至关重要(也许是指出初始月份和缓解mean_brick结构的范围长度，从而节省了第一年的平均值？(。

查看chron的文档，似乎不可能分配MM-YY或类似内容的子集。完成此操作的最佳方法是什么？

这是输入栅格数据(hadsst_raster(的样子，以供参考：

class       : RasterBrick 
dimensions  : 180, 360, 64800, 1766  (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 : ~/Downloads/Hadley/HadISST_sst.nc 
names       : X1870.01.16, X1870.02.14, X1870.03.16, X1870.04.15, X1870.05.16, X1870.06.16, X1870.07.16, X1870.08.16, X1870.09.16, X1870.10.16, X1870.11.16, X1870.12.16, X1871.01.16, X1871.02.15, X1871.03.16, ... 
Date        : 1870-01-16, 2017-02-16 (min, max)
varname     : sst 

以及输出(get_annual_ssts(hadsst_raster, 1966:2007)(的样子：

class       : RasterBrick 
dimensions  : 180, 360, 64800, 42  (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 : in memory
names       :      X1966,      X1967,      X1968,      X1969,      X1970,      X1971,      X1972,      X1973,      X1974,      X1975,      X1976,      X1977,      X1978,      X1979,      X1980, ... 
min values  :  -916.8167,  -916.8167,  -916.8167,  -916.8167,  -916.8167,  -916.8167,  -916.8167,  -916.8167,  -916.8167,  -916.8167,  -916.8167,  -916.8167,  -916.8167, -1000.0000, -1000.0000, ... 
max values  :   29.94996,   29.66276,   29.70941,   30.22522,   29.61913,   29.43723,   29.65050,   29.73929,   29.59117,   29.48381,   29.36425,   29.72932,   29.70908,   29.84216,   29.84868, ... 
Date        : 1966-01-01, 2007-01-01 (min, max)

Answer 1

好吧，我有一些东西。也许您使用它来修改您的功能：

## Generate your layer names (used for indexing later)
nms <- expand.grid(paste0('X',1969:2011),c("01","02","03","04","05","06","07","08","09","10","11","12"),'16')
nms <- apply(nms,1,function(x) paste0(x,collapse = '.'))
nms <- sort(nms)
## Generating fake raster brick
r <- raster()
r[] <- runif(ncell(r))
rst <- lapply(1:length(nms),function(x) r)
rst <- do.call(brick,rst)
names(rst) <- nms

现在，您可以用层名称为砖块索引。在飓风季节中循环(从1 -1年开始(：

for (ix in 1970:2011){

  sel <- rst[[c(grep(paste0(ix-1,'.12'),nms),sapply(paste0(0,1:4),function(x) grep(paste0(ix,'.',x),nms)))]]

  break ## in case you don't want to go through all iterations
  }

对于第一次迭代，我将获得此输出：

> sel
class       : RasterStack 
dimensions  : 180, 360, 64800, 5  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
names       :  X1969.12.16,  X1970.01.16,  X1970.02.16,  X1970.03.16,  X1970.04.16 
min values  : 5.988637e-06, 5.988637e-06, 5.988637e-06, 5.988637e-06, 5.988637e-06 
max values  :    0.9999771,    0.9999771,    0.9999771,    0.9999771,    0.9999771 

让我知道这是否有帮助。

---

编辑：

所以也许是一个更适用的示例：

(该功能假设您的输入砖x的图层名称具有格式Xyyyy.mm.dd(

hadSSTmean <- function(x, years, first.range = 11:12, second.range = 1:4){
  nms <- names(x)
  mts <- c("01","02","03","04","05","06","07","08","09","10","11","12")
  xMeans <- vector(length = length(years)-1,mode='list')
  for (ix in 2:length(years){
    xMeans[[ix-1]] <- mean(x[[c(sapply(first.range,function(x) grep(paste0(years[ix-1],'.',mts[x]),nms)),sapply(1:4,function(x) grep(paste0(years[ix],'.',mts[x]),nms)))]])
  }
  return(do.call(brick,xMeans))
  # you could also return the list instead of a single brick
}