diff --git a/.gitignore b/.gitignore
index adaafc7..110bb27 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,9 +1,11 @@
+.vscode/
 .idea/
 data/
 docs/modules/
 docs/_build/
 docs/auto_examples/
 coverage/
+scratch
 
 # So far, all html files are auto-generated
 *.html
diff --git a/docs/content/target.rst b/docs/content/target.rst
index c3ae0ae..81599b0 100644
--- a/docs/content/target.rst
+++ b/docs/content/target.rst
@@ -11,7 +11,7 @@ for supervised learning. This is achieved with a :py:class:`Target` object:
 
     >>> from sklearn_xarray import wrap, Target
     >>> from sklearn_xarray.datasets import load_digits_dataarray
-    >>> from sklearn.linear_model.logistic import LogisticRegression
+    >>> from sklearn.linear_model import LogisticRegression
     >>>
     >>> X = load_digits_dataarray()
     >>> y = Target(coord='digit')(X)
@@ -61,17 +61,52 @@ Pre-processing
 --------------
 
 In some cases, it is necessary to pre-process the coordinate before it can be
-used as a target. For this, the constructor takes a ``transform_func`` parameter
-which can be used with the ``fit_transform`` method of transformers in
-``sklearn.preprocessing`` (and also any other object implementing the sklearn
-transformer interface):
+used as a target. For this, the constructor takes a ``transformer`` parameter
+which can be used with transformers in ``sklearn.preprocessing`` (and also any
+other object implementing the sklearn transformer interface):
 
 .. doctest::
 
     >>> from sklearn.neural_network import MLPClassifier
     >>> from sklearn.preprocessing import LabelBinarizer
     >>>
-    >>> y = Target(coord='digit', transform_func=LabelBinarizer().fit_transform)(X)
+    >>> y = Target(coord='digit', transformer=LabelBinarizer(), reshapes="feature")
+    >>> wrapper = wrap(MLPClassifier(), reshapes="feature")
+    >>> wrapper.fit(X, y) # doctest:+ELLIPSIS
+    EstimatorWrapper(...)
+
+This approach makes it possible to reverse the pre-processing, e.g. after
+calling ``wrapper.predict``:
+
+.. doctest::
+
+    >>> yp = wrapper.predict(X)
+    >>> yp
+    <xarray.DataArray (sample: 1797, feature: 10)>
+    array([[1, 0, 0, ..., 0, 0, 0],
+           [0, 1, 0, ..., 0, 0, 0],
+           [0, 0, 1, ..., 0, 0, 0],
+           ...,
+           [0, 0, 0, ..., 0, 1, 0],
+           [0, 0, 0, ..., 0, 0, 1],
+           [0, 0, 0, ..., 0, 1, 0]])
+    Coordinates:
+      * sample   (sample) int64 0 1 2 3 4 5 6 ... 1790 1791 1792 1793 1794 1795 1796
+        digit    (sample) int64 0 1 2 3 4 5 6 7 8 9 0 1 ... 7 9 5 4 8 8 4 9 0 8 9 8
+    Dimensions without coordinates: feature
+    >>> y.inverse_transform(yp)
+    <xarray.DataArray (sample: 1797)>
+    array([0, 1, 2, ..., 8, 9, 8])
+    Coordinates:
+      * sample   (sample) int64 0 1 2 3 4 5 6 ... 1790 1791 1792 1793 1794 1795 1796
+        digit    (sample) int64 0 1 2 3 4 5 6 7 8 9 0 1 ... 7 9 5 4 8 8 4 9 0 8 9 8
+
+
+Alternatively, the constructor also accepts a ``transform_func`` parameter:
+
+.. doctest::
+
+    >>> y = Target(coord='digit', transform_func=LabelBinarizer().fit_transform)
     >>> wrapper = wrap(MLPClassifier())
     >>> wrapper.fit(X, y) # doctest:+ELLIPSIS
     EstimatorWrapper(...)
@@ -81,13 +116,13 @@ Indexing
 
 A :py:class:`Target` object can be indexed in the same way as the underlying
 coordinate and interfaces with ``numpy`` by providing an ``__array__``
