从netCDF中快速读取时间序列？

Question

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

每年有360个纬度点，720个经度点和1420个时间点。我有年度文件（12 GB ea）和一个包含110年数据（1.3 TB）的文件存储为netCDF-4（这是1901年数据，1901.nc，其使用策略以及我开始使用的原始公共文件的示例）。

根据我的理解，从一个netCDF文件中读取应该比遍历最初提供的由年份和变量分隔的文件集更快。

我想为每个网格点提取一个时间序列，例如从特定的纬度和经度中提取 10 年或 30 年。但是，我发现这很慢。例如，尽管我可以在 0.002 秒内从单个时间点读取 10000 个值的全局切片（维度的顺序是纬度、lon、时间），但随着时间的推移，我从一个点位置读取 10 个值需要 0.01 秒：

## 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）这将是生成这些文件的气候模型的逻辑顺序，c）因为全球范围是最常见的可视化。

我尝试重新排序变量，以便时间优先:

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 

我如何在一个时间点优化阅读时间序列，而不是在一个时点优化大面积的分层？例如，如果文件写得不同，比如time、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

我已经开始在这里开发一个解决方案。这些片段位于github . com/ebi modeling/model-drivers/tree/master/met/cruncep中的一组脚本中

这些脚本仍然需要一些工作和组织——并不是所有的脚本都有用。但阅读速度非常快。与上述结果不完全相同，但在一天结束时，我可以立即从1.3TB文件（0.5度分辨率，2.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 

（注意：尺寸的顺序已经改变，如这里所述：当使用ncdf4::ncvar_get时，我如何指定尺寸顺序？）

Answer 1

不确定有没有考虑过cdo提取点？

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

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

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

Answer 2

编辑:原始问题有一个错误，但开始阅读时也可能有不同的开销，所以多次重复是公平的。< code>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

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

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

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

在较小的样本数据集上，sst.wkmean.1990-present。nc＜/code＞，我在使用基准测试命令时看到了以下结果：

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

天真的分块简直是瞎猜；我是这样使用< code>nccopy实用程序的:

$ nccopy -c"lat/100，lon/100，time/100，nbnds/" SST . wk mean . 1990-present . NC chunked . NC

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

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

下面由Marco Hermida撰写的博客文章展示了在读取特定netCDF文件的时间序列时，不同的分块策略是如何影响速度的。这可能只能用作一个起点。

• AR-4 3D数据文件的Netcdf-4分块性能结果

关于通过nc复制重新分块显然挂起；该问题似乎与默认的块缓存大小4MB有关。通过将其增加到4GB（或更多），您可以将大文件的复制时间从超过24小时减少到不到11分钟！

• Nccopy非常慢/挂起
• Unidata JIRA Trouble Ticket System:NCF-85，改进nccopy实用程序中块缓存的使用，使其适用于重新存储大型文件

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