加快在python中读取非常大的netcdf文件

Posted 2023-02-16

【中文标题】加快在python中读取非常大的netcdf文件

【英文标题】：Speeding up reading of very large netcdf file in python

【发布时间】：2016-05-27 03:50:04

【问题描述】：

我有一个非常大的 netCDF 文件，我正在 python 中使用 netCDF4 读取它

我无法一次全部读取此文件，因为它的尺寸 (1200 x 720 x 1440) 太大，无法一次将整个文件存储在内存中。第1维代表时间，后2维分别代表纬度和经度。

import netCDF4 
nc_file = netCDF4.Dataset(path_file, 'r', format='NETCDF4')
for yr in years:
    nc_file.variables[variable_name][int(yr), :, :]

但是，一次阅读一年的速度非常慢。对于以下用例，我该如何加快速度？

--编辑

块大小为1

我可以读取一系列年份：nc_file.variables[variable_name][0:100, :, :]

有几个用例：

以年为单位：

numpy.ma.sum(nc_file.variables[variable_name][int(yr), :, :])

---

# Multiply each year by a 2D array of shape (720 x 1440)
for yr in years:
    numpy.ma.sum(nc_file.variables[variable_name][int(yr), :, :] * arr_2d)

---

# Add 2 netcdf files together 
for yr in years:
    numpy.ma.sum(nc_file.variables[variable_name][int(yr), :, :] + 
                 nc_file2.variables[variable_name][int(yr), :, :])

【问题讨论】：

您确定阅读任何其他内容（例如一次读取整个文件）会更快吗？您可以尝试使用裁剪后的文件吗？

任何essential profiling 完成了吗？

您在阅读当年的数据后是否对其进行了处理？你能读懂几年的范围吗，例如[1997:2007,:,:]?

感谢@hapulj，我能读到很多年。有几个用例。编辑问题以反映它们。

【参考方案1】：

我强烈建议您查看xarray 和dask 项目。使用这些强大的工具将允许您轻松地将计算拆分为块。这带来了两个优势：您可以计算无法放入内存的数据，并且您可以使用机器中的所有内核以获得更好的性能。 You can optimize the performance by appropriately choosing the chunk size (see documentation).

您可以通过简单的操作从 netCDF 加载数据

import xarray as xr
ds = xr.open_dataset(path_file)

如果您想沿时间维度以年为单位对数据进行分块，则指定chunks 参数（假设年份坐标名为“年”）：

ds = xr.open_dataset(path_file, chunks='year': 10)

由于其他坐标没有出现在chunks 字典中，因此将使用单个块。 （请参阅文档here 中的更多详细信息。）。这对于您的第一个要求很有用，您希望每年乘以 2D 数组。您只需这样做：

ds['new_var'] = ds['var_name'] * arr_2d

现在，xarray 和 dask 正在懒惰地计算您的结果。为了触发实际计算，您可以简单地要求xarray 将您的结果保存回netCDF：

ds.to_netcdf(new_file)

计算通过dask 触发，该dask 负责将处理分成块，从而可以处理不适合内存的数据。此外，dask 将负责使用您的所有处理器内核来计算块。

xarray 和 dask 项目仍然不能很好地处理块不能很好地“对齐”以进行并行计算的情况。由于在这种情况下，我们仅在“年份”维度中进行了分块，因此我们预计不会有任何问题。

如果你想把两个不同的netCDF文件加在一起，很简单：

ds1 = xr.open_dataset(path_file1, chunks='year': 10)
ds2 = xr.open_dataset(path_file2, chunks='year': 10)
(ds1 + ds2).to_netcdf(new_file)

我提供了一个使用a dataset available online 的完整示例。

In [1]:

import xarray as xr
import numpy as np

# Load sample data and strip out most of it:
ds = xr.open_dataset('ECMWF_ERA-40_subset.nc', chunks = 'time': 4)
ds.attrs = 
ds = ds[['latitude', 'longitude', 'time', 'tcw']]
ds

Out[1]:

<xarray.Dataset>
Dimensions:    (latitude: 73, longitude: 144, time: 62)
Coordinates:
  * latitude   (latitude) float32 90.0 87.5 85.0 82.5 80.0 77.5 75.0 72.5 ...
  * longitude  (longitude) float32 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 ...
  * time       (time) datetime64[ns] 2002-07-01T12:00:00 2002-07-01T18:00:00 ...
Data variables:
    tcw        (time, latitude, longitude) float64 10.15 10.15 10.15 10.15 ...

In [2]:

arr2d = np.ones((73, 144)) * 3.
arr2d.shape

Out[2]:

(73, 144)

In [3]:

myds = ds
myds['new_var'] = ds['tcw'] * arr2d

In [4]:

myds

Out[4]:

<xarray.Dataset>
Dimensions:    (latitude: 73, longitude: 144, time: 62)
Coordinates:
  * latitude   (latitude) float32 90.0 87.5 85.0 82.5 80.0 77.5 75.0 72.5 ...
  * longitude  (longitude) float32 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 ...
  * time       (time) datetime64[ns] 2002-07-01T12:00:00 2002-07-01T18:00:00 ...
Data variables:
    tcw        (time, latitude, longitude) float64 10.15 10.15 10.15 10.15 ...
    new_var    (time, latitude, longitude) float64 30.46 30.46 30.46 30.46 ...

In [5]:

myds.to_netcdf('myds.nc')
xr.open_dataset('myds.nc')

Out[5]:

<xarray.Dataset>
Dimensions:    (latitude: 73, longitude: 144, time: 62)
Coordinates:
  * latitude   (latitude) float32 90.0 87.5 85.0 82.5 80.0 77.5 75.0 72.5 ...
  * longitude  (longitude) float32 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 ...
  * time       (time) datetime64[ns] 2002-07-01T12:00:00 2002-07-01T18:00:00 ...
Data variables:
    tcw        (time, latitude, longitude) float64 10.15 10.15 10.15 10.15 ...
    new_var    (time, latitude, longitude) float64 30.46 30.46 30.46 30.46 ...

In [6]:

(myds + myds).to_netcdf('myds2.nc')
xr.open_dataset('myds2.nc')

Out[6]:

<xarray.Dataset>
Dimensions:    (latitude: 73, longitude: 144, time: 62)
Coordinates:
  * time       (time) datetime64[ns] 2002-07-01T12:00:00 2002-07-01T18:00:00 ...
  * latitude   (latitude) float32 90.0 87.5 85.0 82.5 80.0 77.5 75.0 72.5 ...
  * longitude  (longitude) float32 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 ...
Data variables:
    tcw        (time, latitude, longitude) float64 20.31 20.31 20.31 20.31 ...
    new_var    (time, latitude, longitude) float64 60.92 60.92 60.92 60.92 ...

【参考方案2】：

这有点像 Hacky，但可能是最简单的解决方案：

将文件的子集读入内存，然后将文件 cPickle (https://docs.python.org/3/library/pickle.html) 读回磁盘以供将来使用。从腌制数据结构加载数据可能比每次都解析 netCDF 更快。

【讨论】：

很有可能使用 blosc 压缩写入/读取 hdf5，就像在 PyTables 中一样，实际上比 cPickle 快。更不用说对于未压缩的数字数据可能会变得非常大的文件大小！

【参考方案3】：

检查文件的分块。 ncdump -s <infile> 会给出答案。如果时间维度上的块大小大于 1，您应该一次读取相同数量的年，否则您一次从磁盘读取几年并且一次只使用一个。 慢有多慢？对于这种大小的数组，每个时间步长最多几秒听起来是合理的。 稍后提供有关您如何处理数据的更多信息可能会为我们提供更多关于问题所在的指导。

【讨论】：

时间维度是1，其他维度是什么？您能否也澄清一下您的情况“慢”有多慢！

块大小为 720，其他维度为 1440。循环的每次迭代只需要几分之一秒。但是当你必须迭代超过 1200 年时，它就会加起来

那么您可能已经达到了当前文件和硬件的速度。如果您可以选择重写数据，您可以尝试 PyTables 并将文件转换为 blosc 压缩的 HDF5。这应该比 zlib 压缩的 NetCDF4 更快，尽管文件会稍大一些。由于重写文件不是您的问题中的选项，因此我不会添加它来回答，但由于我最近将 NetCDF 转换为 PyTables，我可以给您一些提示。

感谢@kakk11，重写选项有多慢/快？即，将 netcdf 重写为 hdf5 是否需要很长时间，以至于后续的速度优势无用？

尝试前的时间很难估计，大概15-30分钟吧？但是您每个文件只执行一次，以后的所有分析都可以在 hdf 文件上完成，您的所有分析都会运行得更快。您还可以在转换过程中重新分块数据，这样可以更快地读取空间子集。因此，这实际上取决于您计划读取文件的次数。您也可以尝试并行读取，但同样取决于 IO 速度，它可能无法证明额外的编码工作是合理的。