-attribute which returns ``numpy.array()`` of the (transformed) coordinate.
+attribute which returns ``numpy.array()`` of the (transformed) data.
 
 
 Multi-dimensional coordinates
 -----------------------------
 
-In some cases, the target coordinates span multiple dimensions, but the
+In some cases, the target data spans multiple dimensions, but the
 transformer expects a lower-dimensional input. With  the ``dim`` parameter of
 the :py:class:`Target` class you can specify which of the dimensions to keep.
 You can also specify the callable ``reduce_func`` to perform the reduction of
diff --git a/docs/content/transformers.rst b/docs/content/transformers.rst
index 103d333..8df8b2a 100644
--- a/docs/content/transformers.rst
+++ b/docs/content/transformers.rst
@@ -6,14 +6,83 @@ xarray's powerful array manipulation syntax. Refer to :ref:`API/Pre-processing`
 for a full list.
 
 
+Combining dimensions
+--------------------
+
+scikit-learn's estimators generally assume that data is two-dimensional:
+the first dimension represents the samples, the second dimension the features
+of your data. Since xarray is generally used for higher-dimensional data, it is
+often necessary to perform pre-processing steps that combine multiple
+dimensions to a sample and/or feature dimension, or even combine multiple
+variables of a ``Dataset`` into a single ``DataArray``.
+
+.. py:currentmodule:: sklearn_xarray.datasets
+
+For example, the :py:func:`load_digits_dataarray` method loads a
+three-dimensional array of 8-by-8-pixel grayscale images:
+
+.. doctest::
+
+    >>> from sklearn_xarray.datasets import load_digits_dataarray
+    >>> X = load_digits_dataarray(load_images=True)
+    >>> X # doctest:+ELLIPSIS
+    <xarray.DataArray (sample: 1797, row: 8, col: 8)>
+    array([[[ 0.,  0.,  5., ...,  1.,  0.,  0.],
+            [ 0.,  0., 13., ..., 15.,  5.,  0.],
+            [ 0.,  3., 15., ..., 11.,  8.,  0.],
+            ...,
+            [ 0.,  4., 16., ..., 16.,  6.,  0.],
+            [ 0.,  8., 16., ..., 16.,  8.,  0.],
+            [ 0.,  1.,  8., ..., 12.,  1.,  0.]]])
+    Coordinates:
+      * sample   (sample) int64 0 1 2 3 4 5 6 ... 1790 1791 1792 1793 1794 1795 1796
+      * row      (row) int64 0 1 2 3 4 5 6 7
+      * col      (col) int64 0 1 2 3 4 5 6 7
+        digit    (sample) int64 0 1 2 3 4 5 6 7 8 9 0 1 ... 7 9 5 4 8 8 4 9 0 8 9 8
+
+.. py:currentmodule:: sklearn_xarray.preprocessing
+
+In order to use the individual images as samples to fit an estimator, we need
+to vectorize them first. The :py:class:`Featurizer` combines all dimensions
+of the array except for the sample dimension:
+
+.. doctest::
+
+    >>> from sklearn_xarray.preprocessing import Featurizer
+    >>> Featurizer().fit_transform(X)
+    <xarray.DataArray (sample: 1797, feature: 64)>
+    array([[0., 0., 0., ..., 0., 0., 0.],
+           [0., 0., 0., ..., 0., 0., 0.],
+           [0., 0., 0., ..., 0., 0., 0.],
+           ...,
+           [0., 0., 0., ..., 0., 0., 0.],
+           [0., 0., 0., ..., 0., 0., 0.],
+           [0., 0., 0., ..., 0., 0., 0.]])
+    Coordinates:
+      * sample   (sample) int64 0 1 2 3 4 5 6 ... 1790 1791 1792 1793 1794 1795 1796
+        digit    (sample) int64 0 1 2 3 4 5 6 7 8 9 0 1 ... 7 9 5 4 8 8 4 9 0 8 9 8
+      * feature  (feature) MultiIndex
+      - col      (feature) int64 0 0 0 0 0 0 0 0 1 1 1 1 ... 6 6 6 6 7 7 7 7 7 7 7 7
+      - row      (feature) int64 0 1 2 3 4 5 6 7 0 1 2 3 ... 4 5 6 7 0 1 2 3 4 5 6 7
+
+Other transformers for combining dimensions are:
+
+.. autosummary::
+    :nosignatures:
+
+    Concatenator
+    Featurizer
+    Stacker
+
+Check out the :ref:`examples<examples>` for more use cases.
+
+
 Transformers changing the number of samples
 -------------------------------------------
 
 There are several transformers that change the number of samples in the data,
 namely:
 
