从netCDF更快地读取时间序列？

Question

我有一些大型的netCDF文件，其中包含了地球0.5度分辨率的6小时数据。

每年有360个纬度点、720个经度点和1420个时间点。我有年文件(12 GB ea)和一个具有110年数据(1.3TB)的文件(1.3TB)存储为netCDF-4 (这里是1901年数据、1901.nc、其使用策略和原始的，公开的文件，我开始使用的一个例子)。

据我所理解，从一个netCDF文件读取应该更快，而不是遍历一组按最初提供的年份和变量分列文件。

我想为每个网格点提取一个时间序列，例如从特定的纬度和经度中提取10年或30年的时间序列。然而，我发现这是非常缓慢的。例如，从一个点位置读取10个值需要0.01秒，尽管我可以在0.002秒内从单个时间点读取10000个值的全局切片(维数的顺序是lat、lon、time)：

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

我怀疑这些文件的编写是为了优化空间层的读取，因为( a)变量的顺序是lat，lon，time，b)这将是生成这些文件的气候模型的逻辑顺序，因为全局范围是最常见的可视化。

我试图重新排序变量，以便时间优先：

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

(ncpdq来自NCO (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 

如何在一个点上优化阅读时间序列，而不是在一个时间点上优化一层大面积的时间序列？例如，如果文件以不同的方式编写，比如时间、lat、lon，会更快吗？是否有一种“简单”的方法来转换netCDF-4文件中的维度顺序？

更新

(@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

更新2：

我已经开始在这里开发一个解决方案。这些部分在github.com/ebimodeling/model-drivers/tree/master/met/cruncep中的一组脚本中。

这些脚本仍然需要一些工作和组织--并非所有的脚本都是有用的。但读起来太快了。与上述结果不完全可比，但在一天结束时，我可以立即从一个1.3TB文件(2.5s中的0.5度分辨率)读取一个100年、每小时6小时的时间序列：

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 

(注意:维度的顺序已经改变了，如下所述：到达?)

Answer 1

我认为解决这个问题的方法不是重新排序数据，而是数据块。有关分块netCDF文件的含义的完整讨论，请参阅Unidata首席netCDF开发人员Russ的以下博客文章：

• 分块数据:为什么它重要
• 分块数据:选择形状

其结果是，虽然采用不同的分组策略可以大幅度提高访问速度，但选择正确的策略是非常重要的。

在较小的示例数据集sst.wkmean.1990-present.nc上，我在使用benchmark命令时看到了以下结果：

1)无障碍：

## test replications elapsed relative user.self sys.self user.child sys.child
## 2 spacechunk         1000   0.841    1.000     0.812    0.029          0         0
## 1 timeseries         1000   1.325    1.576     0.944    0.381          0         0

2)天真的块状：

## test replications elapsed relative user.self sys.self user.child sys.child
## 2 spacechunk         1000   0.788    1.000     0.788    0.000          0         0
## 1 timeseries         1000   0.814    1.033     0.814    0.001          0         0

天真的分块只是黑暗中的一次尝试；我这样使用了nccopy实用程序：

$ nccopy -c"lat/100，lon/100，time/100，nbnds/“sst.wkmis.1990-Present.nc chunked.nc

nccopy实用程序的Unidata文档可以找到这里。

我希望我能推荐一个特定的策略来分块你的数据集，但它是高度依赖于数据的。希望上面链接的文章能让您了解如何将您的数据块起来以实现您想要的结果！

更新

下面Marcos Hermida的博客文章展示了不同的分块策略是如何影响特定netCDF文件读取时间序列的速度的。这只应该被用作一个跳伞点。

• Netcdf-4在AR-4三维数据文件上的性能结果

关于通过nccopy进行重分组，显然是挂起的；这个问题似乎与4MB的默认块缓存大小有关。通过将其增加到4GB (或更多)，您可以将大文件的复制时间从超过24小时缩短到11分钟以下！

• 极慢/挂起
• Unidata JIRA故障票系统: NCF-85，改进了nccopy实用程序中块缓存的使用，使其对于重块大文件具有实用价值。

有一点我不确定；在第一个链接中，讨论是关于chunk cache的，但是传递给nccopy的参数-m指定了复制缓冲区中的字节数。nccopy的-m参数控制块缓存的大小。

Answer 2

编辑:最初的问题有一个错误，但也可能有不同的管理费用开始阅读，所以这是公平的多个代表。rbenchmark让这很简单。

这个例子文件有点大，所以我用了一个更小的文件，你能和你的文件做同样的比较吗？

更易访问的示例文件：ftp://ftp.cdc.noaa.gov/Datasets/noaa.oisst.v2/sst.wkmean.1990-present.nc

我得到的时间是时间序列的两倍：

library(ncdf4)

nc <- nc_open("sst.wkmean.1990-present.nc")

library(rbenchmark)
benchmark(timeseries = ncvar_get(nc, "sst", start = c(1, 1, 50), count = c(10, 10, 100)), 
spacechunk = ncvar_get(nc, "sst", start = c(1, 1, 50), count = c(100, 100, 1)),           
replications = 1000)
##        test replications elapsed relative user.self sys.self user.child sys.child
##2 spacechunk         1000    0.47    1.000      0.43     0.03         NA        NA
##1 timeseries         1000    1.04    2.213      0.58     0.47         NA        NA

Answer 3

不确定你是否考虑过cdo来提取要点？

cdo remapnn,lon=x/lat=y in.nc point.nc 

有时CDO内存不足，如果发生这种情况，您可能需要遍历年度文件，然后用

cdo mergetime point_${yyyy}.nc point_series.nc