Repository: pydata/xarray

Issue:
Convolution operation
I would like to perform 2d convolution operation with dask-data stored in xarray, defining the kernel.

Related issue are:
https://github.com/pydata/xarray/issues/1142
https://github.com/pydata/xarray/issues/2010
https://github.com/pydata/xarray/issues/1323

Referring to Dask, they have a specific image-oriented API which allows for convolutions:
http://image.dask.org/en/latest/dask_image.ndfilters.html

And this guy already implemented convolution operation based on xarray and dask:
https://github.com/serazing
https://github.com/serazing/xscale
https://xscale.readthedocs.io/en/latest/api.html

Are there plans to implement this feature?
In the meanwhile I'll try dask-image and xscale.

Discussion:
max-sixty: This looks very promising; I'm surprised it was possible without more code. I'm being slow, but where is the nd-rolling algo? I had thought bottleneck didn't support more than one dimension? https://bottleneck.readthedocs.io/en/latest/bottleneck.move.html, and that we'd have to implement our own in numbagg (which would be very possible)
max-sixty: Re the API, I think the dict is probably the best option, although it does complicate as the arguments become differently typed depending on one vs multiple dimensions.

One alternative is to allow fluent args, like:
```
(
    da
    .rolling(x=3, center=True)
    .rolling(y=5, min_periods=2)
    .mean()
)
```
...but does that then seem like the second `rolling` is operating on the _result_ of the first?
fujiisoup: Hi @max-sixty 

> One alternative is to allow fluent args, like:
> ...but does that then seem like the second rolling is operating on the result of the first?

I couldn't think of it until just now. But yes, it sounds to me like a repeated rolling operation.

> I'm being slow, but where is the nd-rolling algo? I had thought bottleneck didn't support more than one dimension? 

No. With nd-rolling, we need to use numpy reductions.
Its `skipna=True` operation is currently slow, but it can be improved replacing nan *before* the stride-trick.


fujiisoup: Another API concern. We now use `min_periods`, in which we implicitly assume one-dimension cases.

With nd-dimension, I think `min_counts` argument is more appropriate like bottleneck, which will limit the lower bound of the number of missing entries in the n-dimensional window.

Even if we leave it, we may disallow nd-argument of `min_periods` but keep it a scalar.
pep8speaks: Hello @fujiisoup! Thanks for updating this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:





* In the file [`xarray/core/rolling.py`](https://github.com/pydata/xarray/blob/f74b2e8ea9588e2932376fd09f9bfc02e270d34d/xarray/core/rolling.py):

> [Line 93:1](https://github.com/pydata/xarray/blob/f74b2e8ea9588e2932376fd09f9bfc02e270d34d/xarray/core/rolling.py#L93): [W293](https://duckduckgo.com/?q=pep8%20W293) blank line contains whitespace












dcherian: I think the dictionary is OK. We can allow scalars for `center`, `stride` etc. in which case the scalar applies to all dimensions. This would preserve backward compatibility.

scalar for `min_counts` also seems OK to me. We can update it if someone requests it.
fujiisoup: I think now it is ready for review, though I'm sure tests miss a lot of edge cases.
Maybe we can fix them if pointed out.
fujiisoup: A possible improvement will be nan-reduction methods for nd-rolling.
Currently, we just use numpy nan-reductions, which is memory consuming for strided arrays.

This issue can be solved by replacing nan by appropriate values and applying nonnan-reduction methods, 
e.g.,
```python
da.rolling(x=2, y=3).construct(x='xw', y='yw').sum(['xw', 'yw'])
```
should be the same with
```python
da.rolling(x=2, y=3).construct(x='xw', y='yw', fill_value=0).sum(['xw', 'yw'], skipna=False)
```
and the latter is much more memory efficient.

I'd like to leave this improvement to future PR.
fujiisoup: I got an error for typechecking, only in CI but not in local, from the code that I didn't change.
fujiisoup: Could anyone kindly review this?
fujiisoup: Thanks @max-sixty for the review ;)
I'll work for the update in a few days.
max-sixty: I'm still a bit confused. How does the "roll over each dimension in turn" approach equal the "roll over both dimension together" approach with a function like std? Here's a proposed counter example:

```python
import xarray as xr
import numpy as np

da = xr.DataArray(np.asarray([[0,10,0],[0,10,0], [0,10,0]]), dims=list('xy'))
print(da)

<xarray.DataArray (x: 3, y: 3)>
array([[ 0, 10,  0],
       [ 0, 10,  0],
       [ 0, 10,  0]])
Dimensions without coordinates: x, y

x_std = da.rolling(dict(x=2)).std()
print(x_std)

<xarray.DataArray (x: 3, y: 3)>
array([[nan, nan, nan],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])
Dimensions without coordinates: x, y

x_then_y_std = x_std.rolling(dict(y=2)).std()
print(x_then_y_std)

<xarray.DataArray (x: 3, y: 3)>
array([[nan, nan, nan],
       [nan,  0.,  0.],
       [nan,  0.,  0.]])
Dimensions without coordinates: x, y

combined_std = da.rolling(dict(x=2, y=2)).std()
print(combined_std)

<xarray.DataArray (x: 3, y: 3)>
array([[nan, nan, nan],
       [nan,  5.,  5.],
       [nan,  5.,  5.]])
Dimensions without coordinates: x, y

```
dcherian: I think it's more like
```
da.rolling(x=2).construct("xr").rolling(y=2).construct("yr").std(["xr", "yr"])
```

but that gives
```
array([[0., 5., 5.],
       [0., 5., 5.],
       [0., 5., 5.]])
```
fujiisoup: Thanks @max-sixty .

You are completely correct.
As the test pass, I was fooling myself. 

The reason was that the dataset I was using for the test does not have `time` and `x` simultaneously.
So I was not testing the 2d-rolling but just 1d-rolling.

Fixed. Now it correctly fails for `mean` and `std`, but passes with `max` and `sum`

max-sixty: Great! Thanks @fujiisoup ! 

Ready to go from my POV 

Patch:
diff --git a/doc/computation.rst b/doc/computation.rst
index 3660aed9..474c3905 100644
--- a/doc/computation.rst
+++ b/doc/computation.rst
@@ -188,9 +188,16 @@ a value when aggregating:
     r = arr.rolling(y=3, center=True, min_periods=2)
     r.mean()
 
+From version 0.17, xarray supports multidimensional rolling,
+
+.. ipython:: python
+
+    r = arr.rolling(x=2, y=3, min_periods=2)
+    r.mean()
+
 .. tip::
 
-   Note that rolling window aggregations are faster and use less memory when bottleneck_ is installed. This only applies to numpy-backed xarray objects.
+   Note that rolling window aggregations are faster and use less memory when bottleneck_ is installed. This only applies to numpy-backed xarray objects with 1d-rolling.
 
 .. _bottleneck: https://github.com/pydata/bottleneck/
 
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index c4cf931b..9d4261da 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -25,6 +25,9 @@ Breaking changes
 
 New Features
 ~~~~~~~~~~~~
+- :py:meth:`~xarray.DataArray.rolling` and :py:meth:`~xarray.Dataset.rolling`
+  now accept more than 1 dimension.(:pull:`4219`)
+  By `Keisuke Fujii <https://github.com/fujiisoup>`_.
 - Build :py:meth:`CFTimeIndex.__repr__` explicitly as :py:class:`pandas.Index`. Add ``calendar`` as a new
   property for :py:class:`CFTimeIndex` and show ``calendar`` and ``length`` in
   :py:meth:`CFTimeIndex.__repr__` (:issue:`2416`, :pull:`4092`)
diff --git a/xarray/core/common.py b/xarray/core/common.py
index c95df773..bc5035b6 100644
--- a/xarray/core/common.py
+++ b/xarray/core/common.py
@@ -786,7 +786,7 @@ class DataWithCoords(SupportsArithmetic, AttrAccessMixin):
         self,
         dim: Mapping[Hashable, int] = None,
         min_periods: int = None,
-        center: bool = False,
+        center: Union[bool, Mapping[Hashable, bool]] = False,
         keep_attrs: bool = None,
         **window_kwargs: int,
     ):
@@ -802,7 +802,7 @@ class DataWithCoords(SupportsArithmetic, AttrAccessMixin):
             Minimum number of observations in window required to have a value
             (otherwise result is NA). The default, None, is equivalent to
             setting min_periods equal to the size of the window.
-        center : boolean, default False
+        center : boolean, or a mapping, default False
             Set the labels at the center of the window.
         keep_attrs : bool, optional
             If True, the object's attributes (`attrs`) will be copied from
