Оценка климатологии на данных NetCDF с использованием R. Нужен лучший метод - PullRequest
0 голосов
/ 22 февраля 2019

Я работаю над данными о суточной температуре поверхности моря (SST) за 31 год.Данные представлены в формате NetCDF с размерами 28 (длинный) x 40 (широта) x 11686 (дни).Я должен рассчитать среднемесячное климатологическое значение (например, среднее значение за все январь 31 года и т. Д.).Используя библиотеки ncdf4 и chron, я смог получить его в виде массива.

ncin <- nc_open('sstfile.nc')
sst_array <- ncvar_get(ncin, 'sst')

Поскольку переменная времени отделена от данных SST, мне пришлось использовать ее в цикле над массивом.

is.leapyear <- function(year){
return(((year %% 4 == 0) & (year %% 100 != 0)) | (year %% 400 == 0))
}

dateseq <- seq(as.Date("1987-01-01"), as.Date("2018-12-31"), by=1)

Используя библиотеку растров, я конвертирую в растры, а затем выполняю вычисления.

for ( i in seq(11686)) {
dtft <- strsplit(as.character(as.Date(dateseq[i])), split = '-')
y <-  as.integer(dtft[[1]][1])
m <-  as.integer(dtft[[1]][2])
d <-  as.integer(dtft[[1]][3])
while (m == 1){
assign(paste0('r',y,'.',d), raster(matrix(sst_array[1:27, 1:38, i], 
nrow = 27, ncol = 38)))
m = m + 1
}
if (is.leapyear(y) == TRUE) (i = i + 366)
else (i = i + 365)
}

Проблема в том, что он создает слишком много растров и сначала вычисляет среднее значение за месяц, а затем за год.

r87jan <- stack(mget(paste0('r1987.',1:31)))
r87janmean <- calc(r87jan, mean)

Существует ли какая-либо функция / метод, которая может вычислять в течение этого промежутка времени, не создавая столько растров, и вычисления могут оставаться в виде массива или матрицы?Или приведенный выше код может быть лучше рассчитан на все годы сразу?

Ответы [ 2 ]

0 голосов
/ 22 февраля 2019

Если я могу предложить ответ не на R, если у вас есть cdo (установлены операторы климатических данных), вы можете просто сделать это в командной строке linux:

cdo ymonmean sstfile.nc sst_climate.nc 

файл sst_climate.nc будет содержать12 временных шагов, со средним значением за январь, февраль и т. Д. *

Вы можете легко установить cdo, скажем, в ubuntu / mint с

sudo apt-get install cdo 

, а также можетеустановите Ubuntu легко в Windows 10 в эти дни, чтобы легко получить доступ к этим полезным инструментам.Документация доступна здесь https://code.mpimet.mpg.de/projects/cdo/

0 голосов
/ 22 февраля 2019

Вы не предоставляете свои данные, но я думаю, что вы можете сделать что-то вроде этого:

library(raster)
nc <- brick('sstfile.nc')

dates <- getZ(nc)
months <-  as.integer(format(dates, "%m"))

s <- stackApply(nc, months, fun=mean)
...