-.. py:currentmodule:: sklearn_xarray.preprocessing
-
 .. autosummary::
     :nosignatures:
 
@@ -21,6 +90,7 @@ namely:
     Sanitizer
     Segmenter
     Splitter
+    Stacker
 
 These kinds of transformer are usually disallowed by sklearn, because the
 package does not provide any mechanism of also changing the number of samples
@@ -83,25 +153,25 @@ specify the ``groupby`` parameter:
     >>>
     >>> X = load_wisdm_dataarray()
     >>> Xt = segmenter.fit_transform(X)
-    >>> Xt # doctest:+ELLIPSIS doctest:+NORMALIZE_WHITESPACE
+    >>> Xt # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
     <xarray.DataArray 'tmptmp' (sample: 54813, axis: 3, timepoint: 20)>
-    array([[[ -0.15    ,   0.11    , ...,  -2.26    ,  -1.46    ],
-            [  9.15    ,   9.19    , ...,   9.72    ,   9.81    ],
-            [ -0.34    ,   2.76    , ...,   2.03    ,   2.15    ]],
-           [[  0.27    ,  -3.06    , ...,  -2.56    ,  -2.6     ],
-            [ 12.57    ,  13.18    , ...,  14.56    ,   8.96    ],
-            [  5.37    ,   6.47    , ...,   0.31    ,  -3.3     ]],
-           ...,
-           [[ -0.3     ,   0.27    , ...,   0.42    ,   3.17    ],
-            [  8.08    ,   6.63    , ...,  10.5     ,   9.23    ],
-            [  0.994285,   0.994285, ...,  -5.175732,  -4.671779]],
-           [[  5.33    ,   6.44    , ...,  -4.14    ,  -4.9     ],
-            [  8.39    ,   9.04    , ...,   6.21    ,   6.55    ],
-            [ -4.794363,  -2.179256, ...,   5.938472,   3.827318]]])
+    array([[[-0.15    ,   0.11    , ...,  -2.26    ,  -1.46    ],
+            [ 9.15    ,   9.19    , ...,   9.72    ,   9.81    ],
+            [-0.34    ,   2.76    , ...,   2.03    ,   2.15    ]],
+           [[ 0.27    ,  -3.06    , ...,  -2.56    ,  -2.6     ],
+            [12.57    ,  13.18    , ...,  14.56    ,   8.96    ],
+            [ 5.37    ,   6.47    , ...,   0.31    ,  -3.3     ]],
+           ...
+           [[-0.3     ,   0.27    , ...,   0.42    ,   3.17    ],
+            [ 8.08    ,   6.63    , ...,  10.5     ,   9.23    ],
+            [ 0.99... ,   0.99... , ...,  -5.17... ,  -4.67... ]],
+           [[ 5.33    ,   6.44    , ...,  -4.14    ,  -4.9     ],
+            [ 8.39    ,   9.04    , ...,   6.21    ,   6.55    ],
+            [-4.79... ,  -2.17... , ...,   5.93... ,   3.82... ]]])
     Coordinates:
       * axis       (axis) <U1 'x' 'y' 'z'
       * timepoint  (timepoint) int64 0 1 2 3 4 5 6 7 8 ... 12 13 14 15 16 17 18 19
-      * sample     (sample) datetime64[ns] 1970-01-01T13:25:37.050000 ... 1970-01-01T03:12:42.100000
+      * sample     (sample) datetime64[ns] 1970-01-01T13:25:37.050000 ... 1970-01...
         subject    (sample, timepoint) int64 1 1 1 1 1 1 1 ... 36 36 36 36 36 36 36
         activity   (sample, timepoint) object 'Downstairs' ... 'Walking'
 