diff --git a/xarray/core/dask_array_ops.py b/xarray/core/dask_array_ops.py
index 87f64635..74474f43 100644
--- a/xarray/core/dask_array_ops.py
+++ b/xarray/core/dask_array_ops.py
@@ -32,69 +32,80 @@ def rolling_window(a, axis, window, center, fill_value):
     """
     import dask.array as da
 
+    if not hasattr(axis, "__len__"):
+        axis = [axis]
+        window = [window]
+        center = [center]
+
     orig_shape = a.shape
-    if axis < 0:
-        axis = a.ndim + axis
     depth = {d: 0 for d in range(a.ndim)}
-    depth[axis] = int(window / 2)
-    # For evenly sized window, we need to crop the first point of each block.
-    offset = 1 if window % 2 == 0 else 0
-
-    if depth[axis] > min(a.chunks[axis]):
-        raise ValueError(
-            "For window size %d, every chunk should be larger than %d, "
-            "but the smallest chunk size is %d. Rechunk your array\n"
-            "with a larger chunk size or a chunk size that\n"
-            "more evenly divides the shape of your array."
-            % (window, depth[axis], min(a.chunks[axis]))
-        )
-
-    # Although da.overlap pads values to boundaries of the array,
-    # the size of the generated array is smaller than what we want
-    # if center == False.
-    if center:
-        start = int(window / 2)  # 10 -> 5,  9 -> 4
-        end = window - 1 - start
-    else:
-        start, end = window - 1, 0
-    pad_size = max(start, end) + offset - depth[axis]
-    drop_size = 0
-    # pad_size becomes more than 0 when the overlapped array is smaller than
-    # needed. In this case, we need to enlarge the original array by padding
-    # before overlapping.
-    if pad_size > 0:
-        if pad_size < depth[axis]:
-            # overlapping requires each chunk larger than depth. If pad_size is
-            # smaller than the depth, we enlarge this and truncate it later.
-            drop_size = depth[axis] - pad_size
-            pad_size = depth[axis]
-        shape = list(a.shape)
-        shape[axis] = pad_size
-        chunks = list(a.chunks)
-        chunks[axis] = (pad_size,)
-        fill_array = da.full(shape, fill_value, dtype=a.dtype, chunks=chunks)
-        a = da.concatenate([fill_array, a], axis=axis)
-
+    offset = [0] * a.ndim
+    drop_size = [0] * a.ndim
+    pad_size = [0] * a.ndim
+    for ax, win, cent in zip(axis, window, center):
+        if ax < 0:
+            ax = a.ndim + ax
+        depth[ax] = int(win / 2)
+        # For evenly sized window, we need to crop the first point of each block.
+        offset[ax] = 1 if win % 2 == 0 else 0
+
+        if depth[ax] > min(a.chunks[ax]):
+            raise ValueError(
+                "For window size %d, every chunk should be larger than %d, "
+                "but the smallest chunk size is %d. Rechunk your array\n"
+                "with a larger chunk size or a chunk size that\n"
+                "more evenly divides the shape of your array."
+                % (win, depth[ax], min(a.chunks[ax]))
+            )
+
+        # Although da.overlap pads values to boundaries of the array,
+        # the size of the generated array is smaller than what we want
+        # if center == False.
+        if cent:
+            start = int(win / 2)  # 10 -> 5,  9 -> 4
+            end = win - 1 - start
+        else:
+            start, end = win - 1, 0
+        pad_size[ax] = max(start, end) + offset[ax] - depth[ax]
+        drop_size[ax] = 0
+        # pad_size becomes more than 0 when the overlapped array is smaller than
+        # needed. In this case, we need to enlarge the original array by padding
+        # before overlapping.
+        if pad_size[ax] > 0:
+            if pad_size[ax] < depth[ax]:
+                # overlapping requires each chunk larger than depth. If pad_size is
+                # smaller than the depth, we enlarge this and truncate it later.
+                drop_size[ax] = depth[ax] - pad_size[ax]
+                pad_size[ax] = depth[ax]
+
+    # TODO maybe following two lines can be summarized.
+    a = da.pad(
+        a, [(p, 0) for p in pad_size], mode="constant", constant_values=fill_value
+    )
     boundary = {d: fill_value for d in range(a.ndim)}
 
     # create overlap arrays
     ag = da.overlap.overlap(a, depth=depth, boundary=boundary)
 
-    # apply rolling func
-    def func(x, window, axis=-1):
+    def func(x, window, axis):
         x = np.asarray(x)
-        rolling = nputils._rolling_window(x, window, axis)
-        return rolling[(slice(None),) * axis + (slice(offset, None),)]
-
-    chunks = list(a.chunks)
-    chunks.append(window)
+        index = [slice(None)] * x.ndim
+        for ax, win in zip(axis, window):
+            x = nputils._rolling_window(x, win, ax)
+            index[ax] = slice(offset[ax], None)
+        return x[tuple(index)]
+
+    chunks = list(a.chunks) + window
+    new_axis = [a.ndim + i for i in range(len(axis))]
     out = ag.map_blocks(
-        func, dtype=a.dtype, new_axis=a.ndim, chunks=chunks, window=window, axis=axis
+        func, dtype=a.dtype, new_axis=new_axis, chunks=chunks, window=window, axis=axis
     )
 
     # crop boundary.
-    index = (slice(None),) * axis + (slice(drop_size, drop_size + orig_shape[axis]),)
-    return out[index]
+    index = [slice(None)] * a.ndim
+    for ax in axis:
+        index[ax] = slice(drop_size[ax], drop_size[ax] + orig_shape[ax])
+    return out[tuple(index)]
 
 
 def least_squares(lhs, rhs, rcond=None, skipna=False):
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index 9fbaf747..68f76c7a 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -5968,7 +5968,7 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
                     skipna_da = np.any(da.isnull())
 
             dims_to_stack = [dimname for dimname in da.dims if dimname != dim]
-            stacked_coords = {}
+            stacked_coords: Dict[Hashable, DataArray] = {}
             if dims_to_stack:
                 stacked_dim = utils.get_temp_dimname(dims_to_stack, "stacked")
                 rhs = da.transpose(dim, *dims_to_stack).stack(
diff --git a/xarray/core/nputils.py b/xarray/core/nputils.py
index fa6df63e..4f592eb3 100644
--- a/xarray/core/nputils.py
+++ b/xarray/core/nputils.py
@@ -135,14 +135,22 @@ class NumpyVIndexAdapter:
 def rolling_window(a, axis, window, center, fill_value):
     """ rolling window with padding. """
     pads = [(0, 0) for s in a.shape]
-    if center:
-        start = int(window / 2)  # 10 -> 5,  9 -> 4
-        end = window - 1 - start
-        pads[axis] = (start, end)
-    else:
-        pads[axis] = (window - 1, 0)
+    if not hasattr(axis, "__len__"):
+        axis = [axis]
+        window = [window]
+        center = [center]
+
+    for ax, win, cent in zip(axis, window, center):
+        if cent:
+            start = int(win / 2)  # 10 -> 5,  9 -> 4
+            end = win - 1 - start
+            pads[ax] = (start, end)
+        else:
+            pads[ax] = (win - 1, 0)
     a = np.pad(a, pads, mode="constant", constant_values=fill_value)
-    return _rolling_window(a, window, axis)
+    for ax, win in zip(axis, window):
+        a = _rolling_window(a, win, ax)
+    return a
 
 
 def _rolling_window(a, window, axis=-1):
diff --git a/xarray/core/rolling.py b/xarray/core/rolling.py
index ecba5307..5f996565 100644
--- a/xarray/core/rolling.py
+++ b/xarray/core/rolling.py
@@ -75,40 +75,32 @@ class Rolling:
         -------
         rolling : type of input argument
         """
