0

我正在尝试研究海表温度 (SST) 与特定月份范围内的热带气旋活动的相关性。我使用的数据来自Hadley Center（NetCDF 格式），使用hadsstR包中的get_anual_ssts()函数。

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)
4

1 回答 1

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

现在您可以使用图层名称来索引砖块。遍历飓风季节（从 Year1 -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
}
于 2017-05-05T13:44:14.037 回答