根据维度名称提取netcdf4变量切片

Question

我有一个4维的netCDF文件。我想通过给出维度之一

的名称从netCDF文件中提取一个切片

我知道如何按位置来做。例如

from netCDF4 import Dataset
hndl_nc = Dataset(path_to_nc)
# Access by slice
hndl_nc.variables['name_variable'][:,5,:,:]

已知维度的名称，如A, B, C, D。如何按维度名称而不是位置访问?

Answer 1

可以使用xarray的索引功能按维度名称访问netcdf数据。

import xarray as xr
ds = xr.open_dataset('./foo.nc')
var = ds['name_variable']
# Slice var by Dimension "A" between values 0 and 5
var_slice = var.sel(A=slice(0,5))

Answer 2

目前最接近的解决方案似乎是

np.take(nc4_variable[:],dim_ids,axis=dim)

或

nc4_variable[:].take(dim_ids,axis=dim)

，其中dim_ids是切片的列表或元组，dim是想要切片的维度。不幸的是，这似乎首先加载整个数据集，似乎没有办法;[:]是必要的。在第一种方法中忽略它，在不调整add_offset、_FillValue等参数的情况下加载数据;在第二个方法中忽略它会产生错误。

在Ipython中使用%timeit进行测试，确认了正常切片和np.take方法之间的主要差异。

希望有人能想出一个更完整的答案;将是非常有用的不同的数据集。

Answer 3

所以，我可能已经想出了一些可以作为"解决方案"的东西。

numpy数组显然可以用可迭代对象的单例列表进行索引，例如

a = np.reshape(range(0,16),(4,4),order='F')
a = a[ [[0,1], [1]] ]

返回a等于array([4,5])。另一个例子是[[range(3),[1 2],3]]。这些单例列表以*subscripts的方式展开，就好像您直接查询了a[[0,1],1]而不是a[ [[0,1],1] ]。

因此，如果您能够在netCDF变量中查询每个维度的位置和长度(对于nc_fid[var].dimension和nc_fid[var].shape非常容易)，那么您可以根据每个维度的位置简单地排列列表。例如，如果你有形状时间按经度按纬度的数据，并且你想要所有经度，所有纬度和时间索引t=5，你可以使用像

这样的东西

order_want = ['lon', 'lat', 'time'] # must figure out dimension names a priori
nlon = nc_fid[var].shape[nc_fid[var].dimensions.index('lon')]
nlat = nc_fid[var].shape[nc_fid[var].dimensions.index('lat')]
ids = [ range(0,nlon), range(0,nlat), 5 ]
ids_permute = [order_want.index(n) for n in nc_fid[var].dimensions] 
ids_query = [l[i] for l,i in zip(ids,ids_permute)]
sliced_data = nc_fid[var][list_query]

不需要预先知道维度的位置，也不需要加载变量的所有维度。

请注意，在IPython中对%timeit进行了一些测试后，似乎对全整数索引有一些特殊的延迟，例如list_query = [0,0,0,0]将花费80ms，而list_query = [range(1),0,0,0]甚至list_query = [[0,1,2,3,4,5],0,0,0]将花费1ms。非常神秘的;无论如何，显然您应该尝试确保list_query不仅仅是一个整数列表。