-        if len(windows) != 1:
-            raise ValueError("exactly one dim/window should be provided")
-
-        dim, window = next(iter(windows.items()))
-
-        if window <= 0:
-            raise ValueError("window must be > 0")
-
+        self.dim, self.window = [], []
+        for d, w in windows.items():
+            self.dim.append(d)
+            if w <= 0:
+                raise ValueError("window must be > 0")
+            self.window.append(w)
+
+        self.center = self._mapping_to_list(center, default=False)
         self.obj = obj
 
         # attributes
-        self.window = window
         if min_periods is not None and min_periods <= 0:
             raise ValueError("min_periods must be greater than zero or None")
-        self.min_periods = min_periods
 
-        self.center = center
-        self.dim = dim
+        self.min_periods = np.prod(self.window) if min_periods is None else min_periods
 
         if keep_attrs is None:
             keep_attrs = _get_keep_attrs(default=False)
         self.keep_attrs = keep_attrs
 
-    @property
-    def _min_periods(self):
-        return self.min_periods if self.min_periods is not None else self.window
-
     def __repr__(self):
         """provide a nice str repr of our rolling object"""
 
         attrs = [
             "{k}->{v}".format(k=k, v=getattr(self, k))
-            for k in self._attributes
-            if getattr(self, k, None) is not None
+            for k in list(self.dim) + self.window + self.center + [self.min_periods]
         ]
         return "{klass} [{attrs}]".format(
             klass=self.__class__.__name__, attrs=",".join(attrs)
@@ -143,11 +135,29 @@ class Rolling:
 
     def count(self):
         rolling_count = self._counts()
-        enough_periods = rolling_count >= self._min_periods
+        enough_periods = rolling_count >= self.min_periods
         return rolling_count.where(enough_periods)
 
     count.__doc__ = _ROLLING_REDUCE_DOCSTRING_TEMPLATE.format(name="count")
 
+    def _mapping_to_list(
+        self, arg, default=None, allow_default=True, allow_allsame=True
+    ):
+        if utils.is_dict_like(arg):
+            if allow_default:
+                return [arg.get(d, default) for d in self.dim]
+            else:
+                for d in self.dim:
+                    if d not in arg:
+                        raise KeyError("argument has no key {}.".format(d))
+                return [arg[d] for d in self.dim]
+        elif allow_allsame:  # for single argument
+            return [arg] * len(self.dim)
+        elif len(self.dim) == 1:
+            return [arg]
+        else:
+            raise ValueError("Mapping argument is necessary.")
+
 
 class DataArrayRolling(Rolling):
     __slots__ = ("window_labels",)
@@ -196,33 +206,41 @@ class DataArrayRolling(Rolling):
             obj, windows, min_periods=min_periods, center=center, keep_attrs=keep_attrs
         )
 
-        self.window_labels = self.obj[self.dim]
+        # TODO legacy attribute
+        self.window_labels = self.obj[self.dim[0]]
 
     def __iter__(self):
+        if len(self.dim) > 1:
+            raise ValueError("__iter__ is only supported for 1d-rolling")
         stops = np.arange(1, len(self.window_labels) + 1)
-        starts = stops - int(self.window)
-        starts[: int(self.window)] = 0
+        starts = stops - int(self.window[0])
+        starts[: int(self.window[0])] = 0
         for (label, start, stop) in zip(self.window_labels, starts, stops):
-            window = self.obj.isel(**{self.dim: slice(start, stop)})
+            window = self.obj.isel(**{self.dim[0]: slice(start, stop)})
 
-            counts = window.count(dim=self.dim)
-            window = window.where(counts >= self._min_periods)
+            counts = window.count(dim=self.dim[0])
+            window = window.where(counts >= self.min_periods)
 
             yield (label, window)
 
