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

Speeding up reading of very large netcdf file in python

我有一个非常大的 netCDF 文件,我正在 python

中使用 netCDF4 读取它

我无法一次读取此文件,因为它的尺寸 (1200 x 720 x 1440) 太大,整个文件无法一次进入内存。第一个维度代表时间,接下来的两个维度分别代表纬度和经度。

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

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

--编辑

块大小为 1

  1. 我可以读出一个年份范围:nc_file.variables[variable_name][0:100, :, :]

  2. 有几个用例:

    年:

    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), :, :])

检查文件分块。 ncdump -s <infile>会给出答案。如果时间维度中的块大小大于 1,则您应该一次读取相同数量的年份,否则您将一次从磁盘读取若干年并且一次只使用一个年份。 慢到什么程度才算慢?对于这种大小的数组,每个时间步最多几秒听起来很合理。 提供有关您稍后如何处理数据的更多信息可能会给我们更多关于问题可能出在哪里的指导。

这有点 Hacky,但可能是最简单的解决方案:

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

我强烈建议您看一下 xarray and dask projects. Using these powerful tools will allow you to easily split up the computation in chunks. This brings up two advantages: you can compute on data which does not fit in memory, and you can use all of the cores in your machine for better performance. 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 参数(假设年份坐标命名为 'year'):

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

由于其他坐标没有出现在 chunks dict 中,因此它们将使用单个块。 (请参阅文档 here 中的更多详细信息)。这对于您的第一个要求很有用,您希望每年乘以一个二维数组。你只需要做:

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

现在,xarraydask 正在 延迟计算您的结果 。为了触发实际计算,您可以简单地要求 xarray 将结果保存回 netCDF:

ds.to_netcdf(new_file)

计算是通过 dask 触发的,它负责将处理分成块,从而可以处理不适合内存的数据。此外,dask 将负责将所有处理器内核用于计算块。

xarraydask 项目仍然不能很好地处理块不能"align" 并行计算的情况。由于在本例中我们仅在 'year' 维度中分块,因此我们预计不会出现任何问题。

如果你想把两个不同的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 ...