diff --git a/docs/content/whatsnew.rst b/docs/content/whatsnew.rst
index 226a981..7a10bd8 100644
--- a/docs/content/whatsnew.rst
+++ b/docs/content/whatsnew.rst
@@ -2,6 +2,18 @@ What's New
 ==========
 
 
+v0.5.0 (unreleased)
+-------------------------
+
+Enhancements
+~~~~~~~~~~~~
+
+- New ``Stacker`` transformer that provides a transformer interface to
+  xarray's ``stack``/``unstack`` methods (thanks to @mmann1123 for the input).
+- Un-deprecated the ``transformer`` parameter of the ``Target`` class and
+  added an ``inverse_transform`` method that reverses the transformation.
+
+
 v0.4.0 (June 18, 2020)
 -------------------------
 
@@ -17,7 +29,7 @@ Enhancements
 
 - The package can now be installed via conda::
 
-    conda install -c phausamann -c conda-forge sklearn-xarray
+    conda install -c phausamann sklearn-xarray
 
 
 
diff --git a/docs/content/wrappers.rst b/docs/content/wrappers.rst
index 49ba2df..f96ac6e 100644
--- a/docs/content/wrappers.rst
+++ b/docs/content/wrappers.rst
@@ -4,7 +4,7 @@ Wrappers for sklearn estimators
 sklearn-xarray provides wrappers that let you use sklearn estimators on
 xarray DataArrays and Datasets. The goal is to provide a seamless integration
 of both packages by only applying estimator methods on the raw data while
-metadata (coordinates in xarray) remains untouched whereever possible.
+metadata (coordinates in xarray) remains untouched whenever possible.
 
 There are two principal data types in xarray: ``DataArray`` and ``Dataset``.
 The wrappers provided in this package will determine automatically which
diff --git a/environment.yml b/environment.yml
index c570be5..f4c4f64 100644
--- a/environment.yml
+++ b/environment.yml
@@ -2,7 +2,7 @@ name: sklearn-xarray
 channels:
   - conda-forge
 dependencies:
-  - python=3.6
+  - python=3.7
   - numpy
   - scipy
   - scikit-learn
@@ -10,12 +10,10 @@ dependencies:
   - xarray
   - pytest
   - matplotlib
-  - sphinx
-  - pillow
+  - sphinx==2.4.4
   - sphinx-gallery
   - sphinx_rtd_theme
   - numpydoc
-  - bump2version
-  - pre-commit
-  - flake8
-  - black=19.10b0
+  # install sklearn-xarray in development mode through pip
+  - pip:
+    - -e .
diff --git a/examples/README.txt b/examples/README.txt
index a5d244b..d975cbf 100644
--- a/examples/README.txt
+++ b/examples/README.txt
@@ -1,6 +1,4 @@
-.. _general_examples:
+.. _examples:
 
-General examples
-================
-
-Introductory examples.
+Examples
+========
diff --git a/requirements.txt b/requirements.txt
index e444199..bebd695 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,5 @@
-scikit-learn==0.23.1
-xarray==0.15.1
+scikit-learn==0.24.1
+xarray==0.16.2
 pandas==1.0.4
 numpy==1.18.5
 scipy==1.4.1
diff --git a/requirements_dev.txt b/requirements_dev.txt
index ae5b2eb..d74cd58 100644
--- a/requirements_dev.txt
+++ b/requirements_dev.txt
@@ -1,8 +1,6 @@
 pytest==5.4.3
-sphinx=2.4.4
+sphinx==2.4.4
 sphinx_rtd_theme==0.4.3
 sphinx-gallery==0.7.0
 numpydoc==1.0.0
 matplotlib==3.2.1
-pillow==7.1.2
-bump2version==1.0.0
diff --git a/sklearn_xarray/preprocessing.py b/sklearn_xarray/preprocessing.py
index 9bb4616..ed25616 100644
--- a/sklearn_xarray/preprocessing.py
+++ b/sklearn_xarray/preprocessing.py
@@ -49,13 +49,13 @@ def preprocess(X, function, groupby=None, group_dim="sample", **fit_params):
         Xt = X
 
     if groupby is None:
-        Xt = Xt.apply(function, **fit_params)
+        Xt = Xt.map(function, **fit_params)
     else:
         group_idx = get_group_indices(X, groupby, group_dim)
         Xt_list = []
         for i in group_idx:
             x = Xt.isel(**{group_dim: i})
-            Xt_list.append(x.apply(function, **fit_params))
+            Xt_list.append(x.map(function, **fit_params))
         Xt = xr.concat(Xt_list, dim=group_dim)
 
     if was_array:
@@ -412,7 +412,7 @@ def _transform(self, X):
         else:
             raise KeyError("Unrecognized mode for index reduction")
 
-        dim_coord = Xt[self.dim][dim_idx]
+        dim_coord = Xt[self.dim][dim_idx].values
 
         # keep the original coord if desired
         if self.keep_coords_as is not None:
@@ -420,7 +420,7 @@ def _transform(self, X):
 
         # get indices of new dimension
         if self.new_index_func is None:
-            new_dim_coord = Xt[self.dim][: self.new_len]
+            new_dim_coord = Xt[self.dim][: self.new_len].values
         else:
             new_dim_coord = self.new_index_func(self.new_len)
 
@@ -456,7 +456,7 @@ def _inverse_transform(self, X):
 
         if self.keep_coords_as is not None:
             Xt[tmp_dim] = Xt[self.keep_coords_as]
-            Xt = Xt.drop(self.keep_coords_as)
+            Xt = Xt.drop_vars(self.keep_coords_as)
 
         # transpose to original dimensions
         Xt = Xt.rename({tmp_dim: self.dim})
@@ -641,10 +641,12 @@ def _rebuild_array(self, arr, axis):
             else:
                 mg_ord = [0]
             idx = np.vstack(
-                self._segment_array(
-                    np.transpose(g, mg_ord), axis_old, True
-                ).flatten()
-                for g in np.meshgrid(*old_ranges)
+                [
+                    self._segment_array(
+                        np.transpose(g, mg_ord), axis_old, True
+                    ).flatten()
+                    for g in np.meshgrid(*old_ranges)
+                ]
             )
             return npg.aggregate(
                 idx, arr.flatten().T, size=old_shape, func="mean"
@@ -772,7 +774,7 @@ def _restore_coords(self, X):
             )
         }
 
-        X = X.drop(self.keep_coords_as)
+        X = X.drop_vars(self.keep_coords_as)
 
         for c in X.coords:
             if c not in (self.dim, self.new_dim) and self.dim in X[c].dims:
@@ -1298,7 +1300,7 @@ def _transform_var(self, X):
         if self.order is not None:
             stack_dims = self.order
         else:
-            stack_dims = tuple(set(X.dims) - {self.sample_dim})
+            stack_dims = tuple(sorted(set(X.dims) - {self.sample_dim}))
 
         if len(stack_dims) == 0:
             # TODO write a test for this (nothing to stack)
@@ -1593,3 +1595,96 @@ def reduce(X, return_estimator=False, **fit_params):
         return Xt, estimator
     else:
         return Xt