-    def construct(self, window_dim, stride=1, fill_value=dtypes.NA):
+    def construct(
+        self, window_dim=None, stride=1, fill_value=dtypes.NA, **window_dim_kwargs
+    ):
         """
         Convert this rolling object to xr.DataArray,
         where the window dimension is stacked as a new dimension
 
         Parameters
         ----------
-        window_dim: str
-            New name of the window dimension.
-        stride: integer, optional
+        window_dim: str or a mapping, optional
+            A mapping from dimension name to the new window dimension names.
+            Just a string can be used for 1d-rolling.
+        stride: integer or a mapping, optional
             Size of stride for the rolling window.
         fill_value: optional. Default dtypes.NA
             Filling value to match the dimension size.
+        **window_dim_kwargs : {dim: new_name, ...}, optional
+            The keyword arguments form of ``window_dim``.
 
         Returns
         -------
@@ -251,13 +269,27 @@ class DataArrayRolling(Rolling):
 
         from .dataarray import DataArray
 
+        if window_dim is None:
+            if len(window_dim_kwargs) == 0:
+                raise ValueError(
+                    "Either window_dim or window_dim_kwargs need to be specified."
+                )
+            window_dim = {d: window_dim_kwargs[d] for d in self.dim}
+
+        window_dim = self._mapping_to_list(
+            window_dim, allow_default=False, allow_allsame=False
+        )
+        stride = self._mapping_to_list(stride, default=1)
+
         window = self.obj.variable.rolling_window(
             self.dim, self.window, window_dim, self.center, fill_value=fill_value
         )
         result = DataArray(
-            window, dims=self.obj.dims + (window_dim,), coords=self.obj.coords
+            window, dims=self.obj.dims + tuple(window_dim), coords=self.obj.coords
+        )
+        return result.isel(
+            **{d: slice(None, None, s) for d, s in zip(self.dim, stride)}
         )
-        return result.isel(**{self.dim: slice(None, None, stride)})
 
     def reduce(self, func, **kwargs):
         """Reduce the items in this group by applying `func` along some
@@ -300,27 +332,36 @@ class DataArrayRolling(Rolling):
                [ 4.,  9., 15., 18.]])
 
         """
-        rolling_dim = utils.get_temp_dimname(self.obj.dims, "_rolling_dim")
+        rolling_dim = {
+            d: utils.get_temp_dimname(self.obj.dims, "_rolling_dim_{}".format(d))
+            for d in self.dim
+        }
         windows = self.construct(rolling_dim)
-        result = windows.reduce(func, dim=rolling_dim, **kwargs)
+        result = windows.reduce(func, dim=list(rolling_dim.values()), **kwargs)
 
         # Find valid windows based on count.
         counts = self._counts()
-        return result.where(counts >= self._min_periods)
+        return result.where(counts >= self.min_periods)
 
     def _counts(self):
         """ Number of non-nan entries in each rolling window. """
 
-        rolling_dim = utils.get_temp_dimname(self.obj.dims, "_rolling_dim")
+        rolling_dim = {
+            d: utils.get_temp_dimname(self.obj.dims, "_rolling_dim_{}".format(d))
+            for d in self.dim
+        }
         # We use False as the fill_value instead of np.nan, since boolean
         # array is faster to be reduced than object array.
         # The use of skipna==False is also faster since it does not need to
         # copy the strided array.
         counts = (
             self.obj.notnull()
-            .rolling(center=self.center, **{self.dim: self.window})
+            .rolling(
+                center={d: self.center[i] for i, d in enumerate(self.dim)},
+                **{d: w for d, w in zip(self.dim, self.window)},
+            )
             .construct(rolling_dim, fill_value=False)
-            .sum(dim=rolling_dim, skipna=False)
+            .sum(dim=list(rolling_dim.values()), skipna=False)
         )
         return counts
 
@@ -329,39 +370,40 @@ class DataArrayRolling(Rolling):
 
         # bottleneck doesn't allow min_count to be 0, although it should
         # work the same as if min_count = 1
+        # Note bottleneck only works with 1d-rolling.
         if self.min_periods is not None and self.min_periods == 0:
             min_count = 1
         else:
             min_count = self.min_periods
 
-        axis = self.obj.get_axis_num(self.dim)
+        axis = self.obj.get_axis_num(self.dim[0])
 
         padded = self.obj.variable
-        if self.center:
+        if self.center[0]:
             if isinstance(padded.data, dask_array_type):
                 # Workaround to make the padded chunk size is larger than
                 # self.window-1
-                shift = -(self.window + 1) // 2
-                offset = (self.window - 1) // 2
+                shift = -(self.window[0] + 1) // 2
+                offset = (self.window[0] - 1) // 2
                 valid = (slice(None),) * axis + (
                     slice(offset, offset + self.obj.shape[axis]),
                 )
             else:
-                shift = (-self.window // 2) + 1
+                shift = (-self.window[0] // 2) + 1
                 valid = (slice(None),) * axis + (slice(-shift, None),)
-            padded = padded.pad({self.dim: (0, -shift)}, mode="constant")
+            padded = padded.pad({self.dim[0]: (0, -shift)}, mode="constant")
 
         if isinstance(padded.data, dask_array_type):
             raise AssertionError("should not be reachable")
             values = dask_rolling_wrapper(
-                func, padded.data, window=self.window, min_count=min_count, axis=axis
+                func, padded.data, window=self.window[0], min_count=min_count, axis=axis
             )
         else:
             values = func(
-                padded.data, window=self.window, min_count=min_count, axis=axis
+                padded.data, window=self.window[0], min_count=min_count, axis=axis
             )
 
-        if self.center:
+        if self.center[0]:
             values = values[valid]
         result = DataArray(values, self.obj.coords)
 
@@ -378,8 +420,10 @@ class DataArrayRolling(Rolling):
             )
             del kwargs["dim"]
 
-        if bottleneck_move_func is not None and not isinstance(
-            self.obj.data, dask_array_type
+        if (
+            bottleneck_move_func is not None
+            and not isinstance(self.obj.data, dask_array_type)
+            and len(self.dim) == 1
         ):
             # TODO: renable bottleneck with dask after the issues
             # underlying https://github.com/pydata/xarray/issues/2940 are
@@ -412,7 +456,7 @@ class DatasetRolling(Rolling):
             Minimum number of observations in window required to have a value
             (otherwise result is NA). The default, None, is equivalent to
             setting min_periods equal to the size of the window.
-        center : boolean, default False
+        center : boolean, or a mapping from dimension name to boolean, default False
             Set the labels at the center of the window.
         keep_attrs : bool, optional
             If True, the object's attributes (`attrs`) will be copied from
@@ -431,15 +475,22 @@ class DatasetRolling(Rolling):
         DataArray.groupby
         """
         super().__init__(obj, windows, min_periods, center, keep_attrs)
