Tengo algunos archivos grandes de netCDF que contienen 6 datos por hora para la Tierra a una resolución de 0,5 grados.

Hay 360 puntos de latitud, 720 puntos de longitud y 1420 puntos por año. Tengo dos archivos anuales (12 GB ea) y un archivo con 110 años de datos (1.3 TB) almacenados como netCDF-4 (aquí hay un ejemplo de los datos de 1901,1901.nc, espolítica de uso, y elArchivos originales, públicos con los que comencé).

Por lo que entendí, debería ser más rápido leer desde un archivo netCDF en lugar de recorrer el conjunto de archivosSeparado por año y variable originalmente proporcionada..

Quiero extraer una serie de tiempo para cada punto de la cuadrícula, por ejemplo, 10 o 30 años a partir de una latitud y longitud específicas. Sin embargo, estoy encontrando que esto es muy lento. Como ejemplo, me toma 0.01 segundos leer 10 valores a lo largo del tiempo desde una ubicación de punto, aunque puedo leer en una porción global de 10000 valores desde un solo punto de tiempo en 0,002 segundos (el orden de las dimensiones es lat, lon, hora):

## 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 

Sospecho que estos archivos se han escrito para optimizar la lectura de las capas espaciales porque a) el orden de las variables es lat, lon, time, b) ese sería el orden lógico para los modelos climáticos que generaron estos archivos yc) porque las extensiones globales Son la visualización más común.

He intentado reordenar las variables para que el tiempo llegue primero:

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

(ncpdq es de laSoftware 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 

¿Cómo puedo optimizar las series de tiempo de lectura en un punto en lugar de capas de áreas grandes en un punto de tiempo? Por ejemplo, ¿sería más rápido si los archivos se escribieran de manera diferente, como time, lat, lon? ¿Existe una forma "fácil" de transformar el orden de las dimensiones en un archivo netCDF-4?

Actualizar

(puntos de referencia 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

Actualización 2:

He comenzado a desarrollar una solución aquí. Los bits y piezas están en un conjunto de scripts engithub.com/ebimodeling/model-drivers/tree/master/met/cruncep

Los scripts aún necesitan trabajo y organización, no todos los scripts son útiles. Pero las lecturas son muy rápidas. No es exactamente comparable a los resultados anteriores, pero al final del día, puedo leer una serie de tiempo de 100 años y seis horas de un archivo de 1.3TB (resolución de 0.5 grados en 2.5s) al instante:

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: El orden de las dimensiones ha cambiado, como se describe aquí:¿Cómo puedo especificar el orden de las dimensiones cuando uso ncdf4 :: ncvar_get?)