scripy.interpolate.griddata

如何解决scripy.interpolate.griddata

我正在处理来自具有非结构化网格的海洋模型的大文件。例如，具有维度ux 的可变流速(194988,1009,20) 分别用于：'nFlowElement'（命名非结构化网格元素）、'time' 和laydim（深度维度）。我想将这些结果插入到维度为 (600,560,20) 的结构化网格中：纬度、经度、时间和位置。为此，我使用 scipy.interpolate.griddata。由于这些数据数组太大而无法立即加载到您的工作内存中，因此我正在尝试使用“块”（dask）。不幸的是，我在尝试将 apply_ufunc 与设置一起使用时遇到了问题：dask = 'parallelized'。

对于较小的计算域（较小的 nFlowElement 维度），我仍然能够在我的工作内存中加载数据数组。然后，下面的代码给了我想要的结果：

def interp_to_grid(u,xc,yc,xint,yint):
    print(u.shape,xc.shape,xint.shape)
    ug = griddata((xc,yc),u,(xint,yint),method='nearest',fill_value=np.nan)
    return ug 

uxg = xr.apply_ufunc(interp_to_grid,ux,yint,dask = 'allowed',input_core_dims=[['nFlowElem','time','laydim'],['nFlowElem'],['dim_0','dim_1'],'dim_1']],output_core_dims=[['dim_0','dim_1','laydim']],output_dtypes = [xr.DataArray]
                     )

请注意，在函数 interp_to_grid 中，输入变量具有以下维度：

u（即ux，原始流速输出）：(194988,20) for (nFlowElem,time,laydim)
xc,yc（与这些 194988 个元素关联的纬度和经度坐标）所以两者 (194988,)
xint,yint（我想插入数据的结构化网格坐标）：对于 (dim_0,dim_1) 都是 (600,560) 请注意， scipy.interpolate.griddata 不需要我遍历时间和laydim 维度（如上面的代码所示）。为此，以正确的顺序（'time' 和 'laydim' 最后）提供 griddata 维度至关重要。内插结果 uxg 的维度为 (600、560、1009、20) - 符合预期。

但是，对于更大的空间域，需要使用 dask = 'parallelized'，因为这些输入数据数组不能再加载到我的工作内存中。我尝试在时间维度上应用块，但也在 nFlowElement 维度上应用。我知道不可能对核心维度进行分块。

这是我的“并行”尝试之一（随着时间的推移而变暗）：

输入用户体验：

<xarray.DataArray 'ucx' (nFlowElem: 194988,time: 1009,laydim: 20)>
dask.array<transpose,shape=(194988,20),dtype=float64,chunksize=(194988,10,chunktype=numpy.ndarray>
Coordinates:
    FlowElem_xcc  (nFlowElem) float64 dask.array<chunksize=(194988,),meta=np.ndarray>
    FlowElem_ycc  (nFlowElem) float64 dask.array<chunksize=(194988,meta=np.ndarray>
  * time          (time) datetime64[ns] 2014-09-17 ... 2014-10-01
Dimensions without coordinates: nFlowElem,laydim
Attributes:
    standard_name:  eastward_sea_water_velocity
    long_name:      velocity on flow element center,x-component
    units:          m s-1
    grid_mapping:   wgs84

Apply_func：

uxg = xr.apply_ufunc(interp_to_grid,dask = 'parallelized',input_core_dims=[['nFlowElem'],output_dtypes = [xr.DataArray],)

给出错误：

  File "interpnd.pyx",line 78,in scipy.interpolate.interpnd.NDInterpolatorBase.__init__

  File "interpnd.pyx",line 192,in scipy.interpolate.interpnd._check_init_shape

ValueError: different number of values and points

我在更改 apply_ufunc 中的核心维度和要分块的维度方面进行了很多尝试。此外，我还尝试手动更改数据数组 u 的维度顺序，它被“馈送到”griddata（在 interp_to_grid 中）。

非常欢迎任何建议！最良好的祝愿，卢卡

### 更新 1：###

我还尝试将 nFlowElem 维度而不是“时间”分块。然后我定义输入和输出核心尺寸，如下所示：

input_core_dims=[['time',[],

现在我“只”收到关于输出大小的错误（不再分块）。我的输出不再分块，因为插值后 nFlowElem 丢失了。因此，(600,20) 数据数组对于我的工作记忆来说太大了。

有没有办法让 apply_ufunc 沿着任何其他维度分块你的输出？

### 更新 2：###

或者，我尝试对时间维度（50 个时间步）进行分块，并删除了所有输入/输出核心维度。

如果我然后定义 interp_to_grid 如下（以获得正确的输入尺寸）：

def interp_to_grid(u,yint):

    xc = xc[:,0]
    yc = yc[:,0]
    u = u[:,:,0]
    print(u.shape,xint.shape) #input shape

    ug = griddata((xc,fill_value=np.nan)
    print(ug.shape) #output shape

    return ug

我确实为我的 interp_to_grid 输入（见第一行）和输出（见第二行）获得了正确的尺寸。但是，我得到了 ValueError: axes don't match array:

(194988,50,20) (194988,) (600,560)
(600,20)
Traceback (most recent call last):

  File "<ipython-input-11-6b65fe3dba5b>",line 1,in <module>
    uxg.loc[dict(time=np.datetime64('2014-09-18'),laydim=19)].plot()

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\xarray\plot\plot.py",line 444,in __call__
    return plot(self._da,**kwargs)

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\xarray\plot\plot.py",line 160,in plot
    darray = darray.squeeze().compute()

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\xarray\core\dataarray.py",line 899,in compute
    return new.load(**kwargs)

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\xarray\core\dataarray.py",line 873,in load
    ds = self._to_temp_dataset().load(**kwargs)

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\xarray\core\dataset.py",line 798,in load
    evaluated_data = da.compute(*lazy_data.values(),**kwargs)

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\dask\base.py",line 565,in compute
    results = schedule(dsk,keys,**kwargs)

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\dask\threaded.py",line 84,in get
    **kwargs

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\dask\local.py",line 487,in get_async
    raise_exception(exc,tb)

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\dask\local.py",line 317,in reraise
    raise exc

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\dask\local.py",line 222,in execute_task
    result = _execute_task(task,data)

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\dask\core.py",line 121,in _execute_task
    return func(*(_execute_task(a,cache) for a in args))

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\dask\core.py",in <genexpr>
    return func(*(_execute_task(a,cache) for a in args))

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\dask\optimization.py",line 963,in __call__
    return core.get(self.dsk,self.outkey,dict(zip(self.inkeys,args)))

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\dask\core.py",line 151,in get
    result = _execute_task(task,cache)

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\dask\core.py",cache) for a in args))

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\dask\utils.py",line 35,in apply
    return func(*args,**kwargs)

  File "<__array_function__ internals>",line 6,in transpose

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\numpy\core\fromnumeric.py",line 658,in transpose
    return _wrapfunc(a,'transpose',axes)

  File "C:\Users\920507\Anaconda3\envs\dfm_tools_env\lib\site-packages\numpy\core\fromnumeric.py",line 58,in _wrapfunc
    return bound(*args,**kwds)

ValueError: axes don't match array

scripy.interpolate.griddata

如何解决scripy.interpolate.griddata

相关推荐