15

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

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

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

我想为每个网格点提取一个时间序列，例如从特定纬度和经度开始的 10 年或 30 年。但是，我发现这非常慢。例如，尽管我可以在 0.002 秒内从单个时间点读取 10000 个值的全局切片（维度的顺序是纬度、经度、时间）：

## 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）变量的顺序是纬度、经度、时间，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

如何优化某一点的阅读时间序列，而不是一个时间点的大面积层？例如，如果文件以不同的方式编写，例如时间、纬度、经度，会更快吗？有没有一种“简单”的方法来转换 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.5 秒内 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

（注意：维度的顺序已经改变，如下所述：How can I specify dimension order when using ncdf4::ncvar_get?）

4

3 回答 3

11

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

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

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

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.wkmean.1990-present.nc chunked.nc

该nccopy实用程序的 Unidata 文档可在此处找到。

我希望我可以推荐一种特定的策略来分块您的数据集，但它高度依赖于数据。希望上面链接的文章能让您深入了解如何分块数据以实现您正在寻找的结果！

更新

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

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

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

于 2013-11-12T23:30:44.177 回答
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
于 2013-11-12T21:03:02.180 回答
2

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

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

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

cdo mergetime point_${yyyy}.nc point_series.nc
于 2017-06-07T20:20:24.917 回答