+
+
+class Stacker(BaseTransformer):
+    """ Stack multiple dimensions into a single one.
+
+    Parameters
+    ----------
+    stack_dims : iterable of str
+        List (tuple) of the dimensions used to define how the data is sampled.
+        If your sample dim has multiple dimensions, for instance x,y,time
+        these can be passed as a list or tuple. Before stacking, a new
+        multiindex 'sample' will be created for these dimensions.
+
+    new_dim : str
+        Name of multi-index used to stack dims.
+
+    direction : str, optional
+        "stack" or "unstack" defines the direction of transformation.
+        Default is "stack".
+
+    groupby : str or list, optional
+        Name of coordinate or list of coordinates by which the groups are
+        determined.
+
+    group_dim : str, optional
+        Name of dimension along which the groups are indexed.
+    """
+
+    def __init__(
+        self,
+        stack_dims=None,
+        new_dim=None,
+        direction="stack",
+        groupby=None,
+        group_dim="sample",
+    ):
+        if stack_dims is None:
+            raise ValueError("stack_dims cannot be None.")
+
+        if new_dim is None:
+            raise ValueError("new_dim cannot be None.")
+
+        if direction not in ["stack", "unstack"]:
+            raise ValueError("direction must be 'stack' or 'unstack'.")
+
+        self.stack_dims = stack_dims
+        self.direction = direction
+        self.new_dim = new_dim
+        self.groupby = groupby
+        self.group_dim = group_dim
+
+    def _transform(self, X):
+        """ Transform. """
+
+        if self.direction == "stack":
+            return X.stack(**{self.new_dim: self.stack_dims})
+        else:
+            return X.unstack(self.new_dim)
+
+    def _inverse_transform(self, X):
+        """ Reverse transform. """
+
+        if self.direction == "unstack":
+            return X.stack(**{self.new_dim: self.stack_dims})
+        else:
+            return X.unstack(self.new_dim)
+
+
+def stack(X, return_estimator=False, **fit_params):
+    """ Stacks multiple dimensions into a single one.
+
+    Parameters
+    ----------
+    X : xarray DataArray or Dataset
+        The input data.
+
+    return_estimator : bool
+        Whether to return the fitted estimator along with the transformed data.
+
+    Returns
+    -------
+    Xt : xarray DataArray or Dataset
+        The transformed data.
+    """
+
+    estimator = Stacker(**fit_params)
+
+    Xt = estimator.fit_transform(X)
+
+    if return_estimator:
+        return Xt, estimator
+    else:
+        return Xt
diff --git a/sklearn_xarray/target.py b/sklearn_xarray/target.py
index cd1156c..7882d7e 100644
--- a/sklearn_xarray/target.py
+++ b/sklearn_xarray/target.py
@@ -1,6 +1,8 @@
 """``sklearn_xarray.target``"""
 
 import numpy as np
+import xarray as xr
+from sklearn.utils.validation import check_is_fitted, NotFittedError
 
 
 class Target(object):
@@ -13,28 +15,29 @@ class Target(object):
         specified, the target will be the entire DataArray/Dataset.
 
     transform_func : callable, optional
-        A function that transforms the coordinate values to an
-        sklearn-compatible type and shape. If not specified, the coordinate(s)
-        will be used as-is.
+        A function that transforms the data to an sklearn-compatible type and
+        shape. If not specified, the target data will be used as-is. Mutually
+        exclusive with the `transformer` argument.
 
     transformer : sklearn transformer, optional
-        **Deprecated**, use ``transform_func=Transformer().fit_transform``
-        instead.
+        Transforms the data into an sklearn compatible target representation.
+        If not specified, the target data will be used as-is. Mutually
+        exclusive with the `transform_func` argument.
 
-    lazy : bool, optinonal
-        If true, the target coordinate is only transformed by the transformer
+    lazy : bool, default False
+        If True, the target data is only transformed by the transformer
         when needed. The transformer can implement a ``get_transformed_shape``
         method that returns the shape after the transformation of the provided
-        coordinate without actually transforming the data.
+        data without actually transforming it.
 
     dim : str or sequence of str, optional
-        When set, multi-dimensional coordinates will be reduced to this
+        If specified, multi-dimensional data will be reduced to this
         dimension/these dimensions.
 
     reduce_func : callable, optional
