[Help]: Unstacking disp to 2D Grid using suggested optimizations for sparse grids

[SteffanDavies]

I am having difficulty unstacking the cumulative displacement grids using the suggested optimization on Patreon. The velocity works fine and so do the other grids.

# unstack to 2D grid
yxunstack(disp_ps_yxstack, phase_ps_stack['stack'])

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[118], line 2
      1 # unstack to 2D grid
----> 2 yxunstack(disp_ps_yxstack, phase_ps_stack['stack'])

Cell In[24], line 7, in yxunstack(ds, stack_multiindex)
      5 def yxunstack(ds, stack_multiindex):
      6     # x index order is reverted in the stack
----> 7     return ds.squeeze(drop=True).rename(y='stack').assign_coords(stack=stack_multiindex).compute().unstack('stack').sortby('x')

File ~/miniconda3/envs/pygmtsar2/lib/python3.12/site-packages/xarray/core/dataarray.py:1169, in DataArray.compute(self, **kwargs)
   1150 """Manually trigger loading of this array's data from disk or a
   1151 remote source into memory and return a new array. The original is
   1152 left unaltered.
   (...)
   1166 dask.compute
   1167 """
   1168 new = self.copy(deep=False)
-> 1169 return new.load(**kwargs)

File ~/miniconda3/envs/pygmtsar2/lib/python3.12/site-packages/xarray/core/dataarray.py:1143, in DataArray.load(self, **kwargs)
   1125 def load(self, **kwargs) -> Self:
   1126     """Manually trigger loading of this array's data from disk or a
   1127     remote source into memory and return this array.
   1128 
   (...)
   1141     dask.compute
   1142     """
-> 1143     ds = self._to_temp_dataset().load(**kwargs)
   1144     new = self._from_temp_dataset(ds)
   1145     self._variable = new._variable

File ~/miniconda3/envs/pygmtsar2/lib/python3.12/site-packages/xarray/core/dataset.py:850, in Dataset.load(self, **kwargs)
    845     evaluated_data: tuple[np.ndarray[Any, Any], ...] = chunkmanager.compute(
    846         *lazy_data.values(), **kwargs
    847     )
    849     for k, data in zip(lazy_data, evaluated_data):
--> 850         self.variables[k].data = data
    852 # load everything else sequentially
    853 for k, v in self.variables.items():

File ~/miniconda3/envs/pygmtsar2/lib/python3.12/site-packages/xarray/core/variable.py:440, in Variable.data(self, data)
    437 @data.setter
    438 def data(self, data: T_DuckArray | ArrayLike) -> None:
    439     data = as_compatible_data(data)
--> 440     self._check_shape(data)
    441     self._data = data

File ~/miniconda3/envs/pygmtsar2/lib/python3.12/site-packages/xarray/namedarray/core.py:518, in NamedArray._check_shape(self, new_data)
    516 def _check_shape(self, new_data: duckarray[Any, _DType_co]) -> None:
    517     if new_data.shape != self.shape:
--> 518         raise ValueError(
    519             f"replacement data must match the {self.__class__.__name__}'s shape. "
    520             f"replacement data has shape {new_data.shape}; {self.__class__.__name__} has shape {self.shape}"
    521         )

ValueError: replacement data must match the Variable's shape. replacement data has shape (251, 1); Variable has shape (251, 1158725)

I have even tried to do one date at a time but the same issue remains, despite having the same shape as velocity.

[AlexeyPechnikov]

Yes, this workaround is incompatible with some functions, and I remember RMSE calculation also does not work. I plan to support 2D stack-based computations in the same way as 3D ones. For now, I use this hack myself and have shared it with all users who might find it useful.

[AlexeyPechnikov]

@SteffanDavies FYI, you can try the latest GitHub library version where most functions now support stack processing. Convert grids to stack format with stack=grid.stack(['y','x']) and convert them back with grid=stack.unstack(). Note that stack.dropna() is a simple but not lazy way to remove NaNs. I currently utilize stack processing for some of my projects, but I have no intention of showcasing the new features in public examples for the time being.