-        if self.dim not in self.obj.dims:
+        if any(d not in self.obj.dims for d in self.dim):
             raise KeyError(self.dim)
         # Keep each Rolling object as a dictionary
         self.rollings = {}
         for key, da in self.obj.data_vars.items():
             # keeps rollings only for the dataset depending on slf.dim
-            if self.dim in da.dims:
+            dims, center = [], {}
+            for i, d in enumerate(self.dim):
+                if d in da.dims:
+                    dims.append(d)
+                    center[d] = self.center[i]
+
+            if len(dims) > 0:
+                w = {d: windows[d] for d in dims}
                 self.rollings[key] = DataArrayRolling(
-                    da, windows, min_periods, center, keep_attrs
+                    da, w, min_periods, center, keep_attrs
                 )
 
     def _dataset_implementation(self, func, **kwargs):
@@ -447,7 +498,7 @@ class DatasetRolling(Rolling):
 
         reduced = {}
         for key, da in self.obj.data_vars.items():
-            if self.dim in da.dims:
+            if any(d in da.dims for d in self.dim):
                 reduced[key] = func(self.rollings[key], **kwargs)
             else:
                 reduced[key] = self.obj[key]
@@ -491,19 +542,29 @@ class DatasetRolling(Rolling):
             **kwargs,
         )
 