-        A callable that reduces the coordinate(s) to the dimension(s) in
-        ``dim``. If not specified, the values along dimensions not in ``dim``
-        will be reduced to the first element in each of these dimensions.
+        A callable that reduces the data to the dimension(s) in ``dim``. If not
+        specified, the values along dimensions not in ``dim`` will be reduced
+        to the first element in each of these dimensions.
     """
 
     def __init__(
@@ -42,30 +45,31 @@ def __init__(
         coord=None,
         transform_func=None,
         transformer=None,
+        reshapes=None,
         lazy=False,
         dim=None,
         reduce_func=None,
     ):
 
-        self.transform_func = transform_func
         self.coord = coord
+        self.transform_func = transform_func
+        self.transformer = transformer
+        self.reshapes = reshapes
         self.lazy = lazy
         self.reduce_func = reduce_func
         self.dim = dim
 
-        self.transformer = transformer
-        if transformer is not None:
-            import warnings
-
-            warnings.simplefilter("always", DeprecationWarning)
-            warnings.warn(
-                "The transformer argument is deprecated and will be removed "
-                "in a future version. Use "
-                "transform_func=Transformer().fit_transform instead.",
-                DeprecationWarning,
-            )
-            warnings.simplefilter("ignore", DeprecationWarning)
-            self.transform_func = self.transformer.fit_transform
+        if self.transformer is not None:
+            if self.transform_func is not None:
+                raise ValueError(
+                    "You can specify either the 'transformer' or the "
+                    "'transform_func' argument, not both"
+                )
+            try:
+                check_is_fitted(self.transformer)
+                self.transform_func = self.transformer.transform
+            except NotFittedError:
+                self.transform_func = self.transformer.fit_transform
 
         self.values = None
 
@@ -142,6 +146,22 @@ def _reduce(self, values):
             other_dims = [d for d in values.dims if d not in dim]
             return values.reduce(self.reduce_func, dim=other_dims)
 
+    def _update_coords(self, X):
+        """ Update the coordinates of a reshaped DataArray. """
+
+        coords_new = dict()
+
+        # drop all coords along the reshaped dimensions
+        for c in X.coords:
+            if self.reshapes in X[c].dims and c != self.reshapes:
+                c_t = X[c].isel(**{self.reshapes: 0})
+                new_dims = [d for d in X[c].dims if d != self.reshapes]
+                coords_new[c] = (new_dims, c_t.drop(self.reshapes))
+            elif c != self.reshapes:
+                coords_new[c] = X[c]
+
+        return coords_new
+
     @property
     def shape(self):
         """ The shape of the transformed target. """
@@ -195,3 +215,34 @@ def assign_to(self, X):
             self.values = self.transform_func(self.values)
 
         return self
+
+    def inverse_transform(self, X):
+        """ Reverse the transformation performed by the transformer.
+
+        Parameters
+        ----------
+        X : xarray DataArray
+            The input data.
+
+        Returns
+        -------
+        Xt : xarray DataArray
+            The transformed data.
+        """
+        if self.transformer is None:
+            raise ValueError(
+                "Cannot use inverse_transform when the target was constructed "
+                "without the 'transformer' parameter"
+            )
+
+        data = self.transformer.inverse_transform(X)
+        if self.reshapes is not None:
+            coords = self._update_coords(X)
+        else:
+            coords = X.coords
+
+        dims = list(X.dims)
+        if data.ndim < X.ndim and self.reshapes is not None:
+            dims.remove(self.reshapes)
+
+        return xr.DataArray(data, coords, dims, name=X.name, attrs=X.attrs)
diff --git a/tests/test_preprocessing.py b/tests/test_preprocessing.py
index 3d8b405..b0a6331 100644
--- a/tests/test_preprocessing.py
+++ b/tests/test_preprocessing.py
@@ -2,6 +2,7 @@
 import xarray as xr
 import xarray.testing as xrt
 import numpy.testing as npt
+import pytest
 
 from sklearn_xarray.preprocessing import (
     preprocess,
@@ -14,6 +15,7 @@
     select,
     sanitize,
     reduce,
+    stack,
     Splitter,
 )
 
@@ -42,7 +44,7 @@ def test_preprocess():
 
     Xt_ds = preprocess(X_ds, scale)
 
-    xrt.assert_allclose(Xt_ds, X_ds.apply(scale))
+    xrt.assert_allclose(Xt_ds, X_ds.map(scale))
 
 
 def test_groupwise():
@@ -386,3 +388,53 @@ def test_reduce():
     Xt_da = reduce(X_da)
 
     xrt.assert_allclose(Xt_da, X_da.reduce(np.linalg.norm, dim="feature"))
+
+
+def test_stack():
+
+    X_da = xr.DataArray(
+        np.random.random((100, 10, 10)),
+        coords={
+            "sample": range(100),
+            "feat_1": range(10),
+            "feat_2": range(10),
+        },
+        dims=("sample", "feat_1", "feat_2"),
+    )
+
+    Xt_da, stacker = stack(
+        X_da,
+        stack_dims=["feat_1", "feat_2"],
+        new_dim="feature",
+        return_estimator=True,
+    )
+
+    assert Xt_da.shape == (100, 100)
+    assert Xt_da.dims == ("sample", "feature")
+
+    xrt.assert_equal(stacker.inverse_transform(Xt_da), X_da)
+
+    # unstack
+    X_da_act, unstacker = stack(
+        Xt_da,
+        stack_dims=["feat_1", "feat_2"],
+        new_dim="feature",
+        direction="unstack",
+        return_estimator=True,
+    )
+
+    xrt.assert_equal(unstacker.inverse_transform(X_da_act), Xt_da)
+
+    with pytest.raises(ValueError):
+        stack(X_da, new_dim="feature")
+
+    with pytest.raises(ValueError):
+        stack(X_da, stack_dims=["feat_1", "feat_2"])
+
+    with pytest.raises(ValueError):
+        stack(
+            X_da,
+            new_dim="feature",
+            stack_dims=["feat_1", "feat_2"],
+            direction="not_a_direction",
+        )
diff --git a/tests/test_target.py b/tests/test_target.py
index 451a1da..d022663 100644
--- a/tests/test_target.py
+++ b/tests/test_target.py
@@ -1,4 +1,5 @@
 import numpy as np
+import pytest
 import xarray as xr
 import numpy.testing as npt
 
@@ -25,13 +26,14 @@ def test_constructor():
 
     target = Target(transform_func=convert_to_ndarray)
     target.assign_to(X_ds)
-
     npt.assert_equal(target.values, np.array(X_ds.var_1))
 
     target = Target(coord="coord_1", transformer=LabelBinarizer())(X_ds)
-
     npt.assert_equal(target.values, LabelBinarizer().fit_transform(coord_1))
 
+    with pytest.raises(ValueError):
+        Target(transform_func=convert_to_ndarray, transformer=LabelBinarizer())
+
 
 def test_str():
 
@@ -63,11 +65,9 @@ def test_array():
         },
     )
 
-    target = Target(
-        coord="coord_1",
-        transform_func=LabelBinarizer().fit_transform,
-        lazy=True,
-    )(X_ds)
+    target = Target(coord="coord_1", transformer=LabelBinarizer(), lazy=True)(
+        X_ds
+    )
 
     npt.assert_equal(np.array(target), LabelBinarizer().fit_transform(coord_1))
 
@@ -159,9 +159,36 @@ def test_multidim_coord():
         transform_func=LabelBinarizer().fit_transform,
         dim="sample",
     )(X_ds)
+    npt.assert_equal(target_1, LabelBinarizer().fit_transform(coord_1[:, 0]))
+
     target_2 = Target(
         coord="coord_2", dim=["sample", "feat_1"], reduce_func=np.mean
     )(X_ds)
-
-    npt.assert_equal(target_1, LabelBinarizer().fit_transform(coord_1[:, 0]))
     npt.assert_equal(target_2, np.mean(coord_2, 2))
+
+
+def test_inverse_transform():
+
+    X_da = xr.DataArray(
+        np.random.random((100, 10)),
+        coords={
+            "sample": range(100),
+            "feature": range(10),
+            "coord_1": (["sample"], ["a"] * 34 + ["b"] * 33 + ["c"] * 33),
+        },
+        dims=("sample", "feature"),
+    )
+
+    target = Target(
+        coord="coord_1", transformer=LabelBinarizer(), reshapes="feature"
+    )(X_da)
+
+    y_da = xr.DataArray(
+        target.values,
+        coords={"sample": range(100)},
+        dims=("sample", "feature"),
+    )
+
+    yt_da = target.inverse_transform(y_da)
+    yt_exp = X_da.coord_1.drop("coord_1")
+    xr.testing.assert_equal(yt_da, yt_exp)