title: Faster reading of time series from netCDF?
body: <p>I have some large netCDF files that contain 6 hourly data for the earth at 0.5 degree resolution. </p>
  
  <p>There are 360 latitude points, 720 longitude points, and 1420 time points per year. I have both yearly files (12 GB ea) and one file with 110 years of data (1.3 TB) stored as netCDF-4 (here is an example of the 1901 data, <a href="http://file-server.igb.illinois.edu/~dlebauer/met/1901.nc" rel="noreferrer">1901.nc</a>, its <a href="http://nacp.ornl.gov/MsTMIP_datapolicy.shtml" rel="noreferrer">use policy</a>, and the <a href="ftp://nacp.ornl.gov/synthesis/2009/frescati/model_driver/cru_ncep/analysis/" rel="noreferrer">original, public files that I started with</a>).</p>
  
  <p>From what I understood, it should be faster to read from one netCDF file rather than looping over the set of files <a href="ftp://nacp.ornl.gov/synthesis/2009/frescati/model_driver/cru_ncep/analysis/" rel="noreferrer">separated by year and variable originally provided</a>.</p>
  
  <p>I want to extract a time series for each grid point, e.g. 10 or 30 years from a specific latitude and longitude. However, I am finding this to be very slow. As an example, it takes me 0.01 seconds to read in 10 values over time from a point location although I can read in a global slice of 10000 values from a single time point in 0.002 second (the order of the dimensions is lat, lon, time):</p>
  
  <pre><code>## a time series of 10 points from one location:
  library(ncdf4)
  met.nc &lt;- nc_open("1901.nc")
  system.time(a &lt;- 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 &lt;- nc_open("1901.nc"))
  system.time(a &lt;- ncvar_get(met.nc, "lwdown", start = c(100,100,1), 
                                               count = c(100,100,1)))
     user  system elapsed 
    0.002   0.000   0.002 
  </code></pre>
  
  <p>I suspect that these files have been written to optimize reading of spatial layers because a) the order of variables is lat, lon, time, b) that would be the logical order for the climate models that generated these files and c) because global extents are the most common visualization.</p>
  
  <p>I have attempted to reorder variables so that time comes first:</p>
  
  <pre><code>ncpdq -a time,lon,lat 1901.nc 1901_time.nc
  </code></pre>
  
  <p>(<a href="http://www.unidata.ucar.edu/software/netcdf/workshops/2008/utilities/Nco.html" rel="noreferrer">ncpdq</a> is from the <a href="http://nco.sourceforge.net/" rel="noreferrer">NCO (netCDF operators) software</a>)</p>
  
  <pre><code>&gt; library(ncdf4)
  
  ## first with the original data set:
  &gt; system.time(met.nc &lt;- nc_open("test/1901.nc"))
     user  system elapsed 
    0.024   0.045  22.334 
  &gt; system.time(a &lt;- 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:
  &gt; system.time(met_time.nc &lt;- nc_open("test/1901_time.nc"))
     user  system elapsed 
    0.025   0.041  16.704 
  &gt; system.time(a &lt;- 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 
  </code></pre>
  
  <p>How can I optimize reading time series at a point rather than layers of large areas at one time point? For example, would it be faster if the files were written differently, such as time, lat, lon? Is there an "easy" way to transform the order of dimensions in a netCDF-4 file?</p>
  
  <h2>Update</h2>
  
  <p>(benchmarks requested by @mdsumner)</p>
  
  <pre><code>library(rbenchmark)
  library(ncdf4)
  nc &lt;- 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
  </code></pre>
  
  <h2>Update 2:</h2>
  
  <p>I have started developing a solution here. The bits and pieces are in a set of scripts in <a href="https://github.com/ebimodeling/model-drivers/blob/master/met/cruncep/README.md" rel="noreferrer">github.com/ebimodeling/model-drivers/tree/master/met/cruncep</a></p>
  
  <p>The scripts still need some work and organization - not all of the scripts are useful. But the reads are lightning quick. Not exactly comparable to the above results, but at the end of the day, I can read a 100 year, six-hourly time series from a 1.3TB file (0.5 degree resolution in 2.5s) instantly:</p>
  
  <pre><code>system.time(ts &lt;- ncvar_get(met.nc, "lwdown", start = c(50, 1, 1), count = c(160000, 1, 1)))
     user  system elapsed 
    0.004   0.000   0.004 
  </code></pre>
  
  <p>(<em>note:</em> The order of dimensions have changed, as described here: <a href="https://stackoverflow.com/q/22944707/199217">How can I specify dimension order when using ncdf4::ncvar_get?</a>)</p>