Eu tenho alguns arquivos netCDF grandes que contêm 6 dados por hora para a terra em resolução de 0,5 grau.

Existem 360 pontos de latitude, 720 pontos de longitude e 1420 pontos de tempo por ano. Eu tenho dois arquivos anuais (12 GB ea) e um arquivo com 110 anos de dados (1,3 TB) armazenados como netCDF-4 (aqui é um exemplo dos dados de 1901,1901.nc, Estápolítica de uso, e asarquivos públicos originais que eu comecei com).

Pelo que entendi, deve ser mais rápido ler de um arquivo netCDF em vez de fazer um loop sobre o conjunto de arquivosseparados por ano e variável originalmente fornecida.

Desejo extrair uma série temporal para cada ponto da grade, por exemplo, 10 ou 30 anos a partir de uma latitude e longitude específicas. No entanto, estou achando que isso é muito lento. Como exemplo, eu levo 0,01 segundos para ler 10 valores ao longo do tempo a partir de uma localização de ponto, embora eu possa ler em uma fatia global de 10000 valores de um único ponto no tempo em 0,002 segundo (a ordem das dimensões é lat, lon, Tempo):

## a time series of 10 points from one location:
library(ncdf4)
met.nc <- nc_open("1901.nc")
system.time(a <- ncvar_get(met.nc, "lwdown", start = c(100,100,1), 
                                             count = c(1,1,10)))
   user  system elapsed 
  0.001   0.000   0.090 

## close down session

## a global slice of 10k points from one time
library(ncdf4)
system.time(met.nc <- nc_open("1901.nc"))
system.time(a <- ncvar_get(met.nc, "lwdown", start = c(100,100,1), 
                                             count = c(100,100,1)))
   user  system elapsed 
  0.002   0.000   0.002 

Eu suspeito que esses arquivos foram escritos para otimizar a leitura de camadas espaciais porque a) a ordem das variáveis ​​é lat, lon, time, b) que seria a ordem lógica para os modelos climáticos que geraram esses arquivos ec) porque as extensões globais são a visualização mais comum.

Eu tentei reordenar as variáveis ​​para que o tempo chegue primeiro:

ncpdq -a time,lon,lat 1901.nc 1901_time.nc

(ncpdq é doSoftware NCO (operadores netCDF))

> library(ncdf4)

## first with the original data set:
> system.time(met.nc <- nc_open("test/1901.nc"))
   user  system elapsed 
  0.024   0.045  22.334 
> system.time(a <- ncvar_get(met.nc, "lwdown", start = c(100,100,1), count = c(1, 1, 1000))
+ )
   user  system elapsed 
  0.005   0.027  14.958 

## now with the rearranged dimensions:
> system.time(met_time.nc <- nc_open("test/1901_time.nc"))
   user  system elapsed 
  0.025   0.041  16.704 
> system.time(a <- ncvar_get(met_time.nc, "lwdown", start = c(100,100,1), count = c(1, 1, 1000)))
   user  system elapsed 
  0.001   0.019   9.660 

Como posso otimizar a leitura de séries temporais em um ponto em vez de camadas de grandes áreas em um único momento? Por exemplo, seria mais rápido se os arquivos fossem escritos de forma diferente, como time, lat, lon? Existe uma maneira "fácil" de transformar a ordem das dimensões em um arquivo netCDF-4?

Atualizar

(benchmarks solicitados por @mdsumner)

library(rbenchmark)
library(ncdf4)
nc <- nc_open("1901.nc")
benchmark(timeseries = ncvar_get(nc, "lwdown", 
                                 start = c(1, 1, 50), 
                                 count = c(10, 10, 100)), 
          spacechunk = ncvar_get(nc, "lwdown", 
                                  start = c(1, 1, 50), 
                                  count = c(100, 100, 1)),           
          replications = 1000)

        test replications elapsed relative user.self sys.self user.child
2 spacechunk         1000   0.909    1.000     0.843    0.066          0
1 timeseries         1000   2.211    2.432     1.103    1.105          0
  sys.child
2         0
1         0

Atualização 2:

Eu comecei a desenvolver uma solução aqui. Os bits e peças estão em um conjunto de scripts emgithub.com/ebimodeling/model-drivers/tree/master/met/cruncep

Os scripts ainda precisam de algum trabalho e organização - nem todos os scripts são úteis. Mas as leituras são muito rápidas. Não é exatamente comparável aos resultados acima, mas no final do dia, posso ler uma série temporal de 100 anos e seis horas de um arquivo de 1,3 TB (resolução de 0,5 grau em 2,5s) instantaneamente:

system.time(ts <- ncvar_get(met.nc, "lwdown", start = c(50, 1, 1), count = c(160000, 1, 1)))
   user  system elapsed 
  0.004   0.000   0.004 

(Nota: A ordem das dimensões foi alterada, conforme descrito aqui:Como posso especificar a ordem de dimensão ao usar ncdf4 :: ncvar_get?)