regrid данные netcdf в R для интерполяции

Итак, у меня есть некоторые переменные из файла .nc, которые находятся в массивах 4D (x, y, z, t). Дело в том, что координаты z расположены неравномерно, как координаты x и y, то есть z составляет примерно 25 метров, 75 м, 125, 175,..., 500, 600, 700,..., 20000, 21000, 22000. Я пытаюсь линейно интерполировать данные, чтобы получить равномерный интервал 50 м по z. Но функция приблизительного в R работает слишком медленно (мне кажется, массивы слишком большие):

library(ncdf)  
x = get.var.ncdf(nc,'x'); y = get.var.ncdf(nc,'y'); z = get.var.ncdf(nc,'z')  
t = get.var.ncdf(nc,'t')  # time
qc1 = get.var.ncdf(nc,'qc',start=c(1,1,1,1),count=c(-1,-1,-1,-1))  

zlin = seq(z[1],z[length(z)],50)  
qc1_lin = array(0,c(length(x),length(y),length(zlin),length(t)))  
for (i in 1:length(x)) {  
    for (j in 1:length(y)) {  
        for (k in 1:length(t)) {  
            qc1_lin[i,j,,k] = approx(z,qc1[i,j,,k],xout = zlin)  
        }  
    }  
}

Есть ли способ сделать это быстрее? Или кто-то сказал мне изучить данные, чтобы сделать это проще, но я не совсем уверен, что он имеет в виду. Кто-нибудь может мне помочь? Спасибо.


r interpolation netcdf cdo-climate r-grid
person ebonhawkabc    schedule 27.08.2014    source источник
comment
Вам нужно работать со всеми уровнями? Потому что, с моей точки зрения, то, чего вы пытаетесь достичь, не имеет смысла.   -  person    schedule 27.08.2014
comment
Да, мне нужны все уровни. В двух словах, я отслеживаю облака и вершины облаков и сохраняю матрицы данных, которые охватывают 4 км по оси z. Поскольку координаты z более плотно упакованы у земли, это приведет к матрицам разного размера (может потребоваться 60 точек по z для облаков у земли, но только 40 точек для облаков выше по z).   -  person ebonhawkabc    schedule 27.08.2014

Ответы (2)


arrow_upward
3
arrow_downward

Поскольку у меня нет вашего файла ncdf, я использовал в качестве примера набор данных температуры воздуха NOAA:

library(ncdf)
url <- paste("ftp://ftp.cdc.noaa.gov/Datasets/ncep/air.",format(Sys.Date(),"%Y"),".nc",sep="")
download.file(url,destfile="air.nc")
nc <- open.ncdf("air.nc")
x <- get.var.ncdf(nc,'lon')
y <- get.var.ncdf(nc,'lat')
z <- get.var.ncdf(nc,'level')
t <- get.var.ncdf(nc,'time')
qc1 <- get.var.ncdf(nc,'air')

Здесь значение z находится в диапазоне от 1000 до 50, для краткого примера возьмем обычную сетку, разнесенную через каждые 100 уровней (я также ограничу операцию 20 первыми днями набора данных, чтобы пример был относительно небольшим):

zlin <- seq(z[1],z[length(z)],-100)

Используя ваш метод:

qc1_lin <- array(0,dim=c(144,73,10,20))
system.time({
    for (i in 1:length(x)) {  
         for (j in 1:length(y)) {  
             for (k in 1:20) {  
                 # Don't forget that approx outputs a list
                 qc1_lin[i,j,,k] = approx(z,qc1[i,j,,k],xout = zlin)$y
                 }  
             }  
          }
     })
   user  system elapsed 
 26.793   1.196  27.886

Но вы можете использовать apply для выполнения той же операции: аргумент MARGIN также может принимать вектор значений. Здесь мы хотим применить функцию approx к измерениям 1, 2 и 4 (поскольку мы изменяем 3-е измерение):

system.time({
    qc1_lin2 <- apply(qc1[,,,1:20],c(1,2,4),function(X)approx(z,X,xout=zlin)$y)
    })
   user  system elapsed 
 24.413   0.144  24.408

apply к сожалению, выводит новое измерение как первое измерение, поэтому нам нужно изменить результат:

qc1_lin3 <- aperm(qc1_lin2, perm=c(2,3,1,4))

Проверим, что результаты идентичны:

all(qc1_lin3==qc1_lin)
[1] TRUE

Выигрыш во времени относительно небольшой, но, вероятно, он того стоит.

person plannapus    schedule 20.10.2014

arrow_upward
1
arrow_downward

Это не ответ в R, а просто чтобы сказать, что эту задачу можно быстро сделать из командной строки с помощью CDO

 cdo intlevel,`seq -s "," 50 50 22000` in.nc out.nc

команда seq создает список, разделенный запятыми, от 50 до 22000 с интервалом 50 м.

person Adrian Tompkins    schedule 17.09.2019