从“isel”操作创建的Dask xarray无法加载值（速度太慢）

Question

我有一个4d dask阵列，其尺寸对应于（时间、深度、纬度、经度）。仅供参考，这是一个海洋数据集

# Create xarray dataset object for ocean temperature (time x depth x lat x lon)

DS=xr.open_mfdataset([outdir + s for s in flist],combine='by_coords',chunks={'xi_rho':25,'eta_rho':25})

DS.temp

输出：

xarray.DataArray
'temp' (ocean_time: 1456, s_rho: 50, eta_rho: 489, xi_rho: 655)
dask.array<chunksize=(1456, 50, 25, 25), meta=np.ndarray>

当我尝试从这个dask数组加载1d向量时，没有问题

T=DS.temp
%time
T.isel(ocean_time=0,eta_rho=100,xi_rho=500).values

输出（忽略以下值）：

CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 5.96 µs

我现在只想选择那些海底深度超过1000米的区域（纬度、经度）

depth_max=1e3;
deep=np.where(DS.h >=depth_max); # DS.h is the depth of the ocean bottom.

# Number of locations where the ocean is deeper than depth_max
xy_num=len(deep[0])

这给了我一个元组'deep'，它的第一个元素(deep[0])是满足条件的所有'eta_rho'（纬度索引）值的列表。元组(deep[1])的第二个元素是对应的'xi_rho'（经度索引）值的列表。因此，我可以使用(deep[0][0],deep[1][0]), (deep[0][1],deep[1][1]), (deep[0][2],deep[1][2]), (deep[0][3],deep[1][3])等构建(lat,lon)索引对

这很好，因为我想创建一个坐标，它扫过满足上述条件的(lat,lon)对。目标是从(time,depth,lat,lon) -> (time,depth,gthmax)得到gthmax是abvoe描述的新坐标。我是这样做的：

# Picking only those (lat,lon) where the condition is satisfied
T_deep=T.isel(eta_rho=xr.DataArray(deep[0],dims='gthmax'),xi_rho=xr.DataArray(deep[1],dims='gthmax'))

# Explicitly assign the new coordinate
T_deep=T_deep.assign_coords({"gthmax":range(0,xy_num)})

# Create chunks along this new coordinate
T_deep=T_deep.chunk({'gthmax':1000})

T_deep

输出（仅显示尺寸）：

xarray.DataArray 'temp' (ocean_time: 1456, s_rho: 50, gthmax: 133446)
dask.array<chunksize=(1456, 50, 1000), meta=np.ndarray>

这就是问题所在。当我尝试从这个新的3d dask数组中访问值时，即使是在一个点上，下面的命令也永远不会完成，我必须终止内核。我也尝试过使用load()和compute()，但没有效果

T_deep.isel(ocean_time=0,s_rho=46,gthmax=100).values

在从原来的4d dask阵列到3d dask阵列的转换过程中，我在哪里搞砸了

谢谢

Answer 1

我没有数据集进行测试，而且很难说有限制的信息。这是我的第一个猜测

使用xr.open_mfdataset加载信息时，实际上是在为地形变量h指定一个新维度“ocean_time”。“h”的维度应该从[eta_-rho，xi_-rho]转移到[ocean_-time，eta_-rho，xi_-rho]。因此，在工作时

depth_max=1e3  
deep = np.where( DS.h >= depth_max)

深部的尺寸不是[eta_-rho，xi_-rho]；事实上，它们也是[海洋时间，埃塔罗，希罗]。因此，我想问题就出现在这里。通过使用将坐标系从[ocean_time，s_rho，eta_rho，xi_rho]推送到[ocean_time，s_rho，gthmax]

eta_rho=xr.DataArray(deep[0],dims='gthmax')  
xi_rho=xr.DataArray(deep[1],dims='gthmax')  

请注意，维度“海洋时间”的长度为1456，远大于水平维度（eta_-rho:489，xi_-rho:655）。我相信这会混淆xarray/dask，并在需要调用值时减慢进程

---

既然我得到了你的数据集，让我更新我的答案

df = xr.open_mfdataset('./*.nc',  
                       combine='by_coords',  
                       chunks={'ocean_time':10}, 
                       parallel=True, 
)                                                                                                                 

temp = df[['temp']]                                                                                            
temp1 = temp.sel(xi_rho=xi_rho,eta_rho=eta_rho)                                                                   
%time temp1.isel(ocean_time=0,s_rho=46,z=100).compute()                                                          
CPU times: user 2.2 s, sys: 1.11 s, total: 3.31 s
Wall time: 3.05 s