Как фильтровать растровые данные HadISST в диапазоне месяцев в R?

Я пытаюсь изучить корреляцию температуры поверхности моря (ТПМ) с активностью тропических циклонов в определенном диапазоне месяцев. Я использую данные из Центра Хэдли (в формате NetCDF) с использованием функция get_anual_ssts() из пакета hadsstR.

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

Что мне нужно, так это иметь дополнительный параметр, который позволяет мне фильтровать по месяцам ураганной активности вместо расчета среднего ТПМ за весь год.

Например, для юго-западной части Тихого океана я могу назвать get_annual_ssts(hadsst_raster, 12:04, 1966:2007), так как декабрь-апрель — это месяцы активности ураганов. Очень важно установить диапазон месяцев, который включает два разных года (может быть, указать начальный месяц и длину диапазона, чтобы упростить структуру mean_brick, сохранив среднее значение для начального года?).

Глядя на документацию chron, кажется невозможным назначить подмножество мм-гг или что-то подобное. Как это лучше всего сделать?

Вот как выглядят входные растровые данные (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)

r r-raster chron
person Alfredo Hernández    schedule 04.05.2017    source источник
comment
Приведенные выше растровые данные — это то, что вы возвращаете с помощью return(mean_brick)?   -  person Val    schedule 05.05.2017
comment
Это ввод. Я сделаю различие, извините.   -  person Alfredo Hernández    schedule 05.05.2017
comment
Может быть, я вас неправильно понимаю, но мне кажется, что у вас есть информация о месяцах как в names, так и в Date во входном блоке. Вы пробовали индексировать его на основе этого?   -  person Val    schedule 05.05.2017
comment
@Вал, да, это правда. Я беспокоился о фильтрации между месяцами в разные годы. Но, подумав об этом однажды вечером, я понял, что, возможно, вместо этого можно отфильтровать те месяцы, которые мне не нужны. Я опубликую ответ, если придумаю разумный метод.   -  person Alfredo Hernández    schedule 05.05.2017
comment
Правильно ли я предполагаю, что вы хотите отфильтровать определенные месяцы за все годы? Скажем, все слои с апреля за каждый год?   -  person Val    schedule 05.05.2017
comment
@Вал, это правильно. Но учитывая, что мне нужно сравнить данные с активностью ураганов, было бы идеально вычислить среднее значение из того же диапазона активности за месяц. В примере юго-западной части Тихого океана получение среднего значения между январь: апрель и декабрь каждого года приведет к систематической ошибке по сравнению со средним значением декабрь (год_i): апрель (год_i+1).   -  person Alfredo Hernández    schedule 05.05.2017
comment
Хорошо я понял. Таким образом, это не только фильтрация интересующих месяцев, но и правильное их усреднение.   -  person Val    schedule 05.05.2017
comment
@Вал, да. Если бы ураганы были хорошими и происходили в течение одного года, все было бы проще :P   -  person Alfredo Hernández    schedule 05.05.2017

Ответы (1)


arrow_upward
1
arrow_downward

Хорошо, у меня есть кое-что. Возможно, вы используете его для изменения своей функции:

## 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
}
person Val    schedule 05.05.2017
comment
Я предполагаю, что dts будут моими исходными растровыми данными (hadsst_raster), верно? - person Alfredo Hernández; 05.05.2017
comment
Извините, что получится из моего промежуточного продукта... Поправлю. Это должны быть имена слоев nms - person Val; 05.05.2017
comment
sel <-, кажется, переписывается для каждой итерации (комментируя break). Поправьте меня, если я ошибаюсь, но определение sel <- raster() перед циклом, а затем вставка в него (sel <- stack(sel, ... )) поможет, верно? - person Alfredo Hernández; 05.05.2017
comment
Возможно, я должен был быть более ясным: sel — это просто текущий выбор, учитывая год, в котором вы находитесь. Это был просто пример (так сказать, доказательство концепции), отсюда и утверждение break. Вам нужно будет реализовать индексацию в вашей функции. Я посмотрю, есть ли у меня какие-то указатели. Вы хотите выбрать эти растры и выполнить простое среднее? - person Val; 05.05.2017
comment
Итак, first.range означает ноябрь и декабрь предыдущего года и second.range с января по апрель следующего года? - person Val; 05.05.2017
comment
Давайте продолжим обсуждение в чате. - person Alfredo Hernández; 05.05.2017