Mam kilka dużych plików netCDF, które zawierają 6 danych godzinowych dla Ziemi z rozdzielczością 0,5 stopnia.

Istnieją 360 punktów szerokości geograficznej, 720 punktów długości geograficznej i 1420 punktów czasowych rocznie. Mam zarówno pliki roczne (12 GB ea), jak i jeden plik ze 110-letnimi danymi (1,3 TB) przechowywanymi jako netCDF-4 (oto przykład danych z 1901 r.,1901.nc, jegopolityka użytkowaniaioryginalne, publiczne pliki, z którymi zacząłem).

Z tego, co zrozumiałem, powinno być szybsze odczytywanie z jednego pliku netCDF, niż zapętlanie zestawu plikówoddzielone rokiem i zmienną pierwotnie dostarczone.

Chcę wyodrębnić szereg czasowy dla każdego punktu siatki, np. 10 lub 30 lat od określonej szerokości i długości geograficznej. Jednak uważam to za bardzo powolne. Na przykład odczytanie 10 wartości w czasie z lokalizacji punktu zajmuje mi 0,01 sekundy, chociaż mogę odczytać globalny wycinek 10000 wartości z pojedynczego punktu czasowego w 0,002 sekundy (kolejność wymiarów to lat, lon, czas):

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

Podejrzewam, że te pliki zostały napisane w celu zoptymalizowania odczytu warstw przestrzennych, ponieważ a) kolejność zmiennych to lat, lon, czas, b), która byłaby logicznym porządkiem dla modeli klimatycznych, które wygenerowały te pliki c) ponieważ globalne zakresy są najczęstszą wizualizacją.

Próbowałem zmienić kolejność zmiennych, aby najpierw nadszedł czas:

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

(ncpdq jest zOprogramowanie NCO (operatorzy 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 

Jak mogę zoptymalizować odczyt szeregu czasowego w punkcie, a nie warstwy dużych obszarów w jednym punkcie czasowym? Na przykład, czy byłoby szybciej, gdyby pliki były zapisywane inaczej, takie jak czas, lat, lon? Czy istnieje „łatwy” sposób na przekształcenie kolejności wymiarów w pliku netCDF-4?

Aktualizacja

(benchmarki wymagane przez @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

Aktualizacja 2:

Zacząłem tutaj opracowywać rozwiązanie. Bity i elementy są w zestawie skryptówgithub.com/ebimodeling/model-drivers/tree/master/met/cruncep

Skrypty nadal wymagają trochę pracy i organizacji - nie wszystkie skrypty są przydatne. Ale odczyty są błyskawiczne. Nie do końca porównywalne z powyższymi wynikami, ale pod koniec dnia mogę odczytać 100-letni sześciogodzinny szereg czasowy z pliku 1,3 TB (rozdzielczość 0,5 stopnia w 2,5 s) natychmiast:

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 

(Uwaga: Kolejność wymiarów zmieniła się, jak opisano tutaj:Jak mogę określić kolejność wymiarów podczas używania ncdf4 :: ncvar_get?)