-    def construct(self, window_dim, stride=1, fill_value=dtypes.NA, keep_attrs=None):
+    def construct(
+        self,
+        window_dim=None,
+        stride=1,
+        fill_value=dtypes.NA,
+        keep_attrs=None,
+        **window_dim_kwargs,
+    ):
         """
         Convert this rolling object to xr.Dataset,
         where the window dimension is stacked as a new dimension
 
         Parameters
         ----------
-        window_dim: str
-            New name of the window dimension.
+        window_dim: str or a mapping, optional
+            A mapping from dimension name to the new window dimension names.
+            Just a string can be used for 1d-rolling.
         stride: integer, optional
             size of stride for the rolling window.
         fill_value: optional. Default dtypes.NA
             Filling value to match the dimension size.
+        **window_dim_kwargs : {dim: new_name, ...}, optional
+            The keyword arguments form of ``window_dim``.
 
         Returns
         -------
@@ -512,19 +573,35 @@ class DatasetRolling(Rolling):
 
         from .dataset import Dataset
 
+        if window_dim is None:
+            if len(window_dim_kwargs) == 0:
+                raise ValueError(
+                    "Either window_dim or window_dim_kwargs need to be specified."
+                )
+            window_dim = {d: window_dim_kwargs[d] for d in self.dim}
+
+        window_dim = self._mapping_to_list(
+            window_dim, allow_default=False, allow_allsame=False
+        )
+        stride = self._mapping_to_list(stride, default=1)
+
         if keep_attrs is None:
             keep_attrs = _get_keep_attrs(default=True)
 
         dataset = {}
         for key, da in self.obj.data_vars.items():
-            if self.dim in da.dims:
+            # keeps rollings only for the dataset depending on slf.dim
+            dims = [d for d in self.dim if d in da.dims]
+            if len(dims) > 0:
+                wi = {d: window_dim[i] for i, d in enumerate(self.dim) if d in da.dims}
+                st = {d: stride[i] for i, d in enumerate(self.dim) if d in da.dims}
                 dataset[key] = self.rollings[key].construct(
-                    window_dim, fill_value=fill_value
+                    window_dim=wi, fill_value=fill_value, stride=st
                 )
             else:
                 dataset[key] = da
         return Dataset(dataset, coords=self.obj.coords).isel(
-            **{self.dim: slice(None, None, stride)}
+            **{d: slice(None, None, s) for d, s in zip(self.dim, stride)}
         )
 
 
diff --git a/xarray/core/variable.py b/xarray/core/variable.py
index d13de439..1f86a403 100644
--- a/xarray/core/variable.py
+++ b/xarray/core/variable.py
@@ -1883,11 +1883,14 @@ class Variable(
         Parameters
         ----------
         dim: str
-            Dimension over which to compute rolling_window
+            Dimension over which to compute rolling_window.
+            For nd-rolling, should be list of dimensions.
         window: int
             Window size of the rolling
+            For nd-rolling, should be list of integers.
         window_dim: str
             New name of the window dimension.
+            For nd-rolling, should be list of integers.
         center: boolean. default False.
             If True, pad fill_value for both ends. Otherwise, pad in the head
             of the axis.
@@ -1921,15 +1924,21 @@ class Variable(
             dtype = self.dtype
             array = self.data
 
-        new_dims = self.dims + (window_dim,)
+        if isinstance(dim, list):
+            assert len(dim) == len(window)
+            assert len(dim) == len(window_dim)
+            assert len(dim) == len(center)
+        else:
+            dim = [dim]
+            window = [window]
+            window_dim = [window_dim]
+            center = [center]
+        axis = [self.get_axis_num(d) for d in dim]
+        new_dims = self.dims + tuple(window_dim)
         return Variable(
             new_dims,
             duck_array_ops.rolling_window(
-                array,
-                axis=self.get_axis_num(dim),
-                window=window,
-                center=center,
-                fill_value=fill_value,
+                array, axis=axis, window=window, center=center, fill_value=fill_value
             ),
         )
 


Test patch:
diff --git a/xarray/tests/test_dataarray.py b/xarray/tests/test_dataarray.py
index e0da3f15..7ccf1eb1 100644
--- a/xarray/tests/test_dataarray.py
+++ b/xarray/tests/test_dataarray.py
@@ -6193,8 +6193,6 @@ def test_rolling_properties(da):
     assert rolling_obj.obj.get_axis_num("time") == 1
 
     # catching invalid args
-    with pytest.raises(ValueError, match="exactly one dim/window should"):
-        da.rolling(time=7, x=2)
     with pytest.raises(ValueError, match="window must be > 0"):
         da.rolling(time=-2)
     with pytest.raises(ValueError, match="min_periods must be greater than zero"):
@@ -6399,6 +6397,47 @@ def test_rolling_count_correct():
         assert_equal(result, expected)
 
 
+@pytest.mark.parametrize("da", (1,), indirect=True)
+@pytest.mark.parametrize("center", (True, False))
+@pytest.mark.parametrize("min_periods", (None, 1))
+@pytest.mark.parametrize("name", ("sum", "mean", "max"))
+def test_ndrolling_reduce(da, center, min_periods, name):
+    rolling_obj = da.rolling(time=3, x=2, center=center, min_periods=min_periods)
+
+    actual = getattr(rolling_obj, name)()
+    expected = getattr(
+        getattr(
+            da.rolling(time=3, center=center, min_periods=min_periods), name
+        )().rolling(x=2, center=center, min_periods=min_periods),
+        name,
+    )()
+
+    assert_allclose(actual, expected)
+    assert actual.dims == expected.dims
+
+
+@pytest.mark.parametrize("center", (True, False, (True, False)))
+@pytest.mark.parametrize("fill_value", (np.nan, 0.0))
+def test_ndrolling_construct(center, fill_value):
+    da = DataArray(
+        np.arange(5 * 6 * 7).reshape(5, 6, 7).astype(float),
+        dims=["x", "y", "z"],
+        coords={"x": ["a", "b", "c", "d", "e"], "y": np.arange(6)},
+    )
+    actual = da.rolling(x=3, z=2, center=center).construct(
+        x="x1", z="z1", fill_value=fill_value
+    )
+    if not isinstance(center, tuple):
+        center = (center, center)
+    expected = (
+        da.rolling(x=3, center=center[0])
+        .construct(x="x1", fill_value=fill_value)
+        .rolling(z=2, center=center[1])
+        .construct(z="z1", fill_value=fill_value)
+    )
+    assert_allclose(actual, expected)
+
+
 def test_raise_no_warning_for_nan_in_binary_ops():
     with pytest.warns(None) as record:
         xr.DataArray([1, 2, np.NaN]) > 0
diff --git a/xarray/tests/test_dataset.py b/xarray/tests/test_dataset.py
index 9037013c..da7621dc 100644
--- a/xarray/tests/test_dataset.py
+++ b/xarray/tests/test_dataset.py
@@ -5881,8 +5881,6 @@ def test_rolling_keep_attrs():
 
 def test_rolling_properties(ds):
     # catching invalid args
-    with pytest.raises(ValueError, match="exactly one dim/window should"):
-        ds.rolling(time=7, x=2)
     with pytest.raises(ValueError, match="window must be > 0"):
         ds.rolling(time=-2)
     with pytest.raises(ValueError, match="min_periods must be greater than zero"):
@@ -6007,6 +6005,66 @@ def test_rolling_reduce(ds, center, min_periods, window, name):
         assert src_var.dims == actual[key].dims
 
 
+@pytest.mark.parametrize("ds", (2,), indirect=True)
+@pytest.mark.parametrize("center", (True, False))
+@pytest.mark.parametrize("min_periods", (None, 1))
+@pytest.mark.parametrize("name", ("sum", "max"))
+@pytest.mark.parametrize("dask", (True, False))
+def test_ndrolling_reduce(ds, center, min_periods, name, dask):
+    if dask and has_dask:
+        ds = ds.chunk({"x": 4})
+
+    rolling_obj = ds.rolling(time=4, x=3, center=center, min_periods=min_periods)
+
+    actual = getattr(rolling_obj, name)()
+    expected = getattr(
+        getattr(
+            ds.rolling(time=4, center=center, min_periods=min_periods), name
+        )().rolling(x=3, center=center, min_periods=min_periods),
+        name,
+    )()
+    assert_allclose(actual, expected)
+    assert actual.dims == expected.dims
+
+    # Do it in the opposite order
+    expected = getattr(
+        getattr(
+            ds.rolling(x=3, center=center, min_periods=min_periods), name
+        )().rolling(time=4, center=center, min_periods=min_periods),
+        name,
+    )()
+
+    assert_allclose(actual, expected)
+    assert actual.dims == expected.dims
+
+
+@pytest.mark.parametrize("center", (True, False, (True, False)))
+@pytest.mark.parametrize("fill_value", (np.nan, 0.0))
+@pytest.mark.parametrize("dask", (True, False))
+def test_ndrolling_construct(center, fill_value, dask):
+    da = DataArray(
+        np.arange(5 * 6 * 7).reshape(5, 6, 7).astype(float),
+        dims=["x", "y", "z"],
+        coords={"x": ["a", "b", "c", "d", "e"], "y": np.arange(6)},
+    )
+    ds = xr.Dataset({"da": da})
+    if dask and has_dask:
+        ds = ds.chunk({"x": 4})
+
+    actual = ds.rolling(x=3, z=2, center=center).construct(
+        x="x1", z="z1", fill_value=fill_value
+    )
+    if not isinstance(center, tuple):
+        center = (center, center)
+    expected = (
+        ds.rolling(x=3, center=center[0])
+        .construct(x="x1", fill_value=fill_value)
+        .rolling(z=2, center=center[1])
+        .construct(z="z1", fill_value=fill_value)
+    )
+    assert_allclose(actual, expected)
+
+
 def test_raise_no_warning_for_nan_in_binary_ops():
     with pytest.warns(None) as record:
         Dataset(data_vars={"x": ("y", [1, 2, np.NaN])}) > 0
diff --git a/xarray/tests/test_nputils.py b/xarray/tests/test_nputils.py
index 1002a9dd..ccb825dc 100644
--- a/xarray/tests/test_nputils.py
+++ b/xarray/tests/test_nputils.py
@@ -1,4 +1,5 @@
 import numpy as np
+import pytest
 from numpy.testing import assert_array_equal
 
 from xarray.core.nputils import NumpyVIndexAdapter, _is_contiguous, rolling_window
@@ -47,3 +48,19 @@ def test_rolling():
     actual = rolling_window(x, axis=-1, window=3, center=False, fill_value=0.0)
     expected = np.stack([expected, expected * 1.1], axis=0)
     assert_array_equal(actual, expected)
+
+
+@pytest.mark.parametrize("center", [[True, True], [False, False]])
+@pytest.mark.parametrize("axis", [(0, 1), (1, 2), (2, 0)])
+def test_nd_rolling(center, axis):
+    x = np.arange(7 * 6 * 8).reshape(7, 6, 8).astype(float)
+    window = [3, 3]
+    actual = rolling_window(
+        x, axis=axis, window=window, center=center, fill_value=np.nan
+    )
+    expected = x
+    for ax, win, cent in zip(axis, window, center):
+        expected = rolling_window(
+            expected, axis=ax, window=win, center=cent, fill_value=np.nan
+        )
+    assert_array_equal(actual, expected)
diff --git a/xarray/tests/test_testing.py b/xarray/tests/test_testing.py
index 39ad2502..adc29a3c 100644
--- a/xarray/tests/test_testing.py
+++ b/xarray/tests/test_testing.py
@@ -37,7 +37,7 @@ def test_allclose_regression():
     "obj1,obj2",
     (
         pytest.param(
-            xr.Variable("x", [1e-17, 2]), xr.Variable("x", [0, 3]), id="Variable",
+            xr.Variable("x", [1e-17, 2]), xr.Variable("x", [0, 3]), id="Variable"
         ),
         pytest.param(
             xr.DataArray([1e-17, 2], dims="x"),