Repository: pySTEPS/pysteps (language: python)

Issue:
Implement a data model
Pysteps operates on 2d grids that contain metadata (such as the grid definition, units, applied transformations and temporal resolution). Currently this information is carried in a separate metadata dictionary or as separate arguments (see e.g. pysteps.nowcasts.steps.forecast) in a very inconsistent way. It would be a good idea to store all this information into the same object.

Possible solution: use xarray (http://xarray.pydata.org/en/stable) that allows storing metadata into attributes, as well as using named columns and integration with dask.
`dB_transform` ignores existing metadata
Original report
----------------
(by Valentina Gregori)

> I also noticed another small issue (already in the previous version of pysteps). When I do a conversion to reflectivity using the function conversion.to_reflectivity(R, metadata) the output value of metadata['zerovalue'] is not as expected. I think this happens because in the transformation.dB_transform function, which is called by conversion.to_reflectivity, the code that checks if "zerovalue" is a key of metadata is missing:
> ```
> if zerovalue is None:
>     zerovalue = metadata.get("zerovalue", threshold-5)
> ```
> for now, it seems that the "zerovalue" is only replaced by "threshold-5". Do you think it makes sense?
> Thanks in advance!




Discussion / hints:
This is definitely needed and an excellent idea! 

The use of xarray can indeed represent an ideal solution to this issue, at the cost of introducing  an additional external dependency. xarray also offers a natural interface to netCDF files, pysteps output format.

However, this is not an easy task and will require considerable work. @pulkkins how do you plan to introduce such data model?
May I suggest another alternative? Instead of using an xarray object, a class can be used to handle a collection of radar composite images. Internally, this class can store the radar 2D data,  metadata, the array with the observations quality and  the grid information. The benefit over an xarray object is that more method can be added to the class (for example the import, export, or maybe data preprocessing methods, etc.). In case that [xarray](http://xarray.pydata.org/en/stable/) need to be used, the class can include a method to export part of the data to a [DataArray](http://xarray.pydata.org/en/stable/generated/xarray.DataArray.html#xarray.DataArray) object.

This is how the class could look like (with an execution example at the end):
```python
import numpy

from pysteps.io import (import_fmi_pgm, import_bom_rf3,
                        import_mch_gif, import_odim_hdf5)


class RadarComposites(object):
    """
    Class to store a collection of 2D precipitation fields.
    """

    def __init__(self):
        self.times = None  # Not used for now, but this could be useful to retrieve 2D fields at a given times.
        self.data = None  # precip data
        self.quality = None
        self.latitude = None  # Not used for now
        self.longitude = None  # Not used for now
        self.metadata = None

    def add(self, file_paths, file_format='fmi', **kwargs):
        """ Add radar 2D composites to the collection."""

        file_format = file_format.lower()

        available_importers = dict(bom_rf3=import_bom_rf3,
                                   bom=import_bom_rf3,
                                   fmi_pgm=import_fmi_pgm,
                                   fmi=import_fmi_pgm,
                                   mch_gif=import_mch_gif,
                                   odim_hdf5=import_odim_hdf5)

        if file_format not in available_importers:
            raise ValueError("Format not supported.\n"
                             + "The supported formats are: "
                             + str(list(available_importers.keys())))

        importer = available_importers[file_format]
        data = list()
        quality = list()
        metadata = list()

        for file in file_paths:
            _data, _quality, _metadata = importer(file, **kwargs)

            metadata.append(_metadata)

            data.append(_data)
            quality.append(_quality)

        data = numpy.asarray(data)
        quality = numpy.asarray(quality)
        if self.data is None:
            self.data = data
        else:
            self.data = numpy.concatenate((self.data, data), axis=0)

        if self.quality is None:
            self.quality = quality
        else:
            self.quality = numpy.concatenate((self.quality, quality), axis=0)

        # For now check that the metadata is the same for each time
        equal_metadata = numpy.all([metadata[0] == value
                                    for value in metadata])

        if not equal_metadata:
            print("Warning: Metadata differs for between some files!")

        self.metadata = metadata[0]


if __name__ == '__main__':
    #Usage example
    from os.path import join
    import glob

    data_dir = "../data/fmi/20170509"

    files = glob.glob(join(data_dir, '*.pgm.gz'))

    my_composites = RadarComposites()
    #Add composites to the collection
    my_composites.add(files, file_format='fmi', gzipped=True)
```

Another idea in the same line, is to create a sub-class from the [xarray.Dataset](http://xarray.pydata.org/en/stable/generated/xarray.Dataset.html) class with the functionalities that I mention before.





So far we have tried to avoid using classes. The reason was simply that we wanted to keep the code as simple and linear as possible. But as you nicely suggest, classes can offer a way to manage data and metadata, eventually improving the workflow. It is worth considering the option @aperezhortal suggests.

Hi everybody! I'm from the Hong Kong Observatory, currently responsible for the Community SWIRLS nowcast system, which is in operational use in meteorological services in China, India, Malaysia and South Africa. We're also working with those in the Philippines, Vietnam, Myanmar and soon other Southeast Asia countries.

From our (somewhat painful) experience, I fully agree that it is necessary to adopt a single and consistent data model, and that xarray is a good choice. You might be interested to know that we've recently developed some algorithms to read in radar data files and returns xarray.DataArray objects. Currently supported formats include ESRI ASCII, HDF5/ODIM, IRIS (raw and reflectivity), NetCDF (three variants used by Philippines and Vietnam), and UF. The project website is at https://com-swirls.org/. Kindly let me know if you would like to know more.
Hi @wcwoo, many thanks for your comment. We can certainly learn a lot from your experience with the SWIRLS system! 
I think xarray still represents to us the best way to go and you are confirming us this impression.

I wasn't aware of the community version com-SWIRLS, I will certainly have a look at the project website that you've provided. 


Sorry that I missed your reply, as I forgot that I setup an auto email sorting rule some time ago. :sweat_smile: 

I wrote [this page](https://docs.com-swirls.org/data-model-and-xarray.html) to describe Com-SWIRLS' data model. The gist is to make development divisible and modules **swappable**. I hope one day we can run various QPF algorithms simultaneously to benchmark performance or to generate super-ensemble.
Is there some solution for this ahead? I would fully support the xarray approach. There is also no need to subclass Dataset but xarray Accessors can be used instead.

This would greatly simplify interaction with other packages.

Just today I had to go the long way importing/exporting data/results to/from pysteps to xarray to conveniently examine it using eg hvplot and feed it to my further processing chain.

If you consider xarray useful I would try to setup a Pull Request or example code. 
I've created a [simple gist](https://gist.github.com/kmuehlbauer/645e42a53b30752230c08c20a9c964f9), which shows some capabilities of xarray and how it could be used for pysteps.

If there is real interest here, I can try to setup a simple workflow for some LK method etc. which uses xarray.pipe technique.

Hi @kmuehlbauer !
Many thanks for your input and the very interesting example! 

We have discussed adopting xarray in the past and I think the main point and concern was whether we wanted pysteps to remain a general library of nowcasting methods or instead to introduce a "strict" data structure such as xarray.

I think that there are pros and cons for both ways and that probably we are already a bit in between. Something that we currently really miss is a proper way to associate data and metadata in the workflow. 

xarray offers extremely attractive solutions to this and other problems (e.g. interface to netCDF). Moreover, it appears to me that it is becoming somewhat of a standard in the weather radar community. This might further justify its adoption into pysteps. 

This is certainly a topic we need to consider and discuss further and I'm really glad that someone with your experience can help us out to take the right decisions!
@kmuehlbauer  yes. I also support the use of xarray. 
xarray is so good loading and saving netCDF and ploting variables as you are showing in your gist. 
However to wrap the core rutines of pySteps into xarray may require more work. Or maybe no, if xarray.pipe can help. I would be interested on see an example of the xarray.pipe tehcnique.

Also, as soon as I have a bit of spare time I will try to add some working examples on how connect some of pysteps verification rutines with xarray dataarrays. It really simplify a lot the code but on the other hand, it may obfuscate a bit the essence of the method.
Thanks for your immediate feedback @dnerini, @cvelascof. 

I might be a bit biased since I just did the implementation for xarray based polar data handling and processing within wradlib (CfRadial1, 2 and ODIM_H5). (side note: You can also read hdf5 using xarray, I'll add a example for reading HDF5 Odim later next week). And I found it extremely useful. The ability to use xarray Datasets from polar radar data over gridding to feeding into pysteps and other packages is crucial for my (hence the bias) workflow. There is also much momentum around xarray with respect to projection information and how to handle this across multiple projects. 

Neat integration with dask is also crucial if you think of reprocessing years of data.

I would not directly speak of a "strict" structure. From my experience it is easy and versatile. And you get the netcdf output for free.

@cvelascof 
For xarray.pipe the called function need to take a Dataset/Dataarray as first argument and has to return a Dataset/Dataarray. From first glance this is compatible with most of pysteps calling and can be realised using simple checks. Still this will need some thinking.

Update: And, sometimes I forget this, you all have done a great job providing this useful package.


Added simple [LK workflow for GIF images](https://gist.github.com/kmuehlbauer/645e42a53b30752230c08c20a9c964f9#file-xarray_pysteps_gif-ipynb) from meteo swiss. Shows xarray DatasetAccessor feature and some xarray.pipe. Also xarray.rolling capability is shown. 
Trying to advance this issue, I've drafted a PR which aims at including some (partial) support for xarray within pysteps, as well as a tutorial (based on @kmuehlbauer's excellent gist) to show how xarray can be integrated into pysteps.

Any help and feedback is very welcomed!

PR: https://github.com/pySTEPS/pysteps/pull/123 
I've been pondering how to make our modules more inter-operable.  Imagine that we manage to standardize the data interface, then there would be no need for any wrapper.  The user script could be slimmer and cleaner.

For example, wradlib's codes would take in a radar file and output an xarray (DataArray or Dataset), which then can be passed directly to pySTEPS' or SwirlsPy's modules for nowcasting.  The outputs of pySTEP and SwirlsPy would also be in xarray and so could then be directly used for visualization, hydrological model or nowcast warning applications etc.

To achieve this, not only do we need xarray, but we also need to think if and how to standardize the data.  There are questions like:
1. What is stored in the xarray?  e.g. reflectivity, rainfall rate, rainfall depth (over how long?)
1. What is the unit of the above?  Shall we adopt the CF convention?
1. What does the timestamp represent? i.e. the start time or the end time of the observation window, or depending on the data?  Currently, most satellite and radar are marked with the start time, but rain gauge data are marked with the end time, at least in my place.
1.  What is the dimension of the xarray?  For nowcasting it probably starts at 2D+time, but before long we want to do 3D+time, either for better use of radar data or for blending with NWP.  Then for probabilistic nowcast we would need to add one more dimension for "scenarios".  So, is 5D the maximum?
> To achieve this, not only do we need xarray, but we also need to think if and how to standardize the data. There are questions like:
> 
>     1. What is stored in the xarray?  e.g. reflectivity, rainfall rate, rainfall depth (over how long?)

This is something the user should decide. It's just a question of how much magic you want to put in the implementation. Keep it simple might be a good approach.

>     2. What is the unit of the above?  Shall we adopt the CF convention?

I would be in favour of adopting CF conventions. This is standard for most atmospheric science data.

>     3. What does the timestamp represent? i.e. the start time or the end time of the observation window, or depending on the data?  Currently, most satellite and radar are marked with the start time, but rain gauge data are marked with the end time, at least in my place.

This is implementation detail, for `start_time` aligned measurements the data is valid from `start_time` to next `start_time`, for `end_time` aligned measurements it's valid from prior `end_time` to current `end_time`. It depends on measurement frequency, so there need to be some alignment.
> 
>     4. What is the dimension of the xarray?  For nowcasting it probably starts at 2D+time, but before long we want to do 3D+time, either for better use of radar data or for blending with NWP.  Then for probabilistic nowcast we would need to add one more dimension for "scenarios".  So, is 5D the maximum?

As long as we invent clear dimension naming scheme this should also be no big issue.

Just my 2C

> This is something the user should decide. It's just a question of how much magic you want to put in the implementation. Keep it simple might be a good approach.

Agreed that we should keep it simple, at least for the start. 

> I would be in favour of adopting CF conventions. This is standard for most atmospheric science data.

Agreed.

> This is implementation detail, for `start_time` aligned measurements the data is valid from `start_time` to next `start_time`, for `end_time` aligned measurements it's valid from prior `end_time` to current `end_time`. It depends on measurement frequency, so there need to be some alignment.

Indeed it would be clearer to use start_time and end_time, instead of just time.  Some modules would have to deal with the change from start_time to end_time, e.g. when deriving rainfall amount from radar.  Given that it is an already established convention to mark satellite & radar with start_time and rain gauge & lightning with end_time, it is perhaps better to keep things they are, and be careful while performing those steps.

> As long as we invent clear dimension naming scheme this should also be no big issue.

Sure.

Recently we've received a support request from Cemaden (Brazil) to use Rainbow data with SwirlsPy.  We could develop a swirlspy.rad.rainbow module, but may also take this opportunity to try implementing a common data model.  @kmuehlbauer: I wonder if it is possible for you to develop a method in wradlib that takes in Rainbow data files and outputs an xarray.Dataset object?  If so, we can develop a method in swirlspy to performing nowcasting with that object.

I'm tagging @SelokYeung and @cwjchin as they may help implement the changes in swirlspy.



Fix patch:
diff --git a/doc/source/pysteps_reference/utils.rst b/doc/source/pysteps_reference/utils.rst
index 0622e6bbe..ad4c8213e 100644
--- a/doc/source/pysteps_reference/utils.rst
+++ b/doc/source/pysteps_reference/utils.rst
@@ -9,7 +9,6 @@ Implementation of miscellaneous utility functions.
 .. automodule:: pysteps.utils.arrays
 .. automodule:: pysteps.utils.cleansing
 .. automodule:: pysteps.utils.conversion
-.. automodule:: pysteps.utils.dimension
 .. automodule:: pysteps.utils.fft
 .. automodule:: pysteps.utils.images
 .. automodule:: pysteps.utils.interpolate
diff --git a/examples/LK_buffer_mask.py b/examples/LK_buffer_mask.py
index 821f54f19..8f0c628b8 100644
--- a/examples/LK_buffer_mask.py
+++ b/examples/LK_buffer_mask.py
@@ -55,7 +55,7 @@
 
 # Read the radar composites
 importer = io.get_method(importer_name, "importer")
-R, quality, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
+precip, quality, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
 
 del quality  # Not used
 
@@ -64,18 +64,18 @@
 # ~~~~~~~~~~~~~~~~~~~
 
 # Convert to mm/h
-R, metadata = conversion.to_rainrate(R, metadata)
+precip = precip.pysteps.to_rainrate()
 
 # Keep the reference frame in mm/h and its mask (for plotting purposes)
-ref_mm = R[0, :, :].copy()
+ref_mm = precip[0, :, :].copy()
 mask = np.ones(ref_mm.shape)
 mask[~np.isnan(ref_mm)] = np.nan
 
 # Log-transform the data [dBR]
-R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)
+precip = precip.pysteps.db_transform()
 
 # Keep the reference frame in dBR (for plotting purposes)
-ref_dbr = R[0].copy()
+ref_dbr = precip[0].copy()
 ref_dbr[ref_dbr < -10] = np.nan
 
 # Plot the reference field
@@ -106,11 +106,11 @@
 dense_lucaskanade = motion.get_method("LK")
 
 # Mask invalid values
-R = np.ma.masked_invalid(R)
+precip = np.ma.masked_invalid(precip)
 
 # Use no buffering of the radar mask
 fd_kwargs1 = {"buffer_mask": 0}
-xy, uv = dense_lucaskanade(R, dense=False, fd_kwargs=fd_kwargs1)
+xy, uv = dense_lucaskanade(precip, dense=False, fd_kwargs=fd_kwargs1)
 plt.imshow(ref_dbr, cmap=plt.get_cmap("Greys"))
 plt.imshow(mask, cmap=colors.ListedColormap(["black"]), alpha=0.5)
 plt.quiver(
@@ -142,7 +142,7 @@
 # with buffer
 buffer = 10
 fd_kwargs2 = {"buffer_mask": buffer}
-xy, uv = dense_lucaskanade(R, dense=False, fd_kwargs=fd_kwargs2)
+xy, uv = dense_lucaskanade(precip, dense=False, fd_kwargs=fd_kwargs2)
 plt.imshow(ref_dbr, cmap=plt.get_cmap("Greys"))
 plt.imshow(mask, cmap=colors.ListedColormap(["black"]), alpha=0.5)
 plt.quiver(
@@ -170,8 +170,8 @@
 # the negative bias that is introduced by the the erroneous interpretation of
 # velocities near the maximum range of the radars.
 
-UV1 = dense_lucaskanade(R, dense=True, fd_kwargs=fd_kwargs1)
-UV2 = dense_lucaskanade(R, dense=True, fd_kwargs=fd_kwargs2)
+UV1 = dense_lucaskanade(precip, dense=True, fd_kwargs=fd_kwargs1)
+UV2 = dense_lucaskanade(precip, dense=True, fd_kwargs=fd_kwargs2)
 
 V1 = np.sqrt(UV1[0] ** 2 + UV1[1] ** 2)
 V2 = np.sqrt(UV2[0] ** 2 + UV2[1] ** 2)
@@ -196,13 +196,13 @@
 # Get the advection routine and extrapolate the last radar frame by 12 time steps
 # (i.e., 1 hour lead time)
 extrapolate = nowcasts.get_method("extrapolation")
-R[~np.isfinite(R)] = metadata["zerovalue"]
-R_f1 = extrapolate(R[-1], UV1, 12)
-R_f2 = extrapolate(R[-1], UV2, 12)
+precip[~np.isfinite(precip)] = metadata["zerovalue"]
+R_f1 = extrapolate(precip[-1], UV1, 12)
+R_f2 = extrapolate(precip[-1], UV2, 12)
 
 # Back-transform to rain rate
-R_f1 = transformation.dB_transform(R_f1, threshold=-10.0, inverse=True)[0]
-R_f2 = transformation.dB_transform(R_f2, threshold=-10.0, inverse=True)[0]
+R_f1 = R_f1.pysteps.db_transform(inverse=True)
+R_f2 = R_f2.pysteps.db_transform(inverse=True)
 
 # Find the veriyfing observations in the archive
 fns = io.archive.find_by_date(
@@ -211,7 +211,7 @@
 
 # Read and convert the radar composites
 R_o, _, metadata_o = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
-R_o, metadata_o = conversion.to_rainrate(R_o, metadata_o)
+R_o = R_o.pysteps.to_rainrate()
 
 # Compute Spearman correlation
 skill = verification.get_method("corr_s")
diff --git a/examples/advection_correction.py b/examples/advection_correction.py
index 706807dbb..9080b1e03 100644
--- a/examples/advection_correction.py
+++ b/examples/advection_correction.py
@@ -7,7 +7,7 @@
 
 Advection correction is a temporal interpolation procedure that is often used
 when estimating rainfall accumulations to correct for the shift of rainfall patterns
-between consecutive radar rainfall maps. This shift becomes particularly 
+between consecutive radar rainfall maps. This shift becomes particularly
 significant for long radar scanning cycles and in presence of fast moving
 precipitation features.
 
@@ -21,7 +21,6 @@
 import numpy as np
 
 from pysteps import io, motion, rcparams
-from pysteps.utils import conversion, dimension
 from pysteps.visualization import plot_precip_field
 from scipy.ndimage import map_coordinates
 
@@ -60,17 +59,18 @@
 
 # Read the radar composites
 importer = io.get_method(importer_name, "importer")
-R, __, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
+precip = io.read_timeseries(fns, importer, **importer_kwargs)
 
 # Convert to mm/h
-R, metadata = conversion.to_rainrate(R, metadata)
+precip = precip.pysteps.to_rainrate()
 
-# Upscale to 2 km (simply to reduce the memory demand)
-R, metadata = dimension.aggregate_fields_space(R, metadata, 2000)
+# Upscale to double the native resolution (simply to reduce the memory demand)
+precip = precip.coarsen(x=2, y=2).mean()
 
 # Keep only one frame every 10 minutes (i.e., every 2 timesteps)
 # (to highlight the need for advection correction)
-R = R[::2]
+precip = precip.isel(t=slice(0, None, 2))
+
 
 ################################################################################
 # Advection correction
@@ -84,9 +84,9 @@
 # going to use the Lucas-Kanade optical flow routine available in pysteps.
 
 
-def advection_correction(R, T=5, t=1):
+def advection_correction(precip, T=5, t=1):
     """
-    R = np.array([qpe_previous, qpe_current])
+    precip = np.array([qpe_previous, qpe_current])
     T = time between two observations (5 min)
     t = interpolation timestep (1 min)
     """
@@ -94,20 +94,20 @@ def advection_correction(R, T=5, t=1):
     # Evaluate advection
     oflow_method = motion.get_method("LK")
     fd_kwargs = {"buffer_mask": 10}  # avoid edge effects
-    V = oflow_method(np.log(R), fd_kwargs=fd_kwargs)
+    V = oflow_method(np.log(precip), fd_kwargs=fd_kwargs)
 
     # Perform temporal interpolation
-    Rd = np.zeros((R[0].shape))
+    Rd = np.zeros((precip[0].shape))
     x, y = np.meshgrid(
-        np.arange(R[0].shape[1], dtype=float), np.arange(R[0].shape[0], dtype=float)
+        np.arange(precip[0].shape[1], dtype=float), np.arange(precip[0].shape[0], dtype=float)
     )
     for i in range(t, T + t, t):
 
         pos1 = (y - i / T * V[1], x - i / T * V[0])
-        R1 = map_coordinates(R[0], pos1, order=1)
+        R1 = map_coordinates(precip[0], pos1, order=1)
 
         pos2 = (y + (T - i) / T * V[1], x + (T - i) / T * V[0])
-        R2 = map_coordinates(R[1], pos2, order=1)
+        R2 = map_coordinates(precip[1], pos2, order=1)
 
         Rd += (T - i) * R1 + i * R2
 
@@ -118,10 +118,10 @@ def advection_correction(R, T=5, t=1):
 # Finally, we apply the advection correction to the whole sequence of radar
 # images and produce the rainfall accumulation map.
 
-R_ac = R[0].copy()
-for i in range(R.shape[0] - 1):
-    R_ac += advection_correction(R[i : (i + 2)], T=10, t=1)
-R_ac /= R.shape[0]
+R_ac = precip[0].copy()
+for i in range(precip.shape[0] - 1):
+    R_ac += advection_correction(precip[i : (i + 2)], T=10, t=1)
+R_ac /= precip.shape[0]
 
 ###############################################################################
 # Results
@@ -136,7 +136,7 @@ def advection_correction(R, T=5, t=1):
 
 plt.figure(figsize=(9, 4))
 plt.subplot(121)
-plot_precip_field(R.mean(axis=0), title="3-h rainfall accumulation")
+plot_precip_field(precip.mean(axis=0), title="3-h rainfall accumulation")
 plt.subplot(122)
 plot_precip_field(R_ac, title="Same with advection correction")
 plt.tight_layout()
diff --git a/examples/anvil_nowcast.py b/examples/anvil_nowcast.py
index fe0f6ad0a..fc3041aa5 100644
--- a/examples/anvil_nowcast.py
+++ b/examples/anvil_nowcast.py
@@ -5,7 +5,7 @@
 =============
 
 This example demonstrates how to use ANVIL and the advantages compared to
-extrapolation nowcast and S-PROG.
+extrapolation nowcast and S-PROG
 
 Load the libraries.
 """
@@ -17,7 +17,6 @@
 import numpy as np
 from pysteps import motion, io, rcparams, utils
 from pysteps.nowcasts import anvil, extrapolation, sprog
-from pysteps.utils import transformation
 from pysteps.visualization import plot_precip_field
 
 ###############################################################################
@@ -47,12 +46,10 @@
 
 # Read the input time series
 importer = io.get_method(importer_name, "importer")
-rainrate_field, quality, metadata = io.read_timeseries(
-    filenames, importer, legacy=True, **importer_kwargs
-)
+rainrate_field = io.read_timeseries(filenames, importer, **importer_kwargs)
 
 # Convert to rain rate (mm/h)
-rainrate_field, metadata = utils.to_rainrate(rainrate_field, metadata)
+rainrate_field = rainrate_field.pysteps.to_rainrate()
 
 ################################################################################
 # Compute the advection field
@@ -78,9 +75,7 @@
 oflow = motion.get_method("lucaskanade")
 
 # transform the input data to logarithmic scale
-rainrate_field_log, _ = utils.transformation.dB_transform(
-    rainrate_field, metadata=metadata
-)
+rainrate_field_log = rainrate_field.pysteps.db_transform()
 velocity = oflow(rainrate_field_log, **oflow_kwargs)
 
 ###############################################################################
@@ -91,25 +86,18 @@
 )
 forecast_extrap[forecast_extrap < 0.5] = 0.0
 
-# log-transform the data and the threshold value to dBR units for S-PROG
-rainrate_field_db, _ = transformation.dB_transform(
-    rainrate_field, metadata, threshold=0.1, zerovalue=-15.0
-)
-rainrate_thr, _ = transformation.dB_transform(
-    np.array([0.5]), metadata, threshold=0.1, zerovalue=-15.0
-)
+# log-transform the data to dBR units for S-PROG
+rainrate_field_db = rainrate_field.pysteps.db_transform()
 forecast_sprog = sprog.forecast(
     rainrate_field_db[-3:], velocity, 3, n_cascade_levels=8, R_thr=rainrate_thr[0]
 )
-forecast_sprog, _ = transformation.dB_transform(
-    forecast_sprog, threshold=-10.0, inverse=True
-)
-forecast_sprog[forecast_sprog < 0.5] = 0.0
+forecast_sprog = forecast_sprog.pysteps.db_transform(inverse=True)
+forecast_sprog = forecast_sprog.where(forecast_sprog >= 0.5, 0)
 
 forecast_anvil = anvil.forecast(
     rainrate_field[-4:], velocity, 3, ar_window_radius=25, ar_order=2
 )
-forecast_anvil[forecast_anvil < 0.5] = 0.0
+forecast_anvil = forecast_anvil.where(forecast_sprog >= 0.5, 0)
 
 ###############################################################################
 # Read the reference observation field and threshold rain rates below 0.5 mm/h
@@ -122,8 +110,8 @@
     filenames, importer, legacy=True, **importer_kwargs
 )
 
-refobs_field, metadata = utils.to_rainrate(refobs_field[-1], metadata)
-refobs_field[refobs_field < 0.5] = 0.0
+refobs_field = refobs_field.isel(t=-1).pysteps.to_rainrate()
+refobs_field = refobs_field.where(refobs_field >= 0.5, 0.0)
 
 ###############################################################################
 # Plot the extrapolation, S-PROG and ANVIL nowcasts.
diff --git a/examples/plot_cascade_decomposition.py b/examples/cascade_decomposition.py
similarity index 60%
rename from examples/plot_cascade_decomposition.py
rename to examples/cascade_decomposition.py
index 49a3798d8..9c3aa8ed5 100644
--- a/examples/plot_cascade_decomposition.py
+++ b/examples/cascade_decomposition.py
@@ -3,7 +3,7 @@
 Cascade decomposition
 =====================
 
-This example script shows how to compute and plot the cascade decompositon of 
+This example script shows how to compute and plot the cascade decomposition of
 a single radar precipitation field in pysteps.
 
 """
@@ -15,7 +15,6 @@
 from pysteps.cascade.bandpass_filters import filter_gaussian
 from pysteps import io, rcparams
 from pysteps.cascade.decomposition import decomposition_fft
-from pysteps.utils import conversion, transformation
 from pysteps.visualization import plot_precip_field
 
 ###############################################################################
@@ -30,43 +29,48 @@
 filename = os.path.join(
     root_path, "20160928", "201609281600_fmi.radar.composite.lowest_FIN_SUOMI1.pgm.gz"
 )
-R, _, metadata = io.import_fmi_pgm(filename, gzipped=True, legacy=True)
+precip = io.import_fmi_pgm(filename, gzipped=True)
 
 # Convert to rain rate
-R, metadata = conversion.to_rainrate(R, metadata)
+precip = precip.pysteps.to_rainrate()
 
-# Nicely print the metadata
-pprint(metadata)
+# Print the data array
+print(precip)
 
 # Plot the rainfall field
-plot_precip_field(R, geodata=metadata)
+plot_precip_field(precip)
 plt.show()
 
+
 # Log-transform the data
-R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)
+precip = precip.pysteps.db_transform()
 
 ###############################################################################
 # 2D Fourier spectrum
 # --------------------
 #
-# Compute and plot the 2D Fourier power spectrum of the precipitaton field.
+# Compute and plot the 2D Fourier power spectrum of the precipitation field.
 
 # Set Nans as the fill value
-R[~np.isfinite(R)] = metadata["zerovalue"]
+precip = precip.fillna(precip.min())
 
 # Compute the Fourier transform of the input field
-F = abs(np.fft.fftshift(np.fft.fft2(R)))
+power_spectrum = abs(np.fft.fftshift(np.fft.fft2(precip)))
 
 # Plot the power spectrum
-M, N = F.shape
+rows, cols = power_spectrum.shape
 fig, ax = plt.subplots()
 im = ax.imshow(
-    np.log(F ** 2), vmin=4, vmax=24, cmap=cm.jet, extent=(-N / 2, N / 2, -M / 2, M / 2)
+    np.log(power_spectrum ** 2),
+    vmin=4,
+    vmax=24,
+    cmap=cm.jet,
+    extent=(-cols / 2, cols / 2, -rows / 2, rows / 2)
 )
 cb = fig.colorbar(im)
 ax.set_xlabel("Wavenumber $k_x$")
 ax.set_ylabel("Wavenumber $k_y$")
-ax.set_title("Log-power spectrum of R")
+ax.set_title("Log-power spectrum of precip")
 plt.show()
 
 ###############################################################################
@@ -79,23 +83,23 @@
 num_cascade_levels = 7
 
 # Construct the Gaussian bandpass filters
-filter = filter_gaussian(R.shape, num_cascade_levels)
+filter_object = filter_gaussian(precip.shape, num_cascade_levels)
 
 # Plot the bandpass filter weights
-L = max(N, M)
+size = max(cols, rows)
 fig, ax = plt.subplots()
 for k in range(num_cascade_levels):
     ax.semilogx(
-        np.linspace(0, L / 2, len(filter["weights_1d"][k, :])),
-        filter["weights_1d"][k, :],
+        np.linspace(0, size / 2, len(filter_object["weights_1d"][k, :])),
+        filter_object["weights_1d"][k, :],
         "k-",
-        basex=pow(0.5 * L / 3, 1.0 / (num_cascade_levels - 2)),
+        basex=pow(0.5 * size / 3, 1.0 / (num_cascade_levels - 2)),
     )
-ax.set_xlim(1, L / 2)
+ax.set_xlim(1, size / 2)
 ax.set_ylim(0, 1)
-xt = np.hstack([[1.0], filter["central_wavenumbers"][1:]])
+xt = np.hstack([[1.0], filter_object["central_wavenumbers"][1:]])
 ax.set_xticks(xt)
-ax.set_xticklabels(["%.2f" % cf for cf in filter["central_wavenumbers"]])
+ax.set_xticklabels(["%.2f" % cf for cf in filter_object["central_wavenumbers"]])
 ax.set_xlabel("Radial wavenumber $|\mathbf{k}|$")
 ax.set_ylabel("Normalized weight")
 ax.set_title("Bandpass filter weights")
@@ -105,24 +109,24 @@
 # Finally, apply the 2D Gaussian filters to decompose the radar rainfall field
 # into a set of cascade levels of decreasing spatial scale and plot them.
 
-decomp = decomposition_fft(R, filter, compute_stats=True)
+precip_decomp = decomposition_fft(precip, filter_object, compute_stats=True)
 
 # Plot the normalized cascade levels
 for i in range(num_cascade_levels):
-    mu = decomp["means"][i]
-    sigma = decomp["stds"][i]
-    decomp["cascade_levels"][i] = (decomp["cascade_levels"][i] - mu) / sigma
+    mu = precip_decomp["means"][i]
+    sigma = precip_decomp["stds"][i]
+    precip_decomp["cascade_levels"][i] = (precip_decomp["cascade_levels"][i] - mu) / sigma
 
 fig, ax = plt.subplots(nrows=2, ncols=4)
 
-ax[0, 0].imshow(R, cmap=cm.RdBu_r, vmin=-5, vmax=5)
-ax[0, 1].imshow(decomp["cascade_levels"][0], cmap=cm.RdBu_r, vmin=-3, vmax=3)
-ax[0, 2].imshow(decomp["cascade_levels"][1], cmap=cm.RdBu_r, vmin=-3, vmax=3)
-ax[0, 3].imshow(decomp["cascade_levels"][2], cmap=cm.RdBu_r, vmin=-3, vmax=3)
-ax[1, 0].imshow(decomp["cascade_levels"][3], cmap=cm.RdBu_r, vmin=-3, vmax=3)
-ax[1, 1].imshow(decomp["cascade_levels"][4], cmap=cm.RdBu_r, vmin=-3, vmax=3)
-ax[1, 2].imshow(decomp["cascade_levels"][5], cmap=cm.RdBu_r, vmin=-3, vmax=3)
-ax[1, 3].imshow(decomp["cascade_levels"][6], cmap=cm.RdBu_r, vmin=-3, vmax=3)
+ax[0, 0].imshow(precip, cmap=cm.RdBu_r, vmin=-5, vmax=5)
+ax[0, 1].imshow(precip_decomp["cascade_levels"][0], cmap=cm.RdBu_r, vmin=-3, vmax=3)
+ax[0, 2].imshow(precip_decomp["cascade_levels"][1], cmap=cm.RdBu_r, vmin=-3, vmax=3)
+ax[0, 3].imshow(precip_decomp["cascade_levels"][2], cmap=cm.RdBu_r, vmin=-3, vmax=3)
+ax[1, 0].imshow(precip_decomp["cascade_levels"][3], cmap=cm.RdBu_r, vmin=-3, vmax=3)
+ax[1, 1].imshow(precip_decomp["cascade_levels"][4], cmap=cm.RdBu_r, vmin=-3, vmax=3)
+ax[1, 2].imshow(precip_decomp["cascade_levels"][5], cmap=cm.RdBu_r, vmin=-3, vmax=3)
+ax[1, 3].imshow(precip_decomp["cascade_levels"][6], cmap=cm.RdBu_r, vmin=-3, vmax=3)
 
 ax[0, 0].set_title("Observed")
 ax[0, 1].set_title("Level 1")
diff --git a/examples/data_transformations.py b/examples/data_transformations.py
index 8b601ac14..de6e04da8 100644
--- a/examples/data_transformations.py
+++ b/examples/data_transformations.py
@@ -23,7 +23,6 @@
 import matplotlib.pyplot as plt
 import numpy as np
 from pysteps import io, rcparams
-from pysteps.utils import conversion, transformation
 from scipy.stats import skew
 
 ###############################################################################
@@ -58,13 +57,13 @@
 
 # Read the radar composites
 importer = io.get_method(importer_name, "importer")
-Z, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
-
-# Keep only positive rainfall values
-Z = Z[Z > metadata["zerovalue"]].flatten()
+precip = io.read_timeseries(fns, importer, **importer_kwargs)
 
 # Convert to rain rate
-R, metadata = conversion.to_rainrate(Z, metadata)
+precip = precip.pysteps.to_rainrate()
+
+# Keep only positive rainfall values
+precip = precip.where(precip > 0)
 
 ###############################################################################
 # Test data transformations
@@ -74,15 +73,20 @@
 # corresponding skewness
 def plot_distribution(data, labels, skw):
 
-    N = len(data)
-    fig, ax1 = plt.subplots()
+    # Extract and flatten the arrays
+    data = [da.values.flatten() for da in data]
+    # Keep finite values only
+    data = [da[np.isfinite(da)] for da in data]
+
+    n_points = len(data)
+    _, ax1 = plt.subplots()
     ax2 = ax1.twinx()
 
-    ax2.plot(np.arange(N + 2), np.zeros(N + 2), ":r")
+    ax2.plot(np.arange(n_points + 2), np.zeros(n_points + 2), ":r")
     ax1.boxplot(data, labels=labels, sym="", medianprops={"color": "k"})
 
     ymax = []
-    for i in range(N):
+    for i in range(n_points):
         y = skw[i]
         x = i + 1
         ax2.plot(x, y, "*r", ms=10, markeredgecolor="k")
@@ -120,13 +124,13 @@ def plot_distribution(data, labels, skw):
 skw = []
 
 # Test a range of values for the transformation parameter Lambda
-Lambdas = np.linspace(-0.4, 0.4, 11)
-for i, Lambda in enumerate(Lambdas):
-    R_, _ = transformation.boxcox_transform(R, metadata, Lambda)
-    R_ = (R_ - np.mean(R_)) / np.std(R_)
-    data.append(R_)
-    labels.append("{0:.2f}".format(Lambda))
-    skw.append(skew(R_))  # skewness
+lambdas = np.linspace(-0.4, 0.4, 11)
+for i, this_lambda in enumerate(lambdas):
+    precip_ = precip.pysteps.boxcox_transform(boxcox_lambda=this_lambda)
+    precip_ = (precip_ - np.mean(precip_)) / np.std(precip_)
+    data.append(precip_)
+    labels.append("{0:.2f}".format(this_lambda))
+    skw.append(skew(precip_, axis=None, nan_policy="omit"))  # skewness
 
 # Plot the transformed data distribution as a function of lambda
 plot_distribution(data, labels, skw)
@@ -136,9 +140,9 @@ def plot_distribution(data, labels, skw):
 
 # Best lambda
 idx_best = np.argmin(np.abs(skw))
-Lambda = Lambdas[idx_best]
+best_lambda = lambdas[idx_best]
 
-print("Best parameter lambda: %.2f\n(skewness = %.2f)" % (Lambda, skw[idx_best]))
+print("Best parameter lambda: %.2f\n(skewness = %.2f)" % (best_lambda, skw[idx_best]))
 
 ###############################################################################
 # Compare data transformations
@@ -153,39 +157,39 @@ def plot_distribution(data, labels, skw):
 # ~~~~~~~~~~
 # First, let's have a look at the original rain rate values.
 
-data.append((R - np.mean(R)) / np.std(R))
-labels.append("R")
-skw.append(skew(R))
+data.append((precip - np.mean(precip)) / np.std(precip))
+labels.append("precip")
+skw.append(skew(precip, axis=None, nan_policy="omit"))
 
 ###############################################################################
 # dB transform
 # ~~~~~~~~~~~~
-# We transform the rainfall data into dB units: 10*log(R)
+# We transform the rainfall data into dB units: 10*log(precip)
 
-R_, _ = transformation.dB_transform(R, metadata)
-data.append((R_ - np.mean(R_)) / np.std(R_))
+precip_ = precip.pysteps.db_transform()
+data.append((precip_ - np.mean(precip_)) / np.std(precip_))
 labels.append("dB")
-skw.append(skew(R_))
+skw.append(skew(precip_, axis=None, nan_policy="omit"))
 
 ###############################################################################
 # Square-root transform
 # ~~~~~~~~~~~~~~~~~~~~~
-# Transform the data using the square-root: sqrt(R)
+# Transform the data using the square-root: sqrt(precip)
 
-R_, _ = transformation.sqrt_transform(R, metadata)
-data.append((R_ - np.mean(R_)) / np.std(R_))
+precip_ = precip.pysteps.sqrt_transform()
+data.append((precip_ - np.mean(precip_)) / np.std(precip_))
 labels.append("sqrt")
-skw.append(skew(R_))
+skw.append(skew(precip_, axis=None, nan_policy="omit"))
 
 ###############################################################################
 # Box-Cox transform
 # ~~~~~~~~~~~~~~~~~
 # We now apply the Box-Cox transform using the best parameter lambda found above.
 
-R_, _ = transformation.boxcox_transform(R, metadata, Lambda)
-data.append((R_ - np.mean(R_)) / np.std(R_))
-labels.append("Box-Cox\n($\lambda=$%.2f)" % Lambda)
-skw.append(skew(R_))
+precip_ = precip.pysteps.boxcox_transform(best_lambda)
+data.append((precip_ - np.mean(precip_)) / np.std(precip_))
+labels.append("Box-Cox\n($\lambda=$%.2f)" % best_lambda)
+skw.append(skew(precip_, axis=None, nan_policy="omit"))
 
 ###############################################################################
 # Normal quantile transform
@@ -195,14 +199,14 @@ def plot_distribution(data, labels, skw):
 #
 # .. _`Bogner et al (2012)`: http://dx.doi.org/10.5194/hess-16-1085-2012
 
-R_, _ = transformation.NQ_transform(R, metadata)
-data.append((R_ - np.mean(R_)) / np.std(R_))
+precip_ = precip.pysteps.nq_transform()
+data.append((precip_ - np.mean(precip_)) / np.std(precip_))
 labels.append("NQ")
-skw.append(skew(R_))
+skw.append(skew(precip_, axis=None, nan_policy="omit"))
 
 ###############################################################################
 # By plotting all the results, we can notice first of all the strongly asymmetric
-# distribution of the original data (R) and that all transformations manage to
+# distribution of the original data (precip) and that all transformations manage to
 # reduce its skewness. Among these, the Box-Cox transform (using the best parameter
 # lambda) and the normal quantile (NQ) transform provide the best correction.
 # Despite not producing a perfectly symmetric distribution, the square-root (sqrt)
diff --git a/examples/linda_nowcasts.py b/examples/linda_nowcasts.py
index 991bade6c..e551686f6 100644
--- a/examples/linda_nowcasts.py
+++ b/examples/linda_nowcasts.py
@@ -83,9 +83,7 @@
 )[0]
 
 # Compute S-PROG nowcast for comparison
-rainrate_db, _ = transformation.dB_transform(
-    rainrate, metadata, threshold=0.1, zerovalue=-15.0
-)
+rainrate_db = rainrate.pysteps.db_transform()
 nowcast_sprog = sprog.forecast(
     rainrate_db[-3:, :, :],
     advection,
@@ -95,9 +93,7 @@
 )
 
 # Convert reflectivity nowcast to rain rate
-nowcast_sprog = transformation.dB_transform(
-    nowcast_sprog, threshold=-10.0, inverse=True
-)[0]
+nowcast_sprog = nowcast_sprog.pysteps.db_transform(inverse=True)
 
 # Plot the nowcasts
 fig = plt.figure(figsize=(9, 4))
@@ -153,9 +149,7 @@
 )
 
 # Convert reflectivity nowcast to rain rate
-nowcast_steps = transformation.dB_transform(
-    nowcast_steps, threshold=-10.0, inverse=True
-)[0]
+nowcast_steps = nowcast_steps.pysteps.db_transform(inverse=True)
 
 # Plot two ensemble members of both nowcasts
 fig = plt.figure()
diff --git a/examples/my_first_nowcast.ipynb b/examples/my_first_nowcast_pystepsv1.ipynb
similarity index 99%
rename from examples/my_first_nowcast.ipynb
rename to examples/my_first_nowcast_pystepsv1.ipynb
index 3113c4977..945e75ad4 100644
--- a/examples/my_first_nowcast.ipynb
+++ b/examples/my_first_nowcast_pystepsv1.ipynb
@@ -509,7 +509,7 @@
     "\n",
     "# Log-transform the data to dBR. \n",
     "# The threshold of 0.1 mm/h sets the fill value to -15 dBR.\n",
-    "train_precip_dbr, metadata_dbr = transformation.dB_transform(train_precip, metadata, \n",
+    "train_precip_dbr, metadata_dbr = transformation.db_transform(train_precip, metadata, \n",
     "                                                             threshold=0.1, \n",
     "                                                             zerovalue=-15.0)"
    ]
diff --git a/examples/plot_noise_generators.py b/examples/noise_generators.py
similarity index 66%
rename from examples/plot_noise_generators.py
rename to examples/noise_generators.py
index a809eefa1..dfc4b9b96 100644
--- a/examples/plot_noise_generators.py
+++ b/examples/noise_generators.py
@@ -14,12 +14,11 @@
 from matplotlib import cm, pyplot as plt
 import numpy as np
 import os
-from pprint import pprint
 from pysteps import io, rcparams
 from pysteps.noise.fftgenerators import initialize_param_2d_fft_filter
 from pysteps.noise.fftgenerators import initialize_nonparam_2d_fft_filter
 from pysteps.noise.fftgenerators import generate_noise_2d_fft_filter
-from pysteps.utils import conversion, rapsd, transformation
+from pysteps.utils import rapsd
 from pysteps.visualization import plot_precip_field, plot_spectrum1d
 
 ###############################################################################
@@ -34,25 +33,23 @@
 # Import the example radar composite
 root_path = rcparams.data_sources["mch"]["root_path"]
 filename = os.path.join(root_path, "20160711", "AQC161932100V_00005.801.gif")
-R, _, metadata = io.import_mch_gif(
-    filename, product="AQC", unit="mm", accutime=5.0, legacy=True
-)
+precip = io.import_mch_gif(filename, product="AQC", unit="mm", accutime=5.0)
 
 # Convert to mm/h
-R, metadata = conversion.to_rainrate(R, metadata)
+precip = precip.pysteps.to_rainrate()
 
 # Nicely print the metadata
-pprint(metadata)
+print(precip)
 
 # Plot the rainfall field
-plot_precip_field(R, geodata=metadata)
+plot_precip_field(precip)
 plt.show()
 
 # Log-transform the data
-R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)
+precip = precip.pysteps.db_transform()
 
-# Assign the fill value to all the Nans
-R[~np.isfinite(R)] = metadata["zerovalue"]
+# Assign a fill value to all the Nans
+precip = precip.fillna(precip.min())
 
 ###############################################################################
 # Parametric filter
@@ -65,28 +62,30 @@
 # slopes (beta1 and beta2) and one breaking point
 
 # Fit the parametric PSD to the observation
-Fp = initialize_param_2d_fft_filter(R)
+filter_parametric = initialize_param_2d_fft_filter(precip.values)
 
 # Compute the observed and fitted 1D PSD
-L = np.max(Fp["input_shape"])
-if L % 2 == 1:
-    wn = np.arange(0, int(L / 2) + 1)
+size = np.max(filter_parametric["input_shape"])
+if size % 2 == 1:
+    wn = np.arange(0, int(size / 2) + 1)
 else:
-    wn = np.arange(0, int(L / 2))
-R_, freq = rapsd(R, fft_method=np.fft, return_freq=True)
-f = np.exp(Fp["model"](np.log(wn), *Fp["pars"]))
+    wn = np.arange(0, int(size / 2))
+spectrum_1d, freq = rapsd(precip, fft_method=np.fft, return_freq=True)
+spectrum_1d_param = np.exp(
+    filter_parametric["model"](np.log(wn), *filter_parametric["pars"])
+)
 
 # Extract the scaling break in km, beta1 and beta2
-w0 = L / np.exp(Fp["pars"][0])
-b1 = Fp["pars"][2]
-b2 = Fp["pars"][3]
+scaling_break = size / np.exp(filter_parametric["pars"][0])
+beta_1 = filter_parametric["pars"][2]
+beta_2 = filter_parametric["pars"][3]
 
 # Plot the observed power spectrum and the model
 fig, ax = plt.subplots()
 plot_scales = [512, 256, 128, 64, 32, 16, 8, 4]
 plot_spectrum1d(
     freq,
-    R_,
+    spectrum_1d,
     x_units="km",
     y_units="dBR",
     color="k",
@@ -96,7 +95,7 @@
 )
 plot_spectrum1d(
     freq,
-    f,
+    spectrum_1d_param,
     x_units="km",
     y_units="dBR",
     color="r",
@@ -106,8 +105,8 @@
 )
 plt.legend()
 ax.set_title(
-    "Radially averaged log-power spectrum of R\n"
-    r"$\omega_0=%.0f km, \beta_1=%.1f, \beta_2=%.1f$" % (w0, b1, b2)
+    "Radially averaged log-power spectrum of precip\n"
+    f"scale break={scaling_break:.0f} km, beta_1={beta_1:.1f}, beta_2={beta_2:.1f}"
 )
 plt.show()
 
@@ -116,9 +115,9 @@
 # --------------------
 #
 # In the nonparametric approach,  the Fourier filter is obtained directly
-# from the power spectrum of the observed precipitation field R.
+# from the power spectrum of the observed precipitation field precip.
 
-Fnp = initialize_nonparam_2d_fft_filter(R)
+filter_nonparametric = initialize_nonparam_2d_fft_filter(precip.values)
 
 ###############################################################################
 # Noise generator
@@ -131,26 +130,29 @@
 seed = 42
 num_realizations = 3
 
-# Generate noise
-Np = []
+# Generate noise = []
+samples_param = []
+samples_nonparam = []
 Nnp = []
 for k in range(num_realizations):
-    Np.append(generate_noise_2d_fft_filter(Fp, seed=seed + k))
-    Nnp.append(generate_noise_2d_fft_filter(Fnp, seed=seed + k))
+    samples_param.append(generate_noise_2d_fft_filter(filter_parametric, seed=seed + k))
+    samples_nonparam.append(
+        generate_noise_2d_fft_filter(filter_nonparametric, seed=seed + k)
+    )
 
 # Plot the generated noise fields
 
 fig, ax = plt.subplots(nrows=2, ncols=3)
 
 # parametric noise
-ax[0, 0].imshow(Np[0], cmap=cm.RdBu_r, vmin=-3, vmax=3)
-ax[0, 1].imshow(Np[1], cmap=cm.RdBu_r, vmin=-3, vmax=3)
-ax[0, 2].imshow(Np[2], cmap=cm.RdBu_r, vmin=-3, vmax=3)
+ax[0, 0].imshow(samples_param[0], cmap=cm.RdBu_r, vmin=-3, vmax=3)
+ax[0, 1].imshow(samples_param[1], cmap=cm.RdBu_r, vmin=-3, vmax=3)
+ax[0, 2].imshow(samples_param[2], cmap=cm.RdBu_r, vmin=-3, vmax=3)
 
 # nonparametric noise
-ax[1, 0].imshow(Nnp[0], cmap=cm.RdBu_r, vmin=-3, vmax=3)
-ax[1, 1].imshow(Nnp[1], cmap=cm.RdBu_r, vmin=-3, vmax=3)
-ax[1, 2].imshow(Nnp[2], cmap=cm.RdBu_r, vmin=-3, vmax=3)
+ax[1, 0].imshow(samples_nonparam[0], cmap=cm.RdBu_r, vmin=-3, vmax=3)
+ax[1, 1].imshow(samples_nonparam[1], cmap=cm.RdBu_r, vmin=-3, vmax=3)
+ax[1, 2].imshow(samples_nonparam[2], cmap=cm.RdBu_r, vmin=-3, vmax=3)
 
 for i in range(2):
     for j in range(3):
diff --git a/examples/optical_flow_methods_convergence.py b/examples/optical_flow_methods_convergence.py
index f6289450a..f3501b30f 100644
--- a/examples/optical_flow_methods_convergence.py
+++ b/examples/optical_flow_methods_convergence.py
@@ -74,7 +74,7 @@
 # ~~~~~~~~~~~~~~~~~~~
 
 # Convert to mm/h
-reference_field, metadata = stp.utils.to_rainrate(reference_field, metadata)
+reference_field = reference_field.pysteps.to_rainrate()
 
 # Mask invalid values
 reference_field = np.ma.masked_invalid(reference_field)
@@ -84,9 +84,7 @@
 plt.show()
 
 # Log-transform the data [dBR]
-reference_field, metadata = stp.utils.dB_transform(
-    reference_field, metadata, threshold=0.1, zerovalue=-15.0
-)
+reference_field = reference_field.pysteps.db_transform()
 
 print("Precip. pattern shape: " + str(reference_field.shape))
 
@@ -281,7 +279,7 @@ def plot_optflow_method_convergence(input_precip, optflow_method_name, motion_ty
         f"{(time.perf_counter() - elapsed_time):.1f} [s]"
     )
 
-    precip_obs, _ = stp.utils.dB_transform(precip_obs, inverse=True)
+    precip_obs = precip_obs.pysteps.db_transform(inverse=True)
 
     precip_data = precip_obs.max(axis=0)
     precip_data.data[precip_data.mask] = 0
diff --git a/examples/plot_ensemble_verification.py b/examples/plot_ensemble_verification.py
index 799c43899..0e0e9cb3e 100644
--- a/examples/plot_ensemble_verification.py
+++ b/examples/plot_ensemble_verification.py
@@ -3,7 +3,7 @@
 Ensemble verification
 =====================
 
-In this tutorial we perform a verification of a probabilistic extrapolation nowcast 
+In this tutorial we perform a verification of a probabilistic extrapolation nowcast
 using MeteoSwiss radar data.
 
 """
@@ -15,7 +15,7 @@
 from pysteps import io, nowcasts, rcparams, verification
 from pysteps.motion.lucaskanade import dense_lucaskanade
 from pysteps.postprocessing import ensemblestats
-from pysteps.utils import conversion, dimension, transformation
+from pysteps.utils import conversion, transformation
 from pysteps.visualization import plot_precip_field
 
 
@@ -56,24 +56,23 @@
 
 # Read the data from the archive
 importer = io.get_method(importer_name, "importer")
-R, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
+precip, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
 
 # Convert to rain rate
-R, metadata = conversion.to_rainrate(R, metadata)
+precip = precip.pysteps.to_rainrate()
 
 # Upscale data to 2 km
-R, metadata = dimension.aggregate_fields_space(R, metadata, 2000)
+precip = precip.coarsen(x=2, y=2).mean()
 
 # Plot the rainfall field
-plot_precip_field(R[-1, :, :], geodata=metadata)
+plot_precip_field(precip[-1, :, :], geodata=metadata)
 plt.show()
 
-# Log-transform the data to unit of dBR, set the threshold to 0.1 mm/h,
-# set the fill value to -15 dBR
-R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)
+# Log-transform the data to unit of dBR
+precip = precip.pysteps.db_transform()
 
 # Set missing values with the fill value
-R[~np.isfinite(R)] = -15.0
+precip[~np.isfinite(precip)] = -15.0
 
 # Nicely print the metadata
 pprint(metadata)
@@ -85,12 +84,12 @@
 # We use the STEPS approach to produce a ensemble nowcast of precipitation fields.
 
 # Estimate the motion field
-V = dense_lucaskanade(R)
+V = dense_lucaskanade(precip)
 
 # Perform the ensemble nowcast with STEPS
 nowcast_method = nowcasts.get_method("steps")
 R_f = nowcast_method(
-    R[-3:, :, :],
+    precip[-3:, :, :],
     V,
     n_leadtimes,
     n_ens_members,
@@ -107,7 +106,7 @@
 )
 
 # Back-transform to rain rates
-R_f = transformation.dB_transform(R_f, threshold=-10.0, inverse=True)[0]
+R_f = R_f.pysteps.db_transform(inverse=True)
 
 # Plot some of the realizations
 fig = plt.figure()
@@ -145,10 +144,10 @@
 R_o, _, metadata_o = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
 
 # Convert to mm/h
-R_o, metadata_o = conversion.to_rainrate(R_o, metadata_o)
+R_o = R_o.pysteps.to_rainrate()
 
 # Upscale data to 2 km
-R_o, metadata_o = dimension.aggregate_fields_space(R_o, metadata_o, 2000)
+R_o = R_o.coarsen(x=2, y=2).mean()
 
 # Compute the verification for the last lead time
 
diff --git a/examples/plot_extrapolation_nowcast.py b/examples/plot_extrapolation_nowcast.py
index 2e76eeeef..f9f28f3d5 100644
--- a/examples/plot_extrapolation_nowcast.py
+++ b/examples/plot_extrapolation_nowcast.py
@@ -3,7 +3,7 @@
 Extrapolation nowcast
 =====================
 
-This tutorial shows how to compute and plot an extrapolation nowcast using 
+This tutorial shows how to compute and plot an extrapolation nowcast using
 Finnish radar data.
 
 """
@@ -48,21 +48,20 @@
 
 # Read the radar composites
 importer = io.get_method(importer_name, "importer")
-Z, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
+reflectivity, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
 
 # Convert to rain rate
-R, metadata = conversion.to_rainrate(Z, metadata)
+precip = reflectivity.pysteps.to_rainrate()
 
 # Plot the rainfall field
-plot_precip_field(R[-1, :, :], geodata=metadata)
+plot_precip_field(precip[-1, :, :], geodata=metadata)
 plt.show()
 
 # Store the last frame for plotting it later later
-R_ = R[-1, :, :].copy()
+R_ = precip[-1, :, :].copy()
 
-# Log-transform the data to unit of dBR, set the threshold to 0.1 mm/h,
-# set the fill value to -15 dBR
-R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)
+# Log-transform the data to unit of dBR
+precip = precip.pysteps.db_transform()
 
 # Nicely print the metadata
 pprint(metadata)
@@ -78,15 +77,15 @@
 
 # Estimate the motion field with Lucas-Kanade
 oflow_method = motion.get_method("LK")
-V = oflow_method(R[-3:, :, :])
+V = oflow_method(precip[-3:, :, :])
 
 # Extrapolate the last radar observation
 extrapolate = nowcasts.get_method("extrapolation")
-R[~np.isfinite(R)] = metadata["zerovalue"]
-R_f = extrapolate(R[-1, :, :], V, n_leadtimes)
+precip[~np.isfinite(precip)] = metadata["zerovalue"]
+R_f = extrapolate(precip[-1, :, :], V, n_leadtimes)
 
 # Back-transform to rain rate
-R_f = transformation.dB_transform(R_f, threshold=-10.0, inverse=True)[0]
+R_f = R_f.pysteps.db_transform(inverse=True)
 
 # Plot the motion field
 plot_precip_field(R_, geodata=metadata)
@@ -114,7 +113,7 @@
 )
 # Read the radar composites
 R_o, _, metadata_o = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
-R_o, metadata_o = conversion.to_rainrate(R_o, metadata_o, 223.0, 1.53)
+R_o = R_o.pysteps.to_rainrate(zr_a=223.0, zr_b=1.53)
 
 # Compute fractions skill score (FSS) for all lead times, a set of scales and 1 mm/h
 fss = verification.get_method("FSS")
diff --git a/examples/plot_optical_flow.py b/examples/plot_optical_flow.py
index 89ed3b37e..0050d09c6 100644
--- a/examples/plot_optical_flow.py
+++ b/examples/plot_optical_flow.py
@@ -2,8 +2,8 @@
 Optical flow
 ============
 
-This tutorial offers a short overview of the optical flow routines available in 
-pysteps and it will cover how to compute and plot the motion field from a 
+This tutorial offers a short overview of the optical flow routines available in
+pysteps and it will cover how to compute and plot the motion field from a
 sequence of radar images.
 """
 
@@ -47,7 +47,7 @@
 
 # Read the radar composites
 importer = io.get_method(importer_name, "importer")
-R, quality, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
+precip, quality, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
 
 del quality  # Not used
 
@@ -56,13 +56,13 @@
 # ~~~~~~~~~~~~~~~~~~~
 
 # Convert to mm/h
-R, metadata = conversion.to_rainrate(R, metadata)
+precip = precip.pysteps.to_rainrate()
 
 # Store the reference frame
-R_ = R[-1, :, :].copy()
+R_ = precip[-1, :, :].copy()
 
 # Log-transform the data [dBR]
-R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)
+precip = precip.pysteps.db_transform()
 
 # Nicely print the metadata
 pprint(metadata)
@@ -78,7 +78,7 @@
 # field of motion vectors.
 
 oflow_method = motion.get_method("LK")
-V1 = oflow_method(R[-3:, :, :])
+V1 = oflow_method(precip[-3:, :, :])
 
 # Plot the motion field on top of the reference frame
 plot_precip_field(R_, geodata=metadata, title="LK")
@@ -97,7 +97,7 @@
 # at minimizing a cost function between the displaced and the reference image.
 
 oflow_method = motion.get_method("VET")
-V2 = oflow_method(R[-3:, :, :])
+V2 = oflow_method(precip[-3:, :, :])
 
 # Plot the motion field
 plot_precip_field(R_, geodata=metadata, title="VET")
@@ -116,8 +116,8 @@
 # estimating the motion, here we are going to use all the available 10 fields.
 
 oflow_method = motion.get_method("DARTS")
-R[~np.isfinite(R)] = metadata["zerovalue"]
-V3 = oflow_method(R)  # needs longer training sequence
+precip[~np.isfinite(precip)] = metadata["zerovalue"]
+V3 = oflow_method(precip)  # needs longer training sequence
 
 # Plot the motion field
 plot_precip_field(R_, geodata=metadata, title="DARTS")
@@ -133,8 +133,8 @@
 # inconsitency during the solution of the optical flow equations.
 
 oflow_method = motion.get_method("proesmans")
-R[~np.isfinite(R)] = metadata["zerovalue"]
-V4 = oflow_method(R[-2:, :, :])
+precip[~np.isfinite(precip)] = metadata["zerovalue"]
+V4 = oflow_method(precip[-2:, :, :])
 
 # Plot the motion field
 plot_precip_field(R_, geodata=metadata, title="Proesmans")
diff --git a/examples/plot_steps_nowcast.py b/examples/plot_steps_nowcast.py
index 38139b848..7d48639cc 100644
--- a/examples/plot_steps_nowcast.py
+++ b/examples/plot_steps_nowcast.py
@@ -16,7 +16,7 @@
 from pysteps import io, nowcasts, rcparams
 from pysteps.motion.lucaskanade import dense_lucaskanade
 from pysteps.postprocessing.ensemblestats import excprob
-from pysteps.utils import conversion, dimension, transformation
+from pysteps.utils import conversion, transformation
 from pysteps.visualization import plot_precip_field
 
 # Set nowcast parameters
@@ -51,24 +51,22 @@
 
 # Read the data from the archive
 importer = io.get_method(importer_name, "importer")
-R, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
+precip, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
 
 # Convert to rain rate
-R, metadata = conversion.to_rainrate(R, metadata)
+precip = precip.pysteps.to_rainrate()
 
 # Upscale data to 2 km to limit memory usage
-R, metadata = dimension.aggregate_fields_space(R, metadata, 2000)
-
+precip = precip.coarsen(x=2, y=2).mean()
 # Plot the rainfall field
-plot_precip_field(R[-1, :, :], geodata=metadata)
+plot_precip_field(precip[-1, :, :], geodata=metadata)
 plt.show()
 
-# Log-transform the data to unit of dBR, set the threshold to 0.1 mm/h,
-# set the fill value to -15 dBR
-R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)
+# Log-transform the data to unit of dBR
+precip = precip.pysteps.db_transform()
 
 # Set missing values with the fill value
-R[~np.isfinite(R)] = -15.0
+precip[~np.isfinite(precip)] = -15.0
 
 # Nicely print the metadata
 pprint(metadata)
@@ -84,12 +82,12 @@
 # progressively remove the unpredictable spatial scales during the forecast.
 
 # Estimate the motion field
-V = dense_lucaskanade(R)
+V = dense_lucaskanade(precip)
 
 # The S-PROG nowcast
 nowcast_method = nowcasts.get_method("sprog")
 R_f = nowcast_method(
-    R[-3:, :, :],
+    precip[-3:, :, :],
     V,
     n_leadtimes,
     n_cascade_levels=6,
@@ -97,7 +95,7 @@
 )
 
 # Back-transform to rain rate
-R_f = transformation.dB_transform(R_f, threshold=-10.0, inverse=True)[0]
+R_f = R_f.pysteps.db_transform(inverse=True)
 
 # Plot the S-PROG forecast
 plot_precip_field(
@@ -127,7 +125,7 @@
 # The STEPS nowcast
 nowcast_method = nowcasts.get_method("steps")
 R_f = nowcast_method(
-    R[-3:, :, :],
+    precip[-3:, :, :],
     V,
     n_leadtimes,
     n_ens_members,
@@ -142,7 +140,7 @@
 )
 
 # Back-transform to rain rates
-R_f = transformation.dB_transform(R_f, threshold=-10.0, inverse=True)[0]
+R_f = R_f.pysteps.db_transform(inverse=True)
 
 
 # Plot the ensemble mean
diff --git a/examples/probability_forecast.py b/examples/probability_forecast.py
index 059115bc1..2b350b2f0 100644
--- a/examples/probability_forecast.py
+++ b/examples/probability_forecast.py
@@ -75,7 +75,7 @@
 fns = io.find_by_date(date, root, fmt, pattern, ext, timestep, num_prev_files=2)
 importer = io.get_method(importer_name, "importer")
 precip, __, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
-precip, metadata = utils.to_rainrate(precip, metadata)
+precip = precip.pysteps.to_rainrate()
 # precip[np.isnan(precip)] = 0
 
 # motion
@@ -124,8 +124,8 @@
 fns = io.find_by_date(
     date, root, fmt, pattern, ext, timestep, num_next_files=nleadtimes
 )
-obs, __, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
-obs, metadata = utils.to_rainrate(obs, metadata)
+obs = io.read_timeseries(fns, importer, **importer_kwargs)
+obs = obs.pysteps.to_rainrate()
 obs[np.isnan(obs)] = 0
 
 # reliability diagram
diff --git a/examples/rainfarm_downscale.py b/examples/rainfarm_downscale.py
index be4033989..74ebb588d 100644
--- a/examples/rainfarm_downscale.py
+++ b/examples/rainfarm_downscale.py
@@ -18,7 +18,6 @@
 from pprint import pprint
 
 from pysteps import io, rcparams
-from pysteps.utils import aggregate_fields_space, square_domain, to_rainrate
 from pysteps.downscaling import rainfarm
 from pysteps.visualization import plot_precip_field
 
@@ -37,10 +36,7 @@
 )
 
 # Convert to mm/h
-precip, metadata = to_rainrate(precip, metadata)
-
-# Reduce to a square domain
-precip, metadata = square_domain(precip, metadata, "crop")
+precip = precip.pysteps.to_rainrate()
 
 # Nicely print the metadata
 pprint(metadata)
@@ -60,8 +56,7 @@
 # a lower resolution. We are going to use an upscaling factor of 16 x.
 
 upscaling_factor = 16
-upscale_to = metadata["xpixelsize"] * upscaling_factor  # upscaling factor : 16 x
-precip_lr, metadata_lr = aggregate_fields_space(precip, metadata, upscale_to)
+precip_lr = precip.coarsen(x=upscaling_factor, y=upscaling_factor).mean()
 
 # Plot the upscaled rainfall field
 plt.figure()
diff --git a/examples/thunderstorm_detection_and_tracking.py b/examples/thunderstorm_detection_and_tracking.py
index de3cc39e9..628379d15 100644
--- a/examples/thunderstorm_detection_and_tracking.py
+++ b/examples/thunderstorm_detection_and_tracking.py
@@ -34,7 +34,6 @@
 from pysteps import io, rcparams
 from pysteps.feature import tstorm as tstorm_detect
 from pysteps.tracking import tdating as tstorm_dating
-from pysteps.utils import to_reflectivity
 from pysteps.visualization import plot_precip_field, plot_track, plot_cart_contour
 
 ################################################################################
@@ -43,7 +42,7 @@
 #
 # A series of 20 files containing Swiss Cartesian gridded rain rates are imported. Since
 # the algorithm is tuned to Swiss max-reflectivity data, the rain rates are transformed
-# to reflectivity fields using the 'to_reflectivity' utility in pysteps.utils.
+# to reflectivity fields.
 
 # Select the input data
 date = datetime.strptime("201607112100", "%Y%m%d%H%M")
@@ -63,11 +62,11 @@
     date, root_path, path_fmt, fn_pattern, fn_ext, timestep, num_next_files=20
 )
 importer = io.get_method(importer_name, "importer")
-R, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
+precip, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
 
 # Convert to reflectivity (it is possible to give the a- and b- parameters of the
 # Marshall-Palmer relationship here: zr_a = and zr_b =).
-Z, metadata = to_reflectivity(R, metadata)
+reflectivity = precip.pysteps.to_reflectivity()
 
 # Extract the list of timestamps
 timelist = metadata["timestamps"]
@@ -80,7 +79,7 @@
 # The function tstorm_detect.detection requires a 2-D input image, all further inputs are
 # optional.
 
-input_image = Z[2, :, :].copy()
+input_image = reflectivity[2, :, :].copy()
 time = timelist[2]
 cells_id, labels = tstorm_detect.detection(input_image, time=time)
 
@@ -96,7 +95,7 @@
 # flow prediction and are not used to compute tracks.
 
 track_list, cell_list, label_list = tstorm_dating.dating(
-    input_video=Z, timelist=timelist
+    input_video=reflectivity, timelist=timelist
 )
 
 ###############################################################################
@@ -104,7 +103,7 @@
 # ~~~~~~~~~~~~~~~~~~~~
 
 # Plot precipitation field
-plot_precip_field(Z[2, :, :], geodata=metadata, units=metadata["unit"])
+plot_precip_field(reflectivity[2, :, :], geodata=metadata, units=metadata["unit"])
 plt.xlabel("Swiss easting [m]")
 plt.ylabel("Swiss northing [m]")
 
diff --git a/pysteps/accessors.py b/pysteps/accessors.py
new file mode 100644
index 000000000..44d393a01
--- /dev/null
+++ b/pysteps/accessors.py
@@ -0,0 +1,32 @@
+# -*- coding: utf-8 -*-
+"""
+"""
+from functools import partial, wraps
+
+import xarray as xr
+
+
+def _make_accessor(func, accessor_class):
+    """
+    Turn functions for xarray objects into accessors.
+    """
+
+    @wraps(func)
+    def method(accessor, *args, **kwargs):
+        return func(accessor._data, *args, **kwargs)
+
+    setattr(accessor_class, func.__name__, method)
+    return func
+
+
+@xr.register_dataarray_accessor("pysteps")
+class DataArrayUtils:
+    """
+    Xarray accessor class for pysteps utilities on DataArrays.
+    """
+
+    def __init__(self, da: xr.DataArray):
+        self._data = da
+
+
+dataarray_utils = partial(_make_accessor, accessor_class=DataArrayUtils)
\ No newline at end of file
diff --git a/pysteps/datasets.py b/pysteps/datasets.py
index 26b70d6db..f8b53b4d2 100644
--- a/pysteps/datasets.py
+++ b/pysteps/datasets.py
@@ -483,6 +483,6 @@ def load_dataset(case="fmi", frames=14):
     )
 
     # Convert to rain rate
-    precip, metadata = conversion.to_rainrate(reflectivity, metadata)
+    precip = reflectivity.pysteps.to_rainrate()
 
     return precip, metadata, data_source["timestep"]
diff --git a/pysteps/decorators.py b/pysteps/decorators.py
index b9cb401fb..0c81b73d0 100644
--- a/pysteps/decorators.py
+++ b/pysteps/decorators.py
@@ -40,7 +40,7 @@ def _add_extra_kwrds_to_docstrings(target_func, extra_kwargs_doc_text):
     return target_func
 
 
-def postprocess_import(fillna=np.nan, dtype="double"):
+def postprocess_import(fillna=np.nan, dtype="float32"):
     """
     Postprocess the imported precipitation data.
     Operations:
@@ -55,7 +55,7 @@ def postprocess_import(fillna=np.nan, dtype="double"):
     Parameters
     ----------
     dtype: str
-        Default data type for precipitation. Double precision by default.
+        Default data type for precipitation. Float32 precision by default.
     fillna: float or np.nan
         Default value used to represent the missing data ("No Coverage").
         By default, np.nan is used.
@@ -67,7 +67,6 @@ def postprocess_import(fillna=np.nan, dtype="double"):
     def _postprocess_import(importer):
         @wraps(importer)
         def _import_with_postprocessing(*args, **kwargs):
-
             data_array = importer(*args, **kwargs)
 
             if not isinstance(data_array, xr.DataArray):
@@ -161,83 +160,100 @@ def prepare_interpolator(nchunks=4):
 
     def _preamble_interpolation(interpolator):
         @wraps(interpolator)
-        def _interpolator_with_preamble(xy_coord, values, xgrid, ygrid, **kwargs):
+        def _interpolator_with_preamble(sparse_data, xgrid, ygrid, **kwargs):
             nonlocal nchunks  # https://stackoverflow.com/questions/5630409/
 
-            values = values.copy()
-            xy_coord = xy_coord.copy()
+            if not isinstance(sparse_data, (xr.Dataset, xr.DataArray)):
+                raise ValueError(
+                    "sparse_data must be an instance of xarray's Dataset or DataArray"
+                )
 
-            input_ndims = values.ndim
-            input_nvars = 1 if input_ndims == 1 else values.shape[1]
-            input_nsamples = values.shape[0]
+            if "sample" not in sparse_data.dims:
+                raise ValueError("missing dimension 'sample'")
 
-            coord_ndims = xy_coord.ndim
-            coord_nsamples = xy_coord.shape[0]
+            if not all([d in ("sample", "variable") for d in sparse_data.dims]):
+                raise ValueError(
+                    "sparse_data must have either dimension ('sample') or "
+                    f"('sample', 'variable'), but it has {sparse_data.dims}"
+                )
 
-            grid_shape = (ygrid.size, xgrid.size)
+            if not ("x" in sparse_data.coords and "y" in sparse_data.coords):
+                raise ValueError("missing coordinates 'x' and 'y'")
 
-            if np.any(~np.isfinite(values)):
-                raise ValueError("argument 'values' contains non-finite values")
-            if np.any(~np.isfinite(xy_coord)):
-                raise ValueError("argument 'xy_coord' contains non-finite values")
+            if np.isnan(xgrid).any() or np.isnan(ygrid).any():
+                raise ValueError("Target grid coordinates contain missing values")
 
-            if input_ndims > 2:
-                raise ValueError(
-                    "argument 'values' must have 1 (n) or 2 dimensions (n, m), "
-                    f"but it has {input_ndims}"
-                )
-            if not coord_ndims == 2:
-                raise ValueError(
-                    "argument 'xy_coord' must have 2 dimensions (n, 2), "
-                    f"but it has {coord_ndims}"
-                )
+            sparse_data = sparse_data.copy()
 
-            if not input_nsamples == coord_nsamples:
-                raise ValueError(
-                    "the number of samples in argument 'values' does not match the "
-                    f"number of coordinates {input_nsamples}!={coord_nsamples}"
-                )
+            if isinstance(sparse_data, xr.Dataset):
+                sparse_data = sparse_data.to_array(name="convert_to_dataset")
+            elif "variable" not in sparse_data.dims:
+                sparse_data = sparse_data.expand_dims("variable")
+
+            sparse_data = sparse_data.transpose("sample", "variable")
+
+            # drop missing values
+            sparse_data = sparse_data.dropna(dim="sample")
+            sparse_data = sparse_data.drop_isel(sample=np.isnan(sparse_data.x))
+            sparse_data = sparse_data.drop_isel(sample=np.isnan(sparse_data.y))
+
+            input_nvars = sparse_data.sizes["variable"]
+            input_nsamples = sparse_data.sizes["sample"]
+            input_dtype = sparse_data.dtype
+
+            grid_shape = (ygrid.size, xgrid.size)
+            output_grid = xr.DataArray(
+                np.zeros((input_nvars,) + grid_shape, dtype=input_dtype),
+                dims=("variable", "y", "x"),
+                coords={"y": ("y", ygrid), "x": ("x", xgrid)},
+            )
+            output_grid = output_grid.assign_coords(
+                sparse_data.drop_vars(("x", "y", "sample"), errors="ignore").coords
+            )
+            output_grid = output_grid.assign_attrs(sparse_data.attrs)
 
             # only one sample, return uniform output
             if input_nsamples == 1:
-                output_array = np.ones((input_nvars,) + grid_shape)
-                for n, v in enumerate(values[0, ...]):
-                    output_array[n, ...] *= v
-                return output_array.squeeze()
+                for n, v in enumerate(sparse_data.isel(sample=0)):
+                    output_grid[dict(variable=n)] += v
 
             # all equal elements, return uniform output
-            if values.max() == values.min():
-                return np.ones((input_nvars,) + grid_shape) * values.ravel()[0]
-
-            # split grid in n chunks
-            nchunks = int(kwargs.get("nchunks", nchunks) ** 0.5)
-            if nchunks > 1:
-                subxgrids = np.array_split(xgrid, nchunks)
-                subxgrids = [x for x in subxgrids if x.size > 0]
-                subygrids = np.array_split(ygrid, nchunks)
-                subygrids = [y for y in subygrids if y.size > 0]
-
-                # generate a unique identifier to be used for caching
-                # intermediate results
-                kwargs["hkey"] = uuid.uuid1().int
+            elif sparse_data.max() == sparse_data.min():
+                output_grid += sparse_data.values.ravel()[0]
+
+            # actual interpolation
             else:
-                subxgrids = [xgrid]
-                subygrids = [ygrid]
-
-            interpolated = np.zeros((input_nvars,) + grid_shape)
-            indx = 0
-            for subxgrid in subxgrids:
-                deltax = subxgrid.size
-                indy = 0
-                for subygrid in subygrids:
-                    deltay = subygrid.size
-                    interpolated[
-                        :, indy : (indy + deltay), indx : (indx + deltax)
-                    ] = interpolator(xy_coord, values, subxgrid, subygrid, **kwargs)
-                    indy += deltay
-                indx += deltax
-
-            return interpolated.squeeze()
+                # split grid in n chunks
+                nchunks = int(kwargs.get("nchunks", nchunks) ** 0.5)
+                if nchunks > 1:
+                    subxgrids = np.array_split(xgrid, nchunks)
+                    subxgrids = [x for x in subxgrids if x.size > 0]
+                    subygrids = np.array_split(ygrid, nchunks)
+                    subygrids = [y for y in subygrids if y.size > 0]
+
+                    # generate a unique identifier to be used for caching
+                    # intermediate results
+                    kwargs["hkey"] = uuid.uuid1().int
+                else:
+                    subxgrids = [xgrid]
+                    subygrids = [ygrid]
+
+                indx = 0
+                for subxgrid in subxgrids:
+                    deltax = subxgrid.size
+                    indy = 0
+                    for subygrid in subygrids:
+                        deltay = subygrid.size
+                        output_grid[
+                            :, indy : (indy + deltay), indx : (indx + deltax)
+                        ] = interpolator(sparse_data, subxgrid, subygrid, **kwargs)
+                        indy += deltay
+                    indx += deltax
+
+            if "convert_to_dataset" == sparse_data.name:
+                output_grid = output_grid.to_dataset(dim="variable")
+
+            return output_grid.squeeze()
 
         extra_kwargs_doc = """
             nchunks: int, optional
@@ -308,6 +324,7 @@ def _to_xarray(array, metadata):
 
     data_array.attrs.update(
         {
+            # TODO: Revise this list before final 2.0 version ?
             "unit": metadata["unit"],
             "accutime": metadata["accutime"],
             "transform": metadata["transform"],
@@ -315,13 +332,13 @@ def _to_xarray(array, metadata):
             "threshold": metadata["threshold"],
             "zr_a": metadata.get("zr_a", None),
             "zr_b": metadata.get("zr_b", None),
-            "institution": metadata["institution"],
+            "institution": metadata.get("institution", None),
             "projection": metadata["projection"],
-            #
-            # TODO: Remove before final 2.0 version
+            "bounding_box": (x1, x2, y1, y2),
             "yorigin": metadata["yorigin"],
             "xpixelsize": metadata["xpixelsize"],
             "ypixelsize": metadata["ypixelsize"],
+            "cartesian_unit": metadata["cartesian_unit"],
         }
     )
 
@@ -329,10 +346,6 @@ def _to_xarray(array, metadata):
         {
             "standard_name": "projection_x_coordinate",
             "units": metadata["cartesian_unit"],
-            #
-            # TODO: Remove before final 2.0 version
-            "x1": x1,
-            "x2": x2,
         }
     )
 
@@ -340,10 +353,6 @@ def _to_xarray(array, metadata):
         {
             "standard_name": "projection_y_coordinate",
             "units": metadata["cartesian_unit"],
-            #
-            # TODO: Remove before final 2.0 version
-            "y1": y1,
-            "y2": y2,
         }
     )
 
diff --git a/pysteps/feature/blob.py b/pysteps/feature/blob.py
index a51717416..a65212642 100644
--- a/pysteps/feature/blob.py
+++ b/pysteps/feature/blob.py
@@ -93,6 +93,8 @@ def detection(
     else:
         detector = feature.blob_doh
 
+    input_image = np.copy(input_image)
+
     blobs = detector(
         input_image,
         min_sigma=min_sigma,
diff --git a/pysteps/feature/shitomasi.py b/pysteps/feature/shitomasi.py
index 7de1199f0..f1e5d8ea1 100644
--- a/pysteps/feature/shitomasi.py
+++ b/pysteps/feature/shitomasi.py
@@ -119,7 +119,7 @@ def detection(
             "routine but it is not installed"
         )
 
-    input_image = input_image.copy()
+    input_image = np.copy(input_image)
 
     if input_image.ndim != 2:
         raise ValueError("input_image must be a two-dimensional array")
diff --git a/pysteps/feature/tstorm.py b/pysteps/feature/tstorm.py
index 49f830e7a..76d114651 100644
--- a/pysteps/feature/tstorm.py
+++ b/pysteps/feature/tstorm.py
@@ -115,6 +115,7 @@ def detection(
         raise MissingOptionalDependency(
             "pandas is required for thunderstorm DATing " "but it is not installed"
         )
+    input_image = np.copy(input_image)
     filt_image = np.zeros(input_image.shape)
     filt_image[input_image >= minref] = input_image[input_image >= minref]
     filt_image[input_image > maxref] = maxref
diff --git a/pysteps/io/exporters.py b/pysteps/io/exporters.py
index 9de6d2075..78ad1901f 100644
--- a/pysteps/io/exporters.py
+++ b/pysteps/io/exporters.py
@@ -275,7 +275,7 @@ def initialize_forecast_exporter_kineros(
         grids.
 
     metadata: dict
-        Metadata dictionary containing the projection,x1,x2,y1,y2 and unit
+        Metadata dictionary containing the projection,bounding_box and unit
         attributes described in the documentation of
         :py:mod:`pysteps.io.importers`.
 
@@ -424,6 +424,9 @@ def initialize_forecast_exporter_netcdf(
             "pyproj package is required for netcdf " "exporters but it is not installed"
         )
 
+    if not isinstance(startdate, datetime):
+        raise ValueError("argument 'stardate' must be a python datetime object")
+
     if incremental not in [None, "timestep", "member"]:
         raise ValueError(
             f"unknown option {incremental}: incremental must be "
@@ -483,9 +486,9 @@ def initialize_forecast_exporter_netcdf(
     else:
         raise ValueError("unknown unit %s" % metadata["unit"])
 
-    xr = np.linspace(metadata["x1"], metadata["x2"], w + 1)[:-1]
+    xr = np.linspace(*metadata["bounding_box"][:2], w + 1)[:-1]
     xr += 0.5 * (xr[1] - xr[0])
-    yr = np.linspace(metadata["y1"], metadata["y2"], h + 1)[:-1]
+    yr = np.linspace(*metadata["bounding_box"][2:], h + 1)[:-1]
     yr += 0.5 * (yr[1] - yr[0])
 
     # flip yr vector if yorigin is upper
@@ -497,14 +500,14 @@ def initialize_forecast_exporter_netcdf(
     var_xc.axis = "X"
     var_xc.standard_name = "projection_x_coordinate"
     var_xc.long_name = "x-coordinate in Cartesian system"
-    var_xc.units = metadata["cartesian_unit"]
+    var_xc.units = metadata.get("cartesian_unit", "m")
 
     var_yc = ncf.createVariable("y", np.float32, dimensions=("y",))
     var_yc[:] = yr
     var_yc.axis = "Y"
     var_yc.standard_name = "projection_y_coordinate"
     var_yc.long_name = "y-coordinate in Cartesian system"
-    var_yc.units = metadata["cartesian_unit"]
+    var_yc.units = metadata.get("cartesian_unit", "m")
 
     x_2d, y_2d = np.meshgrid(xr, yr)
     pr = pyproj.Proj(metadata["projection"])
@@ -550,8 +553,7 @@ def initialize_forecast_exporter_netcdf(
     if incremental != "timestep":
         var_time[:] = [i * timestep * 60 for i in range(1, n_timesteps + 1)]
     var_time.long_name = "forecast time"
-    startdate_str = datetime.strftime(startdate, "%Y-%m-%d %H:%M:%S")
-    var_time.units = "seconds since %s" % startdate_str
+    var_time.units = f"seconds since {startdate: %Y-%m-%d %H:%M:%S}"
 
     if incremental == "member" or n_ens_gt_one:
         var_f = ncf.createVariable(
@@ -868,10 +870,10 @@ def _create_geotiff_file(outfn, driver, shape, metadata, num_bands):
         gdal.GDT_Float32,
         ["COMPRESS=DEFLATE", "PREDICTOR=3"],
     )
-
-    sx = (metadata["x2"] - metadata["x1"]) / shape[1]
-    sy = (metadata["y2"] - metadata["y1"]) / shape[0]
-    dst.SetGeoTransform([metadata["x1"], sx, 0.0, metadata["y2"], 0.0, -sy])
+    x1, x2, y1, y2 = metadata["bounding_box"]
+    sx = (x2 - x1) / shape[1]
+    sy = (y2 - y1) / shape[0]
+    dst.SetGeoTransform([x1, sx, 0.0, y2, 0.0, -sy])
 
     sr = osr.SpatialReference()
     sr.ImportFromProj4(metadata["projection"])
diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py
index 8c19c1990..848fb308e 100644
--- a/pysteps/io/importers.py
+++ b/pysteps/io/importers.py
@@ -96,7 +96,6 @@
 from pysteps.decorators import postprocess_import
 from pysteps.exceptions import DataModelError
 from pysteps.exceptions import MissingOptionalDependency
-from pysteps.utils import aggregate_fields
 
 try:
     import gdalconst
@@ -234,7 +233,7 @@ def _get_threshold_value(precip):
 
 
 @postprocess_import(dtype="float32")
-def import_mrms_grib(filename, extent=None, window_size=4, **kwargs):
+def import_mrms_grib(filename, extent=None, **kwargs):
     """
     Importer for NSSL's Multi-Radar/Multi-Sensor System
     ([MRMS](https://www.nssl.noaa.gov/projects/mrms/)) rainrate product
@@ -251,12 +250,6 @@ def import_mrms_grib(filename, extent=None, window_size=4, **kwargs):
     precision, doubling the memory footprint.
     To change the precision of the data, use the *dtype* keyword.
 
-    Also, by default, the original data is downscaled by 4
-    (resulting in a ~4 km grid spacing).
-    In case that the original grid spacing is needed, use `window_size=1`.
-    But be aware that a single composite in double precipitation will
-    require 186 Mb of memory.
-
     Finally, if desired, the precipitation data can be extracted over a
     sub region of the full domain using the `extent` keyword.
     By default, the entire domain is returned.
@@ -282,11 +275,6 @@ def import_mrms_grib(filename, extent=None, window_size=4, **kwargs):
         By default (None), the entire domain is retrieved.
         The extent can be in any form that can be converted to a flat array
         of 4 elements array (e.g., lists or tuples).
-    window_size: array_like or int
-        Array containing down-sampling integer factor along each axis.
-        If an integer value is given, the same block shape is used for all the
-        image dimensions.
-        Default: window_size=4.
 
     {extra_kwargs_doc}
 
@@ -314,9 +302,6 @@ def import_mrms_grib(filename, extent=None, window_size=4, **kwargs):
     except OSError:
         raise OSError(f"Error opening NCEP's MRMS file. " f"File Not Found: {filename}")
 
-    if isinstance(window_size, int):
-        window_size = (window_size, window_size)
-
     if extent is not None:
         extent = np.asarray(extent)
         if (extent.ndim != 1) or (extent.size != 4):
@@ -347,27 +332,6 @@ def import_mrms_grib(filename, extent=None, window_size=4, **kwargs):
     precip = grib_msg.values
     no_data_mask = precip == -3  # Missing values
 
-    # Create a function with default arguments for aggregate_fields
-    block_reduce = partial(aggregate_fields, method="mean", trim=True)
-
-    if window_size != (1, 1):
-        # Downscale data
-        lats = block_reduce(lats, window_size[0])
-        lons = block_reduce(lons, window_size[1])
-
-        # Update the limits
-        ul_lat, lr_lat = lats[0], lats[-1]  # Lat from North to south!
-        ul_lon, lr_lon = lons[0], lons[-1]
-
-        precip[no_data_mask] = 0  # block_reduce does not handle nan values
-        precip = block_reduce(precip, window_size, axis=(0, 1))
-
-        # Consider that if a single invalid observation is located in the block,
-        # then mark that value as invalid.
-        no_data_mask = block_reduce(
-            no_data_mask.astype("int"), window_size, axis=(0, 1)
-        ).astype(bool)
-
     lons, lats = np.meshgrid(lons, lats)
     precip[no_data_mask] = np.nan
 
@@ -393,8 +357,8 @@ def import_mrms_grib(filename, extent=None, window_size=4, **kwargs):
     pr = pyproj.Proj(proj_params)
     proj_def = " ".join([f"+{key}={value} " for key, value in proj_params.items()])
 
-    xsize = grib_msg["iDirectionIncrementInDegrees"] * window_size[0]
-    ysize = grib_msg["jDirectionIncrementInDegrees"] * window_size[1]
+    xsize = grib_msg["iDirectionIncrementInDegrees"]
+    ysize = grib_msg["jDirectionIncrementInDegrees"]
 
     x1, y1 = pr(ul_lon, lr_lat)
     x2, y2 = pr(lr_lon, ul_lat)
diff --git a/pysteps/io/nowcast_importers.py b/pysteps/io/nowcast_importers.py
index 5bc7db0c3..3aadf5d98 100755
--- a/pysteps/io/nowcast_importers.py
+++ b/pysteps/io/nowcast_importers.py
@@ -106,6 +106,8 @@ def import_netcdf_pysteps(filename, onerror="warn", **kwargs):
         Precipitation field in mm/h. The dimensions are [latitude, longitude].
         The first grid point (0,0) corresponds to the upper left corner of the
         domain, while (last i, last j) denote the lower right corner.
+    quality: None
+        Not implement.
     metadata: dict
         Associated metadata (pixel sizes, map projections, etc.).
     """
@@ -205,9 +207,13 @@ def import_netcdf_pysteps(filename, onerror="warn", **kwargs):
         metadata["zerovalue"] = np.nanmin(precip)
         metadata["threshold"] = np.nanmin(precip[precip > np.nanmin(precip)])
 
+        metadata["institution"] = ds.getncattr("institution")
+
+        metadata["cartesian_unit"] = ds["x"].getncattr("units")
+
         ds.close()
 
-        return precip, metadata
+        return precip, None, metadata
     except Exception as er:
         if onerror == "warn":
             print("There was an error processing the file", er)
diff --git a/pysteps/scripts/run_vel_pert_analysis.py b/pysteps/scripts/run_vel_pert_analysis.py
index 8ae241553..aa72ad42a 100644
--- a/pysteps/scripts/run_vel_pert_analysis.py
+++ b/pysteps/scripts/run_vel_pert_analysis.py
@@ -96,7 +96,7 @@
     R[~np.isfinite(R)] = metadata["zerovalue"]
     if use_precip_mask:
         MASK = np.any(R < R_min, axis=0)
-    R = transformation.dB_transform(R)[0]
+    R = R.pysteps.db_transform()
 
     if args.oflow == "vet":
         R_ = R[-2:, :, :]
diff --git a/pysteps/utils/__init__.py b/pysteps/utils/__init__.py
index 9594a75ae..c418f7186 100644
--- a/pysteps/utils/__init__.py
+++ b/pysteps/utils/__init__.py
@@ -3,7 +3,6 @@
 from .arrays import *
 from .cleansing import *
 from .conversion import *
-from .dimension import *
 from .images import *
 from .interface import get_method
 from .interpolate import *
diff --git a/pysteps/utils/arrays.py b/pysteps/utils/arrays.py
index 4644ad874..b30632682 100644
--- a/pysteps/utils/arrays.py
+++ b/pysteps/utils/arrays.py
@@ -13,20 +13,20 @@
 import numpy as np
 
 
-def compute_centred_coord_array(M, N):
+def compute_centred_coord_array(height, width):
     """Compute a 2D coordinate array, where the origin is at the center.
 
     Parameters
     ----------
-    M : int
+    height : int
       The height of the array.
-    N : int
+    width : int
       The width of the array.
 
     Returns
     -------
-    out : ndarray
-      The coordinate array.
+    out : tuple of ndarray
+      The two coordinate arrays.
 
     Examples
     --------
@@ -40,16 +40,16 @@ def compute_centred_coord_array(M, N):
 
     """
 
-    if M % 2 == 1:
-        s1 = np.s_[-int(M / 2) : int(M / 2) + 1]
+    if height % 2 == 1:
+        s1 = np.s_[-int(height / 2) : int(height / 2) + 1]
     else:
-        s1 = np.s_[-int(M / 2) : int(M / 2)]
+        s1 = np.s_[-int(height / 2) : int(height / 2)]
 
-    if N % 2 == 1:
-        s2 = np.s_[-int(N / 2) : int(N / 2) + 1]
+    if width % 2 == 1:
+        s2 = np.s_[-int(width / 2) : int(width / 2) + 1]
     else:
-        s2 = np.s_[-int(N / 2) : int(N / 2)]
+        s2 = np.s_[-int(width / 2) : int(width / 2)]
 
-    YC, XC = np.ogrid[s1, s2]
+    y_coordinates, x_coordinates = np.ogrid[s1, s2]
 
-    return YC, XC
+    return y_coordinates, x_coordinates
diff --git a/pysteps/utils/conversion.py b/pysteps/utils/conversion.py
index 08411344a..668d4510f 100644
--- a/pysteps/utils/conversion.py
+++ b/pysteps/utils/conversion.py
@@ -12,308 +12,153 @@
     to_raindepth
     to_reflectivity
 """
+from typing import Any
 
-import warnings
-from . import transformation
+import xarray as xr
 
-# TODO: This should not be done. Instead fix the code so that it doesn't
-# produce the warnings.
-# to deactivate warnings for comparison operators with NaNs
-warnings.filterwarnings("ignore", category=RuntimeWarning)
+from pysteps.accessors import dataarray_utils
 
 
-def to_rainrate(R, metadata, zr_a=None, zr_b=None):
-    """Convert to rain rate [mm/h].
+@dataarray_utils
+def to_rainrate(data_array: xr.DataArray, **kwargs: Any) -> xr.DataArray:
+    """Convert to rain rates [mm/h].
 
     Parameters
     ----------
-    R: array-like
-        Array of any shape to be (back-)transformed.
-    metadata: dict
-        Metadata dictionary containing the accutime, transform, unit, threshold
-        and zerovalue attributes as described in the documentation of
-        :py:mod:`pysteps.io.importers`.
-
-        Additionally, in case of conversion to/from reflectivity units, the
-        zr_a and zr_b attributes are also required,
-        but only if zr_a = zr_b = None.
-        If missing, it defaults to Marshall–Palmer relation,
-        that is, zr_a = 200.0 and zr_b = 1.6.
-    zr_a, zr_b: float, optional
-        The a and b coefficients of the Z-R relationship (Z = a*R^b).
+    data_array : xr.DataArray
+        DataArray to be converted.
+
+    Other parameters
+    ----------------
+    accutime: float
+        The accumulation time in minutes, used if converting from rain depth.
+        If not passed as an argument, the accutime specified in the
+        attributes of the input xr.DataArray is used.
+
+    zr_a, zr_b: float
+        The a and b coefficients of the Z-R relationship (Z = a*R^b), used if
+        converting from Z units.
+        If not passed as arguments nor included as attributes in the input
+        xr.DataArray, the default Marshall-Palmer is used (a=200, b=1.6).
 
     Returns
     -------
-    R: array-like
-        Array of any shape containing the converted units.
-    metadata: dict
-        The metadata with updated attributes.
+    data_array : xr.DataArray
+        DataArray with units converted to mm/h.
     """
+    data_array = data_array.pysteps.back_transform()
+    attrs = data_array.attrs
 
-    R = R.copy()
-    metadata = metadata.copy()
-
-    if metadata["transform"] is not None:
-
-        if metadata["transform"] == "dB":
-
-            R, metadata = transformation.dB_transform(R, metadata, inverse=True)
-
-        elif metadata["transform"] in ["BoxCox", "log"]:
-
-            R, metadata = transformation.boxcox_transform(R, metadata, inverse=True)
-
-        elif metadata["transform"] == "NQT":
-
-            R, metadata = transformation.NQ_transform(R, metadata, inverse=True)
-
-        elif metadata["transform"] == "sqrt":
-
-            R, metadata = transformation.sqrt_transform(R, metadata, inverse=True)
-
-        else:
-
-            raise ValueError("Unknown transformation %s" % metadata["transform"])
-
-    if metadata["unit"] == "mm/h":
+    units = attrs["unit"]
+    accutime = kwargs.get("accutime", data_array.attrs["accutime"])
+    zr_a = kwargs.get("zr_a", data_array.attrs.get("zr_a", 200.0))
+    zr_b = kwargs.get("zr_b", data_array.attrs.get("zr_b", 1.6))
 
+    if units == "mm/h":
         pass
 
-    elif metadata["unit"] == "mm":
-
-        threshold = metadata["threshold"]  # convert the threshold, too
-        zerovalue = metadata["zerovalue"]  # convert the zerovalue, too
-
-        R = R / float(metadata["accutime"]) * 60.0
-        threshold = threshold / float(metadata["accutime"]) * 60.0
-        zerovalue = zerovalue / float(metadata["accutime"]) * 60.0
-
-        metadata["threshold"] = threshold
-        metadata["zerovalue"] = zerovalue
+    elif units == "mm":
+        data_array = data_array / float(accutime) * 60.0
 
-    elif metadata["unit"] == "dBZ":
-
-        threshold = metadata["threshold"]  # convert the threshold, too
-        zerovalue = metadata["zerovalue"]  # convert the zerovalue, too
-
-        # Z to R
-        if zr_a is None:
-            zr_a = metadata.get("zr_a", 200.0)  # default to Marshall–Palmer
-        if zr_b is None:
-            zr_b = metadata.get("zr_b", 1.6)  # default to Marshall–Palmer
-        R = (R / zr_a) ** (1.0 / zr_b)
-        threshold = (threshold / zr_a) ** (1.0 / zr_b)
-        zerovalue = (zerovalue / zr_a) ** (1.0 / zr_b)
-
-        metadata["zr_a"] = zr_a
-        metadata["zr_b"] = zr_b
-        metadata["threshold"] = threshold
-        metadata["zerovalue"] = zerovalue
+    elif units == "Z" or units == "dBZ":
+        data_array = (data_array / zr_a) ** (1.0 / zr_b)
 
     else:
-        raise ValueError(
-            "Cannot convert unit %s and transform %s to mm/h"
-            % (metadata["unit"], metadata["transform"])
-        )
+        raise ValueError(f"Unsupported unit conversion from {units}.")
 
-    metadata["unit"] = "mm/h"
+    attrs.update({"unit": "mm/h"})
+    data_array.attrs = attrs
 
-    return R, metadata
+    return data_array
 
 
-def to_raindepth(R, metadata, zr_a=None, zr_b=None):
+@dataarray_utils
+def to_raindepth(data_array: xr.DataArray, **kwargs: Any) -> xr.DataArray:
     """Convert to rain depth [mm].
 
     Parameters
     ----------
-    R: array-like
-        Array of any shape to be (back-)transformed.
-    metadata: dict
-        Metadata dictionary containing the accutime, transform, unit, threshold
-        and zerovalue attributes as described in the documentation of
-        :py:mod:`pysteps.io.importers`.
-
-        Additionally, in case of conversion to/from reflectivity units, the
-        zr_a and zr_b attributes are also required,
-        but only if zr_a = zr_b = None.
-        If missing, it defaults to Marshall–Palmer relation, that is,
-        zr_a = 200.0 and zr_b = 1.6.
-    zr_a, zr_b: float, optional
-        The a and b coefficients of the Z-R relationship (Z = a*R^b).
+    data_array : xr.DataArray
+        DataArray to be converted.
+
+    Other parameters
+    ----------------
+    accutime: float
+        The accumulation time in minutes.
+        If not passed as an argument, the accutime specified in the
+        attributes of the input xr.DataArray is used.
+
+    zr_a, zr_b: float
+        The a and b coefficients of the Z-R relationship (Z = a*R^b), used if
+        converting from dBZ.
+        If not passed as arguments nor include as attributes in the input
+        xr.DataArray, the default Marshall-Palmer is used (a=200, b=1.6).
 
     Returns
     -------
-    R: array-like
-        Array of any shape containing the converted units.
-    metadata: dict
-        The metadata with updated attributes.
+    data_array : xr.DataArray
+        DataArray with units converted to mm.
     """
+    data_array = data_array.pysteps.to_rainrate(**kwargs)
+    attrs = data_array.attrs
 
-    R = R.copy()
-    metadata = metadata.copy()
-
-    if metadata["transform"] is not None:
-
-        if metadata["transform"] == "dB":
-
-            R, metadata = transformation.dB_transform(R, metadata, inverse=True)
-
-        elif metadata["transform"] in ["BoxCox", "log"]:
+    accutime = kwargs.get("accutime", attrs["accutime"])
+    data_array = data_array / 60.0 * accutime
 
-            R, metadata = transformation.boxcox_transform(R, metadata, inverse=True)
-
-        elif metadata["transform"] == "NQT":
-
-            R, metadata = transformation.NQ_transform(R, metadata, inverse=True)
-
-        elif metadata["transform"] == "sqrt":
-
-            R, metadata = transformation.sqrt_transform(R, metadata, inverse=True)
-
-        else:
-            raise ValueError("Unknown transformation %s" % metadata["transform"])
-
-    if metadata["unit"] == "mm" and metadata["transform"] is None:
-        pass
+    attrs.update({"unit": "mm"})
+    data_array.attrs = attrs
 
-    elif metadata["unit"] == "mm/h":
+    return data_array
 
-        threshold = metadata["threshold"]  # convert the threshold, too
-        zerovalue = metadata["zerovalue"]  # convert the zerovalue, too
 
-        R = R / 60.0 * metadata["accutime"]
-        threshold = threshold / 60.0 * metadata["accutime"]
-        zerovalue = zerovalue / 60.0 * metadata["accutime"]
+@dataarray_utils
+def to_reflectivity(
+    data_array: xr.DataArray, to_decibels: bool = True, offset: float = 0.01, **kwargs: Any
+) -> xr.DataArray:
+    """Convert to linear reflectivity [Z] or dBZ units.
 
-        metadata["threshold"] = threshold
-        metadata["zerovalue"] = zerovalue
-
-    elif metadata["unit"] == "dBZ":
-
-        threshold = metadata["threshold"]  # convert the threshold, too
-        zerovalue = metadata["zerovalue"]  # convert the zerovalue, too
-
-        # Z to R
-        if zr_a is None:
-            zr_a = metadata.get("zr_a", 200.0)  # Default to Marshall–Palmer
-        if zr_b is None:
-            zr_b = metadata.get("zr_b", 1.6)  # Default to Marshall–Palmer
-        R = (R / zr_a) ** (1.0 / zr_b) / 60.0 * metadata["accutime"]
-        threshold = (threshold / zr_a) ** (1.0 / zr_b) / 60.0 * metadata["accutime"]
-        zerovalue = (zerovalue / zr_a) ** (1.0 / zr_b) / 60.0 * metadata["accutime"]
-
-        metadata["zr_a"] = zr_a
-        metadata["zr_b"] = zr_b
-        metadata["threshold"] = threshold
-        metadata["zerovalue"] = zerovalue
-
-    else:
-        raise ValueError(
-            "Cannot convert unit %s and transform %s to mm"
-            % (metadata["unit"], metadata["transform"])
-        )
-
-    metadata["unit"] = "mm"
+    Parameters
+    ----------
+    data_array : xr.DataArray
+        DataArray to be converted.
 
-    return R, metadata
+    to_decibels: bool
+        If true, transform to decibels [dB].
 
+    offset: float
+        A small constant offset which prevents the division by zero or the
+        computation of log(0). Used if to_decibels is True.
 
-def to_reflectivity(R, metadata, zr_a=None, zr_b=None):
-    """Convert to reflectivity [dBZ].
+    Other parameters
+    ----------------
+    accutime: float
+        The accumulation time in minutes, used if converting from rain depth.
+        If not passed as an argument, the accutime specified in the
+        attributes of the input xr.DataArray is used.
 
-    Parameters
-    ----------
-    R: array-like
-        Array of any shape to be (back-)transformed.
-    metadata: dict
-        Metadata dictionary containing the accutime, transform, unit, threshold
-        and zerovalue attributes as described in the documentation of
-        :py:mod:`pysteps.io.importers`.
-
-        Additionally, in case of conversion to/from reflectivity units, the
-        zr_a and zr_b attributes are also required,
-        but only if zr_a = zr_b = None.
-        If missing, it defaults to Marshall–Palmer relation, that is,
-        zr_a = 200.0 and zr_b = 1.6.
-    zr_a, zr_b: float, optional
-        The a and b coefficients of the Z-R relationship (Z = a*R^b).
+    zr_a, zr_b: float
+        The a and b coefficients of the Z-R relationship `Z = a * R ** b`.
+        If not passed as arguments nor include as attributes in the input
+        xr.DataArray, the default Marshall-Palmer is used (a=200, b=1.6).
 
     Returns
     -------
-    R: array-like
-        Array of any shape containing the converted units.
-    metadata: dict
-        The metadata with updated attributes.
+    data_array : xr.DataArray
+        DataArray with units converted to Z or dBZ units.
     """
+    data_array = data_array.pysteps.to_rainrate(**kwargs)
+    attrs = data_array.attrs
 
-    R = R.copy()
-    metadata = metadata.copy()
-
-    if metadata["transform"] is not None:
-
-        if metadata["transform"] == "dB":
-
-            R, metadata = transformation.dB_transform(R, metadata, inverse=True)
+    zr_a = kwargs.get("zr_a", attrs.get("zr_a", 200.0))
+    zr_b = kwargs.get("zr_b", attrs.get("zr_b", 1.6))
+    data_array = zr_a * data_array ** zr_b
 
-        elif metadata["transform"] in ["BoxCox", "log"]:
+    attrs.update({"unit": "Z"})
 
-            R, metadata = transformation.boxcox_transform(R, metadata, inverse=True)
+    if to_decibels:
+        attrs.update({"unit": "dBZ"})
+        data_array = data_array.pysteps.db_transform(offset=offset)
 
-        elif metadata["transform"] == "NQT":
+    data_array.attrs = attrs
 
-            R, metadata = transformation.NQ_transform(R, metadata, inverse=True)
-
-        elif metadata["transform"] == "sqrt":
-
-            R, metadata = transformation.sqrt_transform(R, metadata, inverse=True)
-
-        else:
-
-            raise ValueError("Unknown transformation %s" % metadata["transform"])
-
-    if metadata["unit"] == "mm/h":
-
-        # Z to R
-        if zr_a is None:
-            zr_a = metadata.get("zr_a", 200.0)  # Default to Marshall–Palmer
-        if zr_b is None:
-            zr_b = metadata.get("zr_b", 1.6)  # Default to Marshall–Palmer
-
-        R = zr_a * R ** zr_b
-        metadata["threshold"] = zr_a * metadata["threshold"] ** zr_b
-        metadata["zerovalue"] = zr_a * metadata["zerovalue"] ** zr_b
-        metadata["zr_a"] = zr_a
-        metadata["zr_b"] = zr_b
-
-        # Z to dBZ
-        R, metadata = transformation.dB_transform(R, metadata)
-
-    elif metadata["unit"] == "mm":
-        # depth to rate
-        R, metadata = to_rainrate(R, metadata)
-
-        # Z to R
-        if zr_a is None:
-            zr_a = metadata.get("zr_a", 200.0)  # Default to Marshall-Palmer
-        if zr_b is None:
-            zr_b = metadata.get("zr_b", 1.6)  # Default to Marshall-Palmer
-        R = zr_a * R ** zr_b
-        metadata["threshold"] = zr_a * metadata["threshold"] ** zr_b
-        metadata["zerovalue"] = zr_a * metadata["zerovalue"] ** zr_b
-        metadata["zr_a"] = zr_a
-        metadata["zr_b"] = zr_b
-
-        # Z to dBZ
-        R, metadata = transformation.dB_transform(R, metadata)
-
-    elif metadata["unit"] == "dBZ":
-        # Z to dBZ
-        R, metadata = transformation.dB_transform(R, metadata)
-
-    else:
-        raise ValueError(
-            "Cannot convert unit %s and transform %s to mm/h"
-            % (metadata["unit"], metadata["transform"])
-        )
-    metadata["unit"] = "dBZ"
-    return R, metadata
+    return data_array
diff --git a/pysteps/utils/dimension.py b/pysteps/utils/dimension.py
deleted file mode 100644
index 5e1b1ff9d..000000000
--- a/pysteps/utils/dimension.py
+++ /dev/null
@@ -1,596 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-pysteps.utils.dimension
-=======================
-
-Functions to manipulate array dimensions.
-
-.. autosummary::
-    :toctree: ../generated/
-
-    aggregate_fields
-    aggregate_fields_time
-    aggregate_fields_space
-    clip_domain
-    square_domain
-"""
-
-import numpy as np
-
-_aggregation_methods = dict(
-    sum=np.sum, mean=np.mean, nanmean=np.nanmean, nansum=np.nansum
-)
-
-
-def aggregate_fields_time(R, metadata, time_window_min, ignore_nan=False):
-    """Aggregate fields in time.
-
-    Parameters
-    ----------
-    R: array-like
-        Array of shape (t,m,n) or (l,t,m,n) containing
-        a time series of (ensemble) input fields.
-        They must be evenly spaced in time.
-    metadata: dict
-        Metadata dictionary containing the timestamps and unit attributes as
-        described in the documentation of :py:mod:`pysteps.io.importers`.
-    time_window_min: float or None
-        The length in minutes of the time window that is used to
-        aggregate the fields.
-        The time spanned by the t dimension of R must be a multiple of
-        time_window_min.
-        If set to None, it returns a copy of the original R and metadata.
-    ignore_nan: bool, optional
-        If True, ignore nan values.
-
-    Returns
-    -------
-    outputarray: array-like
-        The new array of aggregated fields of shape (k,m,n) or (l,k,m,n), where
-        k = t*delta/time_window_min and delta is the time interval between two
-        successive timestamps.
-    metadata: dict
-        The metadata with updated attributes.
-
-    See also
-    --------
-    pysteps.utils.dimension.aggregate_fields_space,
-    pysteps.utils.dimension.aggregate_fields
-    """
-
-    R = R.copy()
-    metadata = metadata.copy()
-
-    if time_window_min is None:
-        return R, metadata
-
-    unit = metadata["unit"]
-    timestamps = metadata["timestamps"]
-    if "leadtimes" in metadata:
-        leadtimes = metadata["leadtimes"]
-
-    if len(R.shape) < 3:
-        raise ValueError("The number of dimension must be > 2")
-    if len(R.shape) == 3:
-        axis = 0
-    if len(R.shape) == 4:
-        axis = 1
-    if len(R.shape) > 4:
-        raise ValueError("The number of dimension must be <= 4")
-
-    if R.shape[axis] != len(timestamps):
-        raise ValueError(
-            "The list of timestamps has length %i, " % len(timestamps)
-            + "but R contains %i frames" % R.shape[axis]
-        )
-
-    # assumes that frames are evenly spaced
-    delta = (timestamps[1] - timestamps[0]).seconds / 60
-    if delta == time_window_min:
-        return R, metadata
-    if (R.shape[axis] * delta) % time_window_min:
-        raise ValueError("time_window_size does not equally split R")
-
-    nframes = int(time_window_min / delta)
-
-    # specify the operator to be used to aggregate
-    # the values within the time window
-    if unit == "mm/h":
-        method = "mean"
-    elif unit == "mm":
-        method = "sum"
-    else:
-        raise ValueError(
-            "can only aggregate units of 'mm/h' or 'mm'" + " not %s" % unit
-        )
-
-    if ignore_nan:
-        method = "".join(("nan", method))
-
-    R = aggregate_fields(R, nframes, axis=axis, method=method)
-
-    metadata["accutime"] = time_window_min
-    metadata["timestamps"] = timestamps[nframes - 1 :: nframes]
-    if "leadtimes" in metadata:
-        metadata["leadtimes"] = leadtimes[nframes - 1 :: nframes]
-
-    return R, metadata
-
-
-def aggregate_fields_space(R, metadata, space_window, ignore_nan=False):
-    """Upscale fields in space.
-
-    Parameters
-    ----------
-    R: array-like
-        Array of shape (m,n), (t,m,n) or (l,t,m,n) containing a single field or
-        a time series of (ensemble) input fields.
-    metadata: dict
-        Metadata dictionary containing the xpixelsize, ypixelsize and unit
-        attributes as described in the documentation of
-        :py:mod:`pysteps.io.importers`.
-    space_window: float or None
-        The length of the space window that is used to upscale the fields.
-        The space_window unit is the same used in the geographical projection
-        of R
-        and hence the same as for the xpixelsize and ypixelsize attributes.
-        The space spanned by the m and n dimensions of R must be a multiple of
-        space_window.
-        If set to None, it returns a copy of the original R and metadata.
-    ignore_nan: bool, optional
-        If True, ignore nan values.
-
-    Returns
-    -------
-    outputarray: array-like
-        The new array of aggregated fields of shape (k,j), (t,k,j)
-        or (l,t,k,j),
-        where k = m*ypixelsize/space_window and j = n*xpixelsize/space_window.
-    metadata: dict
-        The metadata with updated attributes.
-
-    See also
-    --------
-    pysteps.utils.dimension.aggregate_fields_time,
-    pysteps.utils.dimension.aggregate_fields
-    """
-
-    R = R.copy()
-    metadata = metadata.copy()
-
-    if space_window is None:
-        return R, metadata
-
-    unit = metadata["unit"]
-    ypixelsize = metadata["ypixelsize"]
-    xpixelsize = metadata["xpixelsize"]
-
-    if len(R.shape) < 2:
-        raise ValueError("The number of dimensions must be >= 2")
-    if len(R.shape) == 2:
-        axes = [0, 1]
-    if len(R.shape) == 3:
-        axes = [1, 2]
-    if len(R.shape) == 4:
-        axes = [2, 3]
-    if len(R.shape) > 4:
-        raise ValueError("The number of dimensions must be <= 4")
-
-    # assumes that frames are evenly spaced
-    if ypixelsize == space_window and xpixelsize == space_window:
-        return R, metadata
-    if (R.shape[axes[0]] * ypixelsize) % space_window or (
-        R.shape[axes[1]] * xpixelsize
-    ) % space_window:
-        raise ValueError("space_window does not equally split R")
-
-    nframes = [int(space_window / ypixelsize), int(space_window / xpixelsize)]
-
-    # specify the operator to be used to aggregate the values
-    # within the space window
-    if unit == "mm/h":
-        method = "mean"
-    elif unit == "mm":
-        method = "sum"
-    else:
-        raise ValueError(
-            "can only aggregate units of 'mm/h' or 'mm' " + "not %s" % unit
-        )
-
-    if ignore_nan:
-        method = "".join(("nan", method))
-
-    R = aggregate_fields(R, nframes[0], axis=axes[0], method=method)
-    R = aggregate_fields(R, nframes[1], axis=axes[1], method=method)
-
-    metadata["ypixelsize"] = space_window
-    metadata["xpixelsize"] = space_window
-
-    return R, metadata
-
-
-def aggregate_fields(data, window_size, axis=0, method="mean", trim=False):
-    """Aggregate fields along a given direction.
-
-    It attempts to aggregate the given R axis in an integer number of sections
-    of length = ``window_size``.
-    If such a aggregation is not possible, an error is raised unless ``trim``
-    set to True, in which case the axis is trimmed (from the end)
-    to make it perfectly divisible".
-
-    Parameters
-    ----------
-    data: array-like
-        Array of any shape containing the input fields.
-    window_size: int or tuple of ints
-        The length of the window that is used to aggregate the fields.
-        If a single integer value is given, the same window is used for
-        all the selected axis.
-
-        If ``window_size`` is a 1D array-like,
-        each element indicates the length of the window that is used
-        to aggregate the fields along each axis. In this case,
-        the number of elements of 'window_size' must be the same as the elements
-        in the ``axis`` argument.
-    axis: int or array-like of ints
-        Axis or axes where to perform the aggregation.
-        If this is a tuple of ints, the aggregation is performed over multiple
-        axes, instead of a single axis
-    method: string, optional
-        Optional argument that specifies the operation to use
-        to aggregate the values within the window.
-        Default to mean operator.
-    trim: bool
-         In case that the ``data`` is not perfectly divisible by
-         ``window_size`` along the selected axis:
-
-         - trim=True: the data will be trimmed (from the end) along that
-           axis to make it perfectly divisible.
-         - trim=False: a ValueError exception is raised.
-
-    Returns
-    -------
-    new_array: array-like
-        The new aggregated array with shape[axis] = k,
-        where k = R.shape[axis] / window_size.
-
-    See also
-    --------
-    pysteps.utils.dimension.aggregate_fields_time,
-    pysteps.utils.dimension.aggregate_fields_space
-    """
-
-    if np.ndim(axis) > 1:
-        raise TypeError(
-            "Only integers or integer 1D arrays can be used for the " "'axis' argument."
-        )
-
-    if np.ndim(axis) == 1:
-        axis = np.asarray(axis)
-        if np.ndim(window_size) == 0:
-            window_size = (window_size,) * axis.size
-
-        window_size = np.asarray(window_size, dtype="int")
-
-        if window_size.shape != axis.shape:
-            raise ValueError(
-                "The 'window_size' and 'axis' shapes are incompatible."
-                f"window_size.shape: {str(window_size.shape)}, "
-                f"axis.shape: {str(axis.shape)}, "
-            )
-
-        new_data = data.copy()
-        for i in range(axis.size):
-            # Recursively call the aggregate_fields function
-            new_data = aggregate_fields(
-                new_data, window_size[i], axis=axis[i], method=method, trim=trim
-            )
-
-        return new_data
-
-    if np.ndim(window_size) != 0:
-        raise TypeError(
-            "A single axis was selected for the aggregation but several"
-            f"of window_sizes were given: {str(window_size)}."
-        )
-
-    data = np.asarray(data).copy()
-    orig_shape = data.shape
-
-    if method not in _aggregation_methods:
-        raise ValueError(
-            "Aggregation method not recognized. "
-            f"Available methods: {str(list(_aggregation_methods.keys()))}"
-        )
-
-    if window_size <= 0:
-        raise ValueError("'window_size' must be strictly positive")
-
-    if (orig_shape[axis] % window_size) and (not trim):
-        raise ValueError(
-            f"Since 'trim' argument was set to False,"
-            f"the 'window_size' {window_size} must exactly divide"
-            f"the dimension along the selected axis:"
-            f"data.shape[axis]={orig_shape[axis]}"
-        )
-
-    new_data = data.swapaxes(axis, 0)
-    if trim:
-        trim_size = data.shape[axis] % window_size
-        if trim_size > 0:
-            new_data = new_data[:-trim_size]
-
-    new_data_shape = list(new_data.shape)
-    new_data_shape[0] //= window_size  # Final shape
-
-    new_data = new_data.reshape(new_data_shape[0], window_size, -1)
-
-    new_data = _aggregation_methods[method](new_data, axis=1)
-
-    new_data = new_data.reshape(new_data_shape).swapaxes(axis, 0)
-
-    return new_data
-
-
-def clip_domain(R, metadata, extent=None):
-    """Clip the field domain by geographical coordinates.
-
-    Parameters
-    ----------
-    R: array-like
-        Array of shape (m,n) or (t,m,n) containing the input fields.
-    metadata: dict
-        Metadata dictionary containing the x1, x2, y1, y2,
-        xpixelsize, ypixelsize,
-        zerovalue and yorigin attributes as described in the documentation of
-        :py:mod:`pysteps.io.importers`.
-    extent: scalars (left, right, bottom, top), optional
-        The extent of the bounding box in data coordinates to be used to clip
-        the data.
-        Note that the direction of the vertical axis and thus the default
-        values for top and bottom depend on origin. We follow the same
-        convention as in the imshow method of matplotlib:
-        https://matplotlib.org/tutorials/intermediate/imshow_extent.html
-
-    Returns
-    -------
-    R: array-like
-        the clipped array
-    metadata: dict
-        the metadata with updated attributes.
-    """
-
-    R = R.copy()
-    R_shape = np.array(R.shape)
-    metadata = metadata.copy()
-
-    if extent is None:
-        return R, metadata
-
-    if len(R.shape) < 2:
-        raise ValueError("The number of dimension must be > 1")
-    if len(R.shape) == 2:
-        R = R[None, None, :, :]
-    if len(R.shape) == 3:
-        R = R[None, :, :, :]
-    if len(R.shape) > 4:
-        raise ValueError("The number of dimension must be <= 4")
-
-    # extract original domain coordinates
-    left = metadata["x1"]
-    right = metadata["x2"]
-    bottom = metadata["y1"]
-    top = metadata["y2"]
-
-    # extract bounding box coordinates
-    left_ = extent[0]
-    right_ = extent[1]
-    bottom_ = extent[2]
-    top_ = extent[3]
-
-    # compute its extent in pixels
-    dim_x_ = int((right_ - left_) / metadata["xpixelsize"])
-    dim_y_ = int((top_ - bottom_) / metadata["ypixelsize"])
-    R_ = np.ones((R.shape[0], R.shape[1], dim_y_, dim_x_)) * metadata["zerovalue"]
-
-    # build set of coordinates for the original domain
-    y_coord = (
-        np.linspace(bottom, top - metadata["ypixelsize"], R.shape[2])
-        + metadata["ypixelsize"] / 2.0
-    )
-    x_coord = (
-        np.linspace(left, right - metadata["xpixelsize"], R.shape[3])
-        + metadata["xpixelsize"] / 2.0
-    )
-
-    # build set of coordinates for the new domain
-    y_coord_ = (
-        np.linspace(bottom_, top_ - metadata["ypixelsize"], R_.shape[2])
-        + metadata["ypixelsize"] / 2.0
-    )
-    x_coord_ = (
-        np.linspace(left_, right_ - metadata["xpixelsize"], R_.shape[3])
-        + metadata["xpixelsize"] / 2.0
-    )
-
-    # origin='upper' reverses the vertical axes direction
-    if metadata["yorigin"] == "upper":
-        y_coord = y_coord[::-1]
-        y_coord_ = y_coord_[::-1]
-
-    # extract original domain
-    idx_y = np.where(np.logical_and(y_coord < top_, y_coord > bottom_))[0]
-    idx_x = np.where(np.logical_and(x_coord < right_, x_coord > left_))[0]
-
-    # extract new domain
-    idx_y_ = np.where(np.logical_and(y_coord_ < top, y_coord_ > bottom))[0]
-    idx_x_ = np.where(np.logical_and(x_coord_ < right, x_coord_ > left))[0]
-
-    # compose the new array
-    R_[:, :, idx_y_[0] : (idx_y_[-1] + 1), idx_x_[0] : (idx_x_[-1] + 1)] = R[
-        :, :, idx_y[0] : (idx_y[-1] + 1), idx_x[0] : (idx_x[-1] + 1)
-    ]
-
-    # update coordinates
-    metadata["y1"] = bottom_
-    metadata["y2"] = top_
-    metadata["x1"] = left_
-    metadata["x2"] = right_
-
-    R_shape[-2] = R_.shape[-2]
-    R_shape[-1] = R_.shape[-1]
-
-    return R_.reshape(R_shape), metadata
-
-
-def square_domain(R, metadata, method="pad", inverse=False):
-    """Either pad or crop a field to obtain a square domain.
-
-    Parameters
-    ----------
-    R: array-like
-        Array of shape (m,n) or (t,m,n) containing the input fields.
-    metadata: dict
-        Metadata dictionary containing the x1, x2, y1, y2,
-        xpixelsize, ypixelsize,
-        attributes as described in the documentation of
-        :py:mod:`pysteps.io.importers`.
-    method: {'pad', 'crop'}, optional
-        Either pad or crop.
-        If pad, an equal number of zeros is added to both ends of its shortest
-        side in order to produce a square domain.
-        If crop, an equal number of pixels is removed
-        to both ends of its longest side in order to produce a square domain.
-        Note that the crop method involves an irreversible loss of data.
-    inverse: bool, optional
-        Perform the inverse method to recover the original domain shape.
-        After a crop, the inverse is performed by padding the field with zeros.
-
-    Returns
-    -------
-    R: array-like
-        the reshape dataset
-    metadata: dict
-        the metadata with updated attributes.
-    """
-
-    R = R.copy()
-    R_shape = np.array(R.shape)
-    metadata = metadata.copy()
-
-    if not inverse:
-
-        if len(R.shape) < 2:
-            raise ValueError("The number of dimension must be > 1")
-        if len(R.shape) == 2:
-            R = R[None, None, :]
-        if len(R.shape) == 3:
-            R = R[None, :]
-        if len(R.shape) > 4:
-            raise ValueError("The number of dimension must be <= 4")
-
-        if R.shape[2] == R.shape[3]:
-            return R.squeeze()
-
-        orig_dim = R.shape
-        orig_dim_n = orig_dim[0]
-        orig_dim_t = orig_dim[1]
-        orig_dim_y = orig_dim[2]
-        orig_dim_x = orig_dim[3]
-
-        if method == "pad":
-
-            new_dim = np.max(orig_dim[2:])
-            R_ = np.ones((orig_dim_n, orig_dim_t, new_dim, new_dim)) * R.min()
-
-            if orig_dim_x < new_dim:
-                idx_buffer = int((new_dim - orig_dim_x) / 2.0)
-                R_[:, :, :, idx_buffer : (idx_buffer + orig_dim_x)] = R
-                metadata["x1"] -= idx_buffer * metadata["xpixelsize"]
-                metadata["x2"] += idx_buffer * metadata["xpixelsize"]
-
-            elif orig_dim_y < new_dim:
-                idx_buffer = int((new_dim - orig_dim_y) / 2.0)
-                R_[:, :, idx_buffer : (idx_buffer + orig_dim_y), :] = R
-                metadata["y1"] -= idx_buffer * metadata["ypixelsize"]
-                metadata["y2"] += idx_buffer * metadata["ypixelsize"]
-
-        elif method == "crop":
-
-            new_dim = np.min(orig_dim[2:])
-            R_ = np.zeros((orig_dim_n, orig_dim_t, new_dim, new_dim))
-
-            if orig_dim_x > new_dim:
-                idx_buffer = int((orig_dim_x - new_dim) / 2.0)
-                R_ = R[:, :, :, idx_buffer : (idx_buffer + new_dim)]
-                metadata["x1"] += idx_buffer * metadata["xpixelsize"]
-                metadata["x2"] -= idx_buffer * metadata["xpixelsize"]
-
-            elif orig_dim_y > new_dim:
-                idx_buffer = int((orig_dim_y - new_dim) / 2.0)
-                R_ = R[:, :, idx_buffer : (idx_buffer + new_dim), :]
-                metadata["y1"] += idx_buffer * metadata["ypixelsize"]
-                metadata["y2"] -= idx_buffer * metadata["ypixelsize"]
-
-        else:
-            raise ValueError("Unknown type")
-
-        metadata["orig_domain"] = (orig_dim_y, orig_dim_x)
-        metadata["square_method"] = method
-
-        R_shape[-2] = R_.shape[-2]
-        R_shape[-1] = R_.shape[-1]
-
-        return R_.reshape(R_shape), metadata
-
-    elif inverse:
-
-        if len(R.shape) < 2:
-            raise ValueError("The number of dimension must be > 2")
-        if len(R.shape) == 2:
-            R = R[None, None, :]
-        if len(R.shape) == 3:
-            R = R[None, :]
-        if len(R.shape) > 4:
-            raise ValueError("The number of dimension must be <= 4")
-
-        method = metadata.pop("square_method")
-        shape = metadata.pop("orig_domain")
-
-        if R.shape[2] == shape[0] and R.shape[3] == shape[1]:
-            return R.squeeze(), metadata
-
-        R_ = np.zeros((R.shape[0], R.shape[1], shape[0], shape[1]))
-
-        if method == "pad":
-
-            if R.shape[2] == shape[0]:
-                idx_buffer = int((R.shape[3] - shape[1]) / 2.0)
-                R_ = R[:, :, :, idx_buffer : (idx_buffer + shape[1])]
-                metadata["x1"] += idx_buffer * metadata["xpixelsize"]
-                metadata["x2"] -= idx_buffer * metadata["xpixelsize"]
-
-            elif R.shape[3] == shape[1]:
-                idx_buffer = int((R.shape[2] - shape[0]) / 2.0)
-                R_ = R[:, :, idx_buffer : (idx_buffer + shape[0]), :]
-                metadata["y1"] += idx_buffer * metadata["ypixelsize"]
-                metadata["y2"] -= idx_buffer * metadata["ypixelsize"]
-
-        elif method == "crop":
-
-            if R.shape[2] == shape[0]:
-                idx_buffer = int((shape[1] - R.shape[3]) / 2.0)
-                R_[:, :, :, idx_buffer : (idx_buffer + R.shape[3])] = R
-                metadata["x1"] -= idx_buffer * metadata["xpixelsize"]
-                metadata["x2"] += idx_buffer * metadata["xpixelsize"]
-
-            elif R.shape[3] == shape[1]:
-                idx_buffer = int((shape[0] - R.shape[2]) / 2.0)
-                R_[:, :, idx_buffer : (idx_buffer + R.shape[2]), :] = R
-                metadata["y1"] -= idx_buffer * metadata["ypixelsize"]
-                metadata["y2"] += idx_buffer * metadata["ypixelsize"]
-
-        R_shape[-2] = R_.shape[-2]
-        R_shape[-1] = R_.shape[-1]
-
-        return R_.reshape(R_shape), metadata
diff --git a/pysteps/utils/images.py b/pysteps/utils/images.py
index ba4ad4e16..d9ebbd4a6 100644
--- a/pysteps/utils/images.py
+++ b/pysteps/utils/images.py
@@ -11,7 +11,7 @@
     morph_opening
 """
 import numpy as np
-from numpy.ma.core import MaskedArray
+import xarray as xr
 
 from pysteps.exceptions import MissingOptionalDependency
 
@@ -23,28 +23,24 @@
     CV2_IMPORTED = False
 
 
-def morph_opening(input_image, thr, n):
-    """Filter out small scale noise on the image by applying a binary
-    morphological opening, that is, erosion followed by dilation.
-
-    .. _MaskedArray:\
-        https://docs.scipy.org/doc/numpy/reference/maskedarray.baseclass.html#numpy.ma.MaskedArray
-
-    .. _ndarray:\
-    https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.html
+def morph_opening(input_array, thr, kernel):
+    """Filter out small scale noise on one image or sequence of images
+    by applying a binary morphological opening, that is, erosion followed
+    by dilation.
 
     Parameters
     ----------
-    input_image: ndarray_ or MaskedArray_
-        Array of shape (m, n) containing the input image.
+    input_array: xr.DataArray
+        Image or sequence of images as an array with dimensions (y, x) or (t, y, x).
     thr: float
-        The threshold used to convert the image into a binary image.
-    n: int
-        The structuring element size [pixels].
+        The threshold used to distinguish features from the background.
+    kernel: int or 2D array_like
+        The size of a rounded structuring element in number of pixels or
+        an arbitrary kernel as a two-dimensional array.
 
     Returns
     -------
-    input_image: ndarray_ or MaskedArray_
+    input_image: xr.DataArray
         Array of shape (m,n) containing the filtered image.
     """
     if not CV2_IMPORTED:
@@ -53,32 +49,39 @@ def morph_opening(input_image, thr, n):
             "routine but it is not installed"
         )
 
-    input_image = input_image.copy()
-
-    # Check if a MaskedArray is used. If not, mask the ndarray
-    to_ndarray = False
-    if not isinstance(input_image, MaskedArray):
-        to_ndarray = True
-        input_image = np.ma.masked_invalid(input_image)
-
-    np.ma.set_fill_value(input_image, input_image.min())
+    input_array = input_array.copy()
+    fill_value = input_array.attrs.get("zerovalue", input_array.min())
 
     # Convert to binary image
-    field_bin = np.ndarray.astype(input_image.filled() > thr, "uint8")
+    field_bin = (input_array > thr).astype(np.uint8)
 
-    # Build a structuring element of size n
-    kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (n, n))
+    # Build a structuring element
+    if isinstance(kernel, int):
+        kernel = cv2.getStructuringElement(
+            cv2.MORPH_ELLIPSE,
+            (kernel, kernel)
+        )
+    else:
+        kernel = np.array(kernel, dtype=np.uint8)
+
+    if "t" not in field_bin.dims:
+        field_bin = field_bin.expand_dims("t")
 
     # Apply morphological opening (i.e. erosion then dilation)
-    field_bin_out = cv2.morphologyEx(field_bin, cv2.MORPH_OPEN, kernel)
+    field_bin_out = xr.apply_ufunc(
+        cv2.morphologyEx,
+        field_bin.groupby("t"),
+        cv2.MORPH_OPEN,
+        kernel,
+    )
 
     # Build mask to be applied on the original image
-    mask = (field_bin - field_bin_out) > 0
+    mask_remove = (field_bin - field_bin_out) > 0
 
-    # Filter out small isolated pixels based on mask
-    input_image[mask] = np.nanmin(input_image)
+    if mask_remove.t.size == 1:
+        mask_remove = mask_remove.squeeze("t")
 
-    if to_ndarray:
-        input_image = np.array(input_image)
+    # Filter out small isolated pixels based on mask
+    input_image = input_array.where(~mask_remove, fill_value)
 
     return input_image
diff --git a/pysteps/utils/interface.py b/pysteps/utils/interface.py
index eebf84b1d..3972ed433 100644
--- a/pysteps/utils/interface.py
+++ b/pysteps/utils/interface.py
@@ -14,7 +14,6 @@
 from . import arrays
 from . import cleansing
 from . import conversion
-from . import dimension
 from . import fft
 from . import images
 from . import interpolate
@@ -59,20 +58,6 @@ def get_method(name, **kwargs):
     | reflectivity      |                                                     |
     +-------------------+-----------------------------------------------------+
 
-    Dimension methods:
-
-    +-------------------+-----------------------------------------------------+
-    |     Name          |              Description                            |
-    +===================+=====================================================+
-    |  accumulate       | aggregate fields in time                            |
-    +-------------------+-----------------------------------------------------+
-    |  clip             | resize the field domain by geographical coordinates |
-    +-------------------+-----------------------------------------------------+
-    |  square           | either pad or crop the data to get a square domain  |
-    +-------------------+-----------------------------------------------------+
-    |  upscale          | upscale the field                                   |
-    +-------------------+-----------------------------------------------------+
-
     FFT methods (wrappers to different implementations):
 
     +-------------------+-----------------------------------------------------+
@@ -154,8 +139,8 @@ def get_method(name, **kwargs):
 
     name = name.lower()
 
-    def donothing(R, metadata=None, *args, **kwargs):
-        return R.copy(), {} if metadata is None else metadata.copy()
+    def donothing(data_array, *args, **kwargs):
+        return data_array.copy()
 
     methods_objects = dict()
     methods_objects["none"] = donothing
@@ -175,12 +160,6 @@ def donothing(R, metadata=None, *args, **kwargs):
     methods_objects["dbz"] = conversion.to_reflectivity
     methods_objects["reflectivity"] = conversion.to_reflectivity
 
-    # dimension methods
-    methods_objects["accumulate"] = dimension.aggregate_fields_time
-    methods_objects["clip"] = dimension.clip_domain
-    methods_objects["square"] = dimension.square_domain
-    methods_objects["upscale"] = dimension.aggregate_fields_space
-
     # image processing methods
     methods_objects["morph_opening"] = images.morph_opening
 
@@ -198,10 +177,10 @@ def donothing(R, metadata=None, *args, **kwargs):
     # transformation methods
     methods_objects["boxcox"] = transformation.boxcox_transform
     methods_objects["box-cox"] = transformation.boxcox_transform
-    methods_objects["db"] = transformation.dB_transform
-    methods_objects["decibel"] = transformation.dB_transform
+    methods_objects["db"] = transformation.db_transform
+    methods_objects["decibel"] = transformation.db_transform
     methods_objects["log"] = transformation.boxcox_transform
-    methods_objects["nqt"] = transformation.NQ_transform
+    methods_objects["nqt"] = transformation.nq_transform
     methods_objects["sqrt"] = transformation.sqrt_transform
 
     # FFT methods
diff --git a/pysteps/utils/interpolate.py b/pysteps/utils/interpolate.py
index 7f087e90b..7e936d3b3 100644
--- a/pysteps/utils/interpolate.py
+++ b/pysteps/utils/interpolate.py
@@ -12,37 +12,38 @@
     rbfinterp2d
 """
 
-import warnings
-
 import numpy as np
-from scipy.interpolate import Rbf
+import xarray as xr
+from pandas import MultiIndex
 from scipy.spatial import cKDTree
 from scipy.spatial.distance import cdist
 
-
 from pysteps.decorators import memoize, prepare_interpolator
+from pysteps.exceptions import MissingOptionalDependency
 
+# only scipy>=1.7
+try:
+    from scipy.interpolate import RBFInterpolator
 
-@prepare_interpolator()
-def idwinterp2d(
-    xy_coord, values, xgrid, ygrid, power=0.5, k=20, dist_offset=0.5, **kwargs
-):
-    """Inverse distance weighting interpolation of a sparse (multivariate) array.
+    RBF_IMPORTED = True
+except ImportError:
+    from scipy import __version__ as scipy_version
+    RBF_IMPORTED = False
 
-    .. _ndarray:\
-    https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.html
+
+@prepare_interpolator()
+def idwinterp2d(sparse_data, xgrid, ygrid, power=0.5, k=20, dist_offset=0.5, **kwargs):
+    """
+    Gridding of sparse data points with inverse distance weighting (IDW).
+    Samples with missing values or coordinates are dropped.
 
     Parameters
     ----------
-    xy_coord: ndarray_
-        Array of shape (n, 2) containing the coordinates of the data points in
-        a 2-dimensional space.
-    values: ndarray_
-        Array of shape (n) or (n, m) containing the values of the data points,
-        where *n* is the number of data points and *m* the number of co-located
-        variables. All elements in ``values`` are required to be finite.
-    xgrid, ygrid: ndarray_
-        1-D arrays representing the coordinates of the 2-D output grid.
+    sparse_data: xr.DataArray
+        The sparse dataset with dimension "sample" and coordinates ("x", "y").
+    xgrid, ygrid: array_like
+        1-D arrays representing the coordinates of the 2-D target grid at which
+        to interpolate data.
     power: positive float, optional
         The power parameter used to compute the distance weights as
         ``weight = distance ** (-power)``.
@@ -59,18 +60,11 @@ def idwinterp2d(
 
     Returns
     -------
-    output_array: ndarray_
-        The interpolated field(s) having shape (``ygrid.size``, ``xgrid.size``)
-        or (*m*, ``ygrid.size``, ``xgrid.size``).
+    output_grid: xr.DataArray
+        The dataset interpolated on the target grid.
     """
-    if values.ndim == 1:
-        nvar = 1
-        values = values[:, None]
-
-    elif values.ndim == 2:
-        nvar = values.shape[1]
-
-    npoints = values.shape[0]
+    npoints = sparse_data.sizes["sample"]
+    xy_coord = np.column_stack((sparse_data.x.values, sparse_data.y.values))
 
     # generate the target grid
     xgridv, ygridv = np.meshgrid(xgrid, ygrid)
@@ -89,6 +83,7 @@ def idwinterp2d(
         inds = np.arange(npoints)[None, :] * np.ones((gridv.shape[0], npoints)).astype(
             int
         )
+    inds = xr.DataArray(inds, dims=("grid", "neighbor"))
 
     # convert geographical distances to number of pixels
     x_res = np.gradient(xgrid)
@@ -100,41 +95,36 @@ def idwinterp2d(
     dist += dist_offset  # avoid zero distances
     weights = 1 / np.power(dist, power)
     weights = weights / np.sum(weights, axis=1, keepdims=True)
+    weights = xr.DataArray(weights, dims=("grid", "neighbor"))
 
     # interpolate
-    output_array = np.sum(
-        values[inds, :] * weights[..., None],
-        axis=1,
-    )
+    output_array = (sparse_data[inds] * weights).sum("neighbor")
+    output_array = output_array.astype(sparse_data.dtype)
 
-    # reshape to final grid size
-    output_array = output_array.reshape(ygrid.size, xgrid.size, nvar)
+    # assign multi-index coordinate for the grid
+    grid = MultiIndex.from_product([ygrid, xgrid], names=("y", "x"))
+    output_array = output_array.assign_coords({"grid": grid})
 
-    return np.moveaxis(output_array, -1, 0).squeeze()
+    return output_array.unstack("grid")  # reshape as grid
 
 
 @prepare_interpolator()
-def rbfinterp2d(xy_coord, values, xgrid, ygrid, **kwargs):
-    """Radial basis function interpolation of a sparse (multivariate) array.
+def rbfinterp2d(sparse_data, xgrid, ygrid, **kwargs):
+    """Gridding of sparse data points with radial basis functions (RBF).
 
-    .. _ndarray:\
-    https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.html
-    .. _`scipy.interpolate.Rbf`:\
-    https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.Rbf.html
+    .. _`scipy.interpolate.RBFInterpolator`:\
+    https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RBFInterpolator.html
 
-    This method wraps the `scipy.interpolate.Rbf`_ class.
+    This method wraps the `scipy.interpolate.RBFInterpolator`_ class available
+    from scipy's version 1.7 onwards.
 
     Parameters
     ----------
-    xy_coord: ndarray_
-        Array of shape (n, 2) containing the coordinates of the data points in
-        a 2-dimensional space.
-    values: ndarray_
-        Array of shape (n) or (n, m) containing the values of the data points,
-        where *n* is the number of data points and *m* the number of co-located
-        variables. All values in ``values`` are required to be finite.
-    xgrid, ygrid: ndarray_
-        1-D arrays representing the coordinates of the 2-D output grid.
+    sparse_data: xr.DataArray
+        The sparse dataset with dimension "sample" and coordinates ("x", "y").
+    xgrid, ygrid: array_like
+        1-D arrays representing the coordinates of the 2-D target grid at which
+        to interpolate data.
 
     Other Parameters
     ----------------
@@ -143,38 +133,46 @@ def rbfinterp2d(xy_coord, values, xgrid, ygrid, **kwargs):
 
     Returns
     -------
-    output_array: ndarray_
-        The interpolated field(s) having shape (``ygrid.size``, ``xgrid.size``)
-        or (*m*, ``ygrid.size``, ``xgrid.size``).
+    output_grid: xr.DataArray
+        The dataset interpolated on the target grid.
     """
-    deprecated_args = ["rbfunction", "k"]
-    deprecated_args = [arg for arg in deprecated_args if arg in list(kwargs.keys())]
-    if deprecated_args:
-        warnings.warn(
-            "rbfinterp2d: The following keyword arguments are deprecated:\n"
-            + str(deprecated_args),
-            DeprecationWarning,
+    if not RBF_IMPORTED:
+        raise MissingOptionalDependency(
+            "Scipy's RBFInterpolate could not be imported. "
+            "Check your scipy installation: "
+            f"we found version {scipy_version} (should be >=1.7)"
         )
 
-    if values.ndim == 1:
-        kwargs["mode"] = "1-D"
-    else:
-        kwargs["mode"] = "N-D"
-
+    # generate the target grid
     xgridv, ygridv = np.meshgrid(xgrid, ygrid)
-    rbfi = _Rbf_cached(*np.split(xy_coord, xy_coord.shape[1], 1), values, **kwargs)
-    output_array = rbfi(xgridv, ygridv)
+    gridv = np.column_stack((xgridv.ravel(), ygridv.ravel()))
+
+    # interpolate
+    rbfi = _Rbf_cached(sparse_data, **kwargs)
+    output_array = xr.DataArray(
+        rbfi(gridv),
+        dims=("grid", "variable"),
+        coords=sparse_data.drop_vars(("x", "y", "sample"), errors="ignore").coords,
+        attrs=sparse_data.attrs,
+    )
+    output_array = output_array.astype(sparse_data.dtype)
+
+    # assign multi-index coordinate for the grid
+    grid = MultiIndex.from_product([ygrid, xgrid], names=("y", "x"))
+    output_array = output_array.assign_coords({"grid": grid})
 
-    return np.moveaxis(output_array, -1, 0).squeeze()
+    return output_array.unstack("grid")  # reshape as grid
 
 
 @memoize()
 def _cKDTree_cached(*args, **kwargs):
     """Add LRU cache to cKDTree class."""
-    return cKDTree(*args)
+    return cKDTree(*args, **kwargs)
 
 
 @memoize()
-def _Rbf_cached(*args, **kwargs):
+def _Rbf_cached(sparse_data, **kwargs):
     """Add LRU cache to Rbf class."""
-    return Rbf(*args, **kwargs)
+    data_coords = np.column_stack((sparse_data.x, sparse_data.y))
+    data_values = sparse_data.values
+    return RBFInterpolator(data_coords, data_values, **kwargs)
diff --git a/pysteps/utils/transformation.py b/pysteps/utils/transformation.py
index 1c389b748..155d4f9c2 100644
--- a/pysteps/utils/transformation.py
+++ b/pysteps/utils/transformation.py
@@ -8,65 +8,61 @@
 .. autosummary::
     :toctree: ../generated/
 
+    back_transform
     boxcox_transform
-    dB_transform
-    NQ_transform
+    db_transform
+    nq_transform
     sqrt_transform
 """
+from functools import partial
+from typing import Any, Optional
 
 import numpy as np
-import scipy.stats as scipy_stats
-import warnings
-from scipy.interpolate import interp1d
+import xarray as xr
+from scipy.stats import norm
 
-warnings.filterwarnings(
-    "ignore", category=RuntimeWarning
-)  # To deactivate warnings for comparison operators with NaNs
+from pysteps.accessors import dataarray_utils
 
 
+@dataarray_utils
 def boxcox_transform(
-    R, metadata=None, Lambda=None, threshold=None, zerovalue=None, inverse=False
-):
+    data_array: xr.DataArray,
+    boxcox_lambda: float = 0.0,
+    offset: float = 0.01,
+    inverse: bool = False,
+) -> xr.DataArray:
     """The one-parameter Box-Cox transformation.
 
     The Box-Cox transform is a well-known power transformation introduced by
     Box and Cox (1964). In its one-parameter version, the Box-Cox transform
-    takes the form T(x) = ln(x) for Lambda = 0,
-    or T(x) = (x**Lambda - 1)/Lambda otherwise.
+    takes the form T(x) = ln(x) for boxcox_lambda = 0,
+    or T(x) = (x ** boxcox_lambda - 1) / boxcox_lambda otherwise.
 
-    Default parameters will produce a log transform (i.e. Lambda=0).
+    Default parameters will produce a log transform (i.e. boxcox_lambda=0).
 
     Parameters
     ----------
-    R: array-like
-        Array of any shape to be transformed.
-    metadata: dict, optional
-        Metadata dictionary containing the transform, zerovalue and threshold
-        attributes as described in the documentation of
-        :py:mod:`pysteps.io.importers`.
-    Lambda: float, optional
-        Parameter Lambda of the Box-Cox transformation.
+    data_array: xr.DataArray
+        Array of any shape to be (back-)transformed.
+
+    boxcox_lambda: float
+        Parameter lambda of the Box-Cox transformation.
         It is 0 by default, which produces the log transformation.
 
-        Choose Lambda < 1 for positively skewed data, Lambda > 1 for negatively
-        skewed data.
-    threshold: float, optional
-        The value that is used for thresholding with the same units as R.
-        If None, the threshold contained in metadata is used.
-        If no threshold is found in the metadata,
-        a value of 0.1 is used as default.
-    zerovalue: float, optional
-        The value to be assigned to no rain pixels as defined by the threshold.
-        It is equal to the threshold - 1 by default.
-    inverse: bool, optional
+        Choose boxcox_lambda < 1 for positively skewed data, boxcox_lambda > 1
+        for negatively skewed data.
+
+    offset: float
+        A small constant offset which prevents the division by zero or the
+        computation of log(0).
+
+    inverse: bool
         If set to True, it performs the inverse transform. False by default.
 
     Returns
     -------
-    R: array-like
-        Array of any shape containing the (back-)transformed units.
-    metadata: dict
-        The metadata with updated attributes.
+    data_array: xr.DataArray
+        DataArray containing the (back-)transformed units.
 
     References
     ----------
@@ -74,198 +70,132 @@ def boxcox_transform(
     of the Royal Statistical Society: Series B (Methodological), 26: 211-243.
     doi:10.1111/j.2517-6161.1964.tb00553.x
     """
-
-    R = R.copy()
-
-    if metadata is None:
-        if inverse:
-            metadata = {"transform": "BoxCox"}
-        else:
-            metadata = {"transform": None}
-
-    else:
-        metadata = metadata.copy()
+    data_array = data_array.copy()
+    attrs = data_array.attrs
 
     if not inverse:
 
-        if metadata["transform"] == "BoxCox":
-            return R, metadata
-
-        if Lambda is None:
-            Lambda = metadata.get("BoxCox_lambda", 0.0)
-
-        if threshold is None:
-            threshold = metadata.get("threshold", 0.1)
-
-        zeros = R < threshold
-
-        # Apply Box-Cox transform
-        if Lambda == 0.0:
-            R[~zeros] = np.log(R[~zeros])
-            threshold = np.log(threshold)
-
+        data_array = data_array.pysteps.back_transform()
+        data_array += offset
+        if boxcox_lambda == 0.0:
+            data_array = np.log(data_array)
         else:
-            R[~zeros] = (R[~zeros] ** Lambda - 1) / Lambda
-            threshold = (threshold ** Lambda - 1) / Lambda
-
-        # Set value for zeros
-        if zerovalue is None:
-            zerovalue = threshold - 1  # TODO: set to a more meaningful value
-        R[zeros] = zerovalue
-
-        metadata["transform"] = "BoxCox"
-        metadata["BoxCox_lambda"] = Lambda
-        metadata["zerovalue"] = zerovalue
-        metadata["threshold"] = threshold
-
-    elif inverse:
+            data_array = (data_array ** boxcox_lambda - 1) / boxcox_lambda
+
+        attrs.update(
+            {
+                "transform": "BoxCox",
+                "boxcox_lambda": boxcox_lambda,
+                "offset": offset,
+            }
+        )
+        data_array.attrs = attrs
 
-        if metadata["transform"] not in ["BoxCox", "log"]:
-            return R, metadata
+    else:
 
-        if Lambda is None:
-            Lambda = metadata.pop("BoxCox_lambda", 0.0)
-        if threshold is None:
-            threshold = metadata.get("threshold", -10.0)
-        if zerovalue is None:
-            zerovalue = 0.0
+        if attrs["transform"].lower() not in ["boxcox", "log"]:
+            raise ValueError(f"cannot back-transform from {attrs['transform']}")
 
-        # Apply inverse Box-Cox transform
-        if Lambda == 0.0:
-            R = np.exp(R)
-            threshold = np.exp(threshold)
+        boxcox_lambda = attrs.pop("boxcox_lambda", boxcox_lambda)
 
+        if boxcox_lambda == 0.0:
+            data_array = np.exp(data_array)
         else:
-            R = np.exp(np.log(Lambda * R + 1) / Lambda)
-            threshold = np.exp(np.log(Lambda * threshold + 1) / Lambda)
+            data_array = np.log(boxcox_lambda * data_array + 1) / boxcox_lambda
+            data_array = np.exp(data_array)
 
-        R[R < threshold] = zerovalue
+        data_array -= attrs.pop("offset", offset)
+        attrs.update({"transform": None})
+        data_array.attrs = attrs
 
-        metadata["transform"] = None
-        metadata["zerovalue"] = zerovalue
-        metadata["threshold"] = threshold
+    return data_array
 
-    return R, metadata
 
-
-def dB_transform(R, metadata=None, threshold=None, zerovalue=None, inverse=False):
-    """Methods to transform precipitation intensities to/from dB units.
+@dataarray_utils
+def db_transform(
+    data_array: xr.DataArray, offset: float = 0.01, inverse: bool = False
+) -> xr.DataArray:
+    """Transform to/from decibel (dB).
 
     Parameters
     ----------
-    R: array-like
+    data_array: xr.DataArray
         Array of any shape to be (back-)transformed.
-    metadata: dict, optional
-        Metadata dictionary containing the transform, zerovalue and threshold
-        attributes as described in the documentation of
-        :py:mod:`pysteps.io.importers`.
-    threshold: float, optional
-        Optional value that is used for thresholding with the same units as R.
-        If None, the threshold contained in metadata is used.
-        If no threshold is found in the metadata,
-        a value of 0.1 is used as default.
-    zerovalue: float, optional
-        The value to be assigned to no rain pixels as defined by the threshold.
-        It is equal to the threshold - 1 by default.
-    inverse: bool, optional
+
+    offset: float
+        A small constant offset which prevents the division by zero or the
+        computation of log(0).
+
+    inverse: bool
         If set to True, it performs the inverse transform. False by default.
 
     Returns
     -------
-    R: array-like
-        Array of any shape containing the (back-)transformed units.
-    metadata: dict
-        The metadata with updated attributes.
+    data_array: xr.DataArray
+        DataArray containing the (back-)transformed units.
     """
+    data_array = data_array.copy()
+    attrs = data_array.attrs
 
-    R = R.copy()
-
-    if metadata is None:
-        if inverse:
-            metadata = {"transform": "dB"}
-        else:
-            metadata = {"transform": None}
-
-    else:
-        metadata = metadata.copy()
-
-    # to dB units
     if not inverse:
+        data_array = data_array.pysteps.back_transform()
+        data_array += offset
+        data_array = 10.0 * np.log10(data_array)
+
+        attrs.update(
+            {
+                "transform": "dB",
+                "offset": offset,
+            }
+        )
+        data_array.attrs = attrs
 
-        if metadata["transform"] == "dB":
-            return R, metadata
-
-        if threshold is None:
-            threshold = metadata.get("threshold", 0.1)
-
-        zeros = R < threshold
-
-        # Convert to dB
-        R[~zeros] = 10.0 * np.log10(R[~zeros])
-        threshold = 10.0 * np.log10(threshold)
-
-        # Set value for zeros
-        if zerovalue is None:
-            zerovalue = threshold - 5  # TODO: set to a more meaningful value
-        R[zeros] = zerovalue
-
-        metadata["transform"] = "dB"
-        metadata["zerovalue"] = zerovalue
-        metadata["threshold"] = threshold
-
-        return R, metadata
-
-    # from dB units
-    elif inverse:
-
-        if metadata["transform"] != "dB":
-            return R, metadata
-
-        if threshold is None:
-            threshold = metadata.get("threshold", -10.0)
-        if zerovalue is None:
-            zerovalue = 0.0
+    else:
 
-        R = 10.0 ** (R / 10.0)
-        threshold = 10.0 ** (threshold / 10.0)
-        R[R < threshold] = zerovalue
+        if not attrs["transform"].lower() == "db":
+            raise ValueError(f"cannot back-transform from {attrs['transform']}")
 
-        metadata["transform"] = None
-        metadata["threshold"] = threshold
-        metadata["zerovalue"] = zerovalue
+        data_array = 10.0 ** (data_array / 10.0)
+        data_array -= attrs.pop("offset", offset)
+        attrs.update({"transform": None})
+        data_array.attrs = attrs
 
-        return R, metadata
+    return data_array
 
 
-def NQ_transform(R, metadata=None, inverse=False, **kwargs):
+@dataarray_utils
+def nq_transform(
+    data_array: xr.DataArray,
+    nq_a: float = 0.0,
+    template: Optional[Any] = None,
+    inverse: bool = False,
+) -> xr.DataArray:
     """The normal quantile transformation as in Bogner et al (2012).
     Zero rain vales are set to zero in norm space.
 
     Parameters
     ----------
-    R: array-like
-        Array of any shape to be transformed.
-    metadata: dict, optional
-        Metadata dictionary containing the transform, zerovalue and threshold
-        attributes as described in the documentation of
-        :py:mod:`pysteps.io.importers`.
-    inverse: bool, optional
-        If set to True, it performs the inverse transform. False by default.
+    data_array: xr.DataArray
+        Array of any shape to be (back-)transformed.
+
+    nq_a: float
+        The offset fraction to be used for plotting positions; typically in (0, 1).
+        The default is 0, that is, it spaces the points evenly in the uniform
+        distribution.
 
-    Other Parameters
-    ----------------
-    a: float, optional
-        The offset fraction to be used for plotting positions;
-        typically in (0,1).
-        The default is 0., that is, it spaces the points evenly in the uniform
+    template: array_like, optional
+        Array of any shape containing the samples to build the target empirical
         distribution.
+        Required if inverse is True.
+
+    inverse: bool
+        If set to True, it performs the inverse transform. False by default.
+
 
     Returns
     -------
-    R: array-like
-        Array of any shape containing the (back-)transformed units.
-    metadata: dict
-        The metadata with updated attributes.
+    data_array: xr.DataArray
+        DataArray containing the (back-)transformed units.
 
     References
     ----------
@@ -274,108 +204,139 @@ def NQ_transform(R, metadata=None, inverse=False, **kwargs):
     Hydrol. Earth Syst. Sci., 16, 1085-1094,
     https://doi.org/10.5194/hess-16-1085-2012, 2012.
     """
+    data_array = data_array.copy()
+    attrs = data_array.attrs
+
+    if not inverse:
 
-    # defaults
-    a = kwargs.get("a", 0.0)
+        data_array = data_array.pysteps.back_transform()
 
-    R = R.copy()
-    shape0 = R.shape
-    R = R.ravel().astype(float)
-    idxNan = np.isnan(R)
-    R_ = R[~idxNan]
+        # Resolve ties at random
+        data_array += np.random.random(data_array.shape) / 1e10
 
-    if metadata is None:
-        if inverse:
-            metadata = {"transform": "NQT"}
-        else:
-            metadata = {"transform": None}
-        metadata["zerovalue"] = np.min(R_)
+        # Compute quantiles
+        data_array = data_array.stack(dummy=(data_array.dims))
+        n = int(data_array.notnull().sum())
+        data_quantiles = (data_array.rank("dummy") - nq_a) / (n + 1 - 2 * nq_a)
+        data_quantiles = data_quantiles.unstack("dummy")
 
-    else:
-        metadata = metadata.copy()
+        # Normal quantile transform
+        data_array = xr.apply_ufunc(norm.ppf, data_quantiles)
 
-    if not inverse:
-        # Plotting positions
-        # https://en.wikipedia.org/wiki/Q%E2%80%93Q_plot#Plotting_position
-        n = R_.size
-        Rpp = ((np.arange(n) + 1 - a) / (n + 1 - 2 * a)).reshape(R_.shape)
+        attrs.update(
+            {
+                "transform": "NQ",
+                "unit": None,
+            }
+        )
+        data_array.attrs = attrs
 
-        # NQ transform
-        Rqn = scipy_stats.norm.ppf(Rpp)
-        R__ = np.interp(R_, R_[np.argsort(R_)], Rqn)
+    else:
 
-        # set zero rain to 0 in norm space
-        R__[R[~idxNan] == metadata["zerovalue"]] = 0
+        if not attrs["transform"].lower() == "nq":
+            raise ValueError(f"cannot back-transform from {attrs['transform']}")
 
-        # build inverse transform
-        metadata["inqt"] = interp1d(
-            Rqn, R_[np.argsort(R_)], bounds_error=False, fill_value=(R_.min(), R_.max())
-        )
+        if template is None:
+            raise ValueError("inverse NQ transform needs a template")
 
-        metadata["transform"] = "NQT"
-        metadata["zerovalue"] = 0
-        metadata["threshold"] = R__[R__ > 0].min()
+        try:
+            template_unit = template.attrs.get("unit", None)
+        except AttributeError:
+            template_unit = None
 
-    else:
-        f = metadata.pop("inqt")
-        R__ = f(R_)
-        metadata["transform"] = None
-        metadata["zerovalue"] = R__.min()
-        metadata["threshold"] = R__[R__ > R__.min()].min()
+        # Flatten and sort template
+        template = np.array(template).flatten()
+        template.sort()
+        n = int(np.isfinite(template).sum())
+        template = np.concatenate((template, np.array([np.nan])))
 
-    R[~idxNan] = R__
+        # Inverse of ECDF
+        data_array = data_array.stack(dummy=(data_array.dims))
+        data_quantiles = xr.apply_ufunc(norm.cdf, data_array)
+        index = np.floor((n - 1) * data_quantiles)
+        index = index.fillna(n).astype("uint32")
+        data_array = xr.DataArray(
+            data=template[index], coords=data_array.coords, dims=data_array.dims
+        ).unstack("dummy")
 
-    return R.reshape(shape0), metadata
+        attrs.update({"transform": None, "unit": template_unit})
+        data_array.attrs = attrs
 
+    return data_array
 
-def sqrt_transform(R, metadata=None, inverse=False, **kwargs):
+
+@dataarray_utils
+def sqrt_transform(data_array: xr.DataArray, inverse: bool = False) -> xr.DataArray:
     """Square-root transform.
 
     Parameters
     ----------
-    R: array-like
-        Array of any shape to be transformed.
-    metadata: dict, optional
-        Metadata dictionary containing the transform, zerovalue and threshold
-        attributes as described in the documentation of
-        :py:mod:`pysteps.io.importers`.
+    data_array: xr.DataArray
+        Array of any shape to be (back-)transformed.
+
     inverse: bool, optional
         If set to True, it performs the inverse transform. False by default.
 
     Returns
     -------
-    R: array-like
-        Array of any shape containing the (back-)transformed units.
-    metadata: dict
-        The metadata with updated attributes.
-
+    data_array: xr.DataArray
+        DataArray containing the (back-)transformed units.
     """
+    data_array = data_array.copy()
+    attrs = data_array.attrs
 
-    R = R.copy()
+    if not inverse:
+        data_array = data_array.pysteps.back_transform()
+        data_array = np.sqrt(data_array)
+
+        attrs.update(
+            {
+                "transform": "sqrt",
+            }
+        )
+        data_array.attrs = attrs
 
-    if metadata is None:
-        if inverse:
-            metadata = {"transform": "sqrt"}
-        else:
-            metadata = {"transform": None}
-        metadata["zerovalue"] = np.nan
-        metadata["threshold"] = np.nan
     else:
-        metadata = metadata.copy()
+        if not attrs["transform"].lower() == "sqrt":
+            raise ValueError(f"cannot back-transform from {attrs['transform']}")
 
-    if not inverse:
-        # sqrt transform
-        R = np.sqrt(R)
+        data_array = data_array ** 2.0
+        attrs.update({"transform": None})
+        data_array.attrs = attrs
 
-        metadata["transform"] = "sqrt"
-        metadata["zerovalue"] = np.sqrt(metadata["zerovalue"])
-        metadata["threshold"] = np.sqrt(metadata["threshold"])
-    else:
-        # inverse sqrt transform
-        R = R ** 2
+    return data_array
+
+
+_methods = dict()
+_methods[None] = lambda x: x.copy()
+_methods["dB"] = partial(db_transform, inverse=True)
+_methods["boxcox"] = partial(boxcox_transform, inverse=True)
+_methods["log"] = partial(boxcox_transform, inverse=True, boxcox_lambda=0.0)
+_methods["nq"] = partial(nq_transform, inverse=True)
+_methods["sqrt"] = partial(sqrt_transform, inverse=True)
 
-        metadata["transform"] = None
-        metadata["zerovalue"] = metadata["zerovalue"] ** 2
-        metadata["threshold"] = metadata["threshold"] ** 2
 
-    return R, metadata
+@dataarray_utils
+def back_transform(data_array: xr.DataArray) -> xr.DataArray:
+    """
+    Remove any existing data transformation.
+
+    Parameters
+    ----------
+    data_array: xr.DataArray
+        Array of any shape to be back-transformed.
+
+    Returns
+    -------
+    data_array: xr.DataArray
+        DataArray containing back-transformed units.
+    """
+    transform = data_array.attrs.get("transform")
+    try:
+        inverse_method = _methods[transform]
+    except KeyError:
+        raise ValueError(
+            f"Unknown transformation {transform}. "
+            f"The available methods are: {list(_methods.keys())}"
+        )
+    return inverse_method(data_array)
diff --git a/pysteps/verification/ensscores.py b/pysteps/verification/ensscores.py
index 7dae71df9..8836f7252 100644
--- a/pysteps/verification/ensscores.py
+++ b/pysteps/verification/ensscores.py
@@ -143,8 +143,8 @@ def rankhist(X_f, X_o, X_min=None, normalize=True):
         the bin counts sum to one.
     """
 
-    X_f = X_f.copy()
-    X_o = X_o.copy()
+    X_f = np.copy(X_f)
+    X_o = np.copy(X_o)
     num_ens_members = X_f.shape[0]
     rhist = rankhist_init(num_ens_members, X_min)
     rankhist_accum(rhist, X_f, X_o)
diff --git a/pysteps/verification/spatialscores.py b/pysteps/verification/spatialscores.py
index 809c7569c..67d9a57b6 100644
--- a/pysteps/verification/spatialscores.py
+++ b/pysteps/verification/spatialscores.py
@@ -579,9 +579,9 @@ def fss_accum(fss, X_f, X_o):
         message += " having the same shape"
         raise ValueError(message)
 
-    X_f = X_f.copy()
+    X_f = np.copy(X_f)
     X_f[~np.isfinite(X_f)] = fss["thr"] - 1
-    X_o = X_o.copy()
+    X_o = np.copy(X_o)
     X_o[~np.isfinite(X_o)] = fss["thr"] - 1
 
     # Convert to binary fields with the given intensity threshold


Test patch:
diff --git a/pysteps/tests/helpers.py b/pysteps/tests/helpers.py
index 86c59b663..30b04e5be 100644
--- a/pysteps/tests/helpers.py
+++ b/pysteps/tests/helpers.py
@@ -9,10 +9,8 @@
 import numpy as np
 import pytest
 
-import pysteps as stp
-from pysteps import io, rcparams
-from pysteps.decorators import _to_xarray, _xarray2legacy
-from pysteps.utils import aggregate_fields_space
+from pysteps import io, utils, rcparams
+from pysteps.decorators import _xarray2legacy
 
 _reference_dates = dict()
 _reference_dates["bom"] = datetime(2018, 6, 16, 10, 0)
@@ -27,86 +25,70 @@
 def get_precipitation_fields(
     num_prev_files=0,
     num_next_files=0,
-    return_raw=False,
-    metadata=False,
-    upscale=None,
     source="mch",
-    legacy=True,
-    log_transform=True,
+    convert_to=None,
+    transform_to=None,
+    filled=False,
+    coarsen=None,
     clip=None,
-    **importer_kwargs,
+    importer_kwargs=None,
+    **kwargs,
 ):
     """
-        Get a precipitation field from the archive to be used as reference.
-
-        Source: bom
-        Reference time: 2018/06/16 10000 UTC
-
-        Source: fmi
-        Reference time: 2016/09/28 1600 UTC
-
-        Source: knmi
-        Reference time: 2010/08/26 0000 UTC
-
-        Source: mch
-        Reference time: 2015/05/15 1630 UTC
-
-        Source: opera
-        Reference time: 2018/08/24 1800 UTC
-
-        Source: saf
-        Reference time: 2018/06/01 0700 UTC
-
-        Source: mrms
-        Reference time: 2019/06/10 0000 UTC
-
-        Parameters
-        ----------
-
-        num_prev_files: int, optional
-            Number of previous times (files) to return with respect to the
-            reference time.
-
-        num_next_files: int, optional
-            Number of future times (files) to return with respect to the
-            reference time.
-
-        return_raw: bool, optional
-            Do not preprocess the precipitation fields. False by default.
-            The pre-processing steps are: 1) Convert to mm/h,
-            2) Mask invalid values, 3) Log-transform the data [dBR].
-
-    <<<<<<< HEAD
-    =======
-        metadata: bool, optional
-            If True, also return file metadata.
-
-        clip: scalars (left, right, bottom, top), optional
-            The extent of the bounding box in data coordinates to be used to clip
-            the data.
-
-    >>>>>>> master
-        upscale: float or None, optional
-            Upscale fields in space during the pre-processing steps.
-            If it is None, the precipitation field is not modified.
-            If it is a float, represents the length of the space window that is
-            used to upscale the fields.
-
-        source: {"bom", "fmi" , "knmi", "mch", "opera", "saf", "mrms"}, optional
-            Name of the data source to be used.
-
-        log_transform: bool
-            Whether to transform the output to dB.
-
-        Other Parameters
-        ----------------
-
-        importer_kwargs : dict
-            Additional keyword arguments passed to the importer.
-
-        Returns
-        -------
-        data_array : xr.DataArray
+    Get a precipitation field from the archive to be used as reference.
+
+    Source: bom
+    Reference time: 2018/06/16 10000 UTC
+
+    Source: fmi
+    Reference time: 2016/09/28 1600 UTC
+
+    Source: knmi
+    Reference time: 2010/08/26 0000 UTC
+
+    Source: mch
+    Reference time: 2015/05/15 1630 UTC
+
+    Source: opera
+    Reference time: 2018/08/24 1800 UTC
+
+    Source: saf
+    Reference time: 2018/06/01 0700 UTC
+
+    Source: mrms
+    Reference time: 2019/06/10 0000 UTC
+
+    Parameters
+    ----------
+
+    num_prev_files: int, optional
+        Number of previous times (files) to return with respect to the
+        reference time.
+    num_next_files: int, optional
+        Number of future times (files) to return with respect to the
+        reference time.
+    source: {"bom", "fmi" , "knmi", "mch", "opera", "saf", "mrms"}
+        Name of the data source to be used.
+    convert_to: {"mm/h", "mm", "Z", "dBZ"}, optional
+        Convert to given units.
+    transform_to: {"boxcox", "log", "db", "nq", "sqrt"}, optional
+        Tranform data.
+    filled: bool
+        Whether to fill all the missing values.
+    clip: scalars (left, right, bottom, top), optional
+        The extent of the bounding box in data coordinates to be used to clip
+        the data.
+    coarsen: float, optional
+        Upscale fields in space during the pre-processing steps.
+        If it is None, the precipitation field is not modified.
+        If it is a float, represents the the window size that is used to
+        upscale the fields.
+    importer_kwargs : dict, optional
+        Additional keyword arguments passed to the importer.
+
+    Returns
+    -------
+    data_array : xr.DataArray
     """
 
     if source == "bom":
@@ -139,6 +121,9 @@ def get_precipitation_fields(
             f"{str(list(_reference_dates.keys()))}"
         )
 
+    if importer_kwargs is None:
+        importer_kwargs = {}
+
     data_source = rcparams.data_sources[source]
     root_path = data_source["root_path"]
     path_fmt = data_source["path_fmt"]
@@ -163,52 +148,53 @@ def get_precipitation_fields(
 
     # Read the radar composites
     importer = io.get_method(importer_name, "importer")
-
     reference_field = io.read_timeseries(fns, importer, **_importer_kwargs)
 
-    reference_field, _, ref_metadata = _xarray2legacy(reference_field)
-
-    if not return_raw:
-
-        if (num_prev_files == 0) and (num_next_files == 0):
-            # Remove time dimension
-            reference_field = np.squeeze(reference_field)
+    # Squeeze single fields
+    if (num_prev_files == 0) and (num_next_files == 0):
+        reference_field = reference_field.squeeze("t")
 
-        # Convert to mm/h
-        reference_field, ref_metadata = stp.utils.to_rainrate(
-            reference_field, ref_metadata
+    # Clip domain
+    if clip:
+        reference_field = reference_field.sel(
+            x=slice(clip[0], clip[1]),
+            y=slice(clip[2], clip[3]),
         )
 
-        # Clip domain
-        reference_field, ref_metadata = stp.utils.clip_domain(
-            reference_field, ref_metadata, clip
-        )
-
-        # Upscale data
-        reference_field, ref_metadata = aggregate_fields_space(
-            reference_field, ref_metadata, upscale
-        )
-
-        # Mask invalid values
-        reference_field = np.ma.masked_invalid(reference_field)
-
-        if log_transform:
-            # Log-transform the data [dBR]
-            reference_field, ref_metadata = stp.utils.dB_transform(
-                reference_field, ref_metadata, threshold=0.1, zerovalue=-15.0
-            )
-
-            # Set missing values with the fill value
-            np.ma.set_fill_value(reference_field, -15.0)
-            reference_field.data[reference_field.mask] = -15.0
-
-    if legacy or metadata:
-        if metadata:
+    # Fill nans
+    if filled:
+        units = reference_field.attrs.get("unit")
+        reference_field = reference_field.pysteps.to_rainrate()
+        reference_field = reference_field.fillna(0)
+        converter = utils.get_method(units)
+        reference_field = converter(reference_field)
+
+    # Coarsen data
+    if coarsen:
+        units = reference_field.attrs.get("unit")
+        reference_field = reference_field.pysteps.to_rainrate()
+        reference_field = reference_field.coarsen(
+            x=coarsen, y=coarsen, boundary="trim"
+        ).mean()
+        converter = utils.get_method(units)
+        reference_field = converter(reference_field)
+
+    # Convert
+    converter = utils.get_method(convert_to)
+    reference_field = converter(reference_field)
+
+    # Transform
+    transformer = utils.get_method(transform_to)
+    reference_field = transformer(reference_field)
+
+    if "legacy" in kwargs or "metadata" in kwargs:
+        reference_field, _, ref_metadata = _xarray2legacy(reference_field)
+        if kwargs.get("metadata", False):
             return reference_field, ref_metadata
         else:
             return reference_field
 
-    return _to_xarray(reference_field, ref_metadata)
+    return reference_field
 
 
 def smart_assert(actual_value, expected, tolerance=None):
diff --git a/pysteps/tests/test_cascade.py b/pysteps/tests/test_cascade.py
index 77951922c..f9e5b0d5b 100644
--- a/pysteps/tests/test_cascade.py
+++ b/pysteps/tests/test_cascade.py
@@ -11,46 +11,42 @@
 from pysteps.cascade.bandpass_filters import filter_gaussian
 from pysteps.cascade.bandpass_filters import filter_uniform
 from pysteps.cascade.decomposition import decomposition_fft, recompose_fft
-from pysteps.tests.helpers import smart_assert
+from pysteps.tests.helpers import smart_assert, get_precipitation_fields
 
-# TODO: Fix tests for xarray fields
-pytestmark = pytest.mark.skip("Needs migration to xarray")
+
+PRECIP = get_precipitation_fields(
+    source="bom",
+    convert_to="mm/h",
+    transform_to="db",
+    coarsen=2,
+    filled=True,
+)
 
 
 def test_decompose_recompose():
     """Tests cascade decomposition."""
 
-    pytest.importorskip("netCDF4")
-
-    root_path = pysteps.rcparams.data_sources["bom"]["root_path"]
-    rel_path = os.path.join("prcp-cscn", "2", "2018", "06", "16")
-    filename = os.path.join(root_path, rel_path, "2_20180616_120000.prcp-cscn.nc")
-    precip, _, metadata = pysteps.io.import_bom_rf3(filename)
-
-    # Convert to rain rate from mm
-    precip, metadata = pysteps.utils.to_rainrate(precip, metadata)
-
-    # Log-transform the data
-    precip, metadata = pysteps.utils.dB_transform(
-        precip, metadata, threshold=0.1, zerovalue=-15.0
-    )
-
-    # Set Nans as the fill value
-    precip[~np.isfinite(precip)] = metadata["zerovalue"]
-
     # Set number of cascade levels
     num_cascade_levels = 9
 
     # Construct the Gaussian bandpass filters
-    _filter = filter_gaussian(precip.shape, num_cascade_levels)
+    _filter = filter_gaussian(PRECIP.shape, num_cascade_levels)
 
     # Decompose precip
-    decomp = decomposition_fft(precip, _filter)
+    decomp = decomposition_fft(PRECIP, _filter)
+
+    assert all(key in decomp for key in ("cascade_levels", "means", "stds"))
+    assert decomp["cascade_levels"].shape[0] == num_cascade_levels
+    assert decomp["cascade_levels"].shape[1:] == PRECIP.shape
+    assert isinstance(decomp["means"], list)
+    assert isinstance(decomp["stds"], list)
+    assert len(decomp["means"]) == num_cascade_levels
+    assert len(decomp["stds"]) == num_cascade_levels
 
     # Recomposed precip from decomp
     recomposed = recompose_fft(decomp)
-    # Assert
-    assert_array_almost_equal(recomposed.squeeze(), precip)
+
+    assert_array_almost_equal(recomposed, PRECIP)
 
 
 test_metadata_filter = [
diff --git a/pysteps/tests/test_datasets.py b/pysteps/tests/test_datasets.py
index e524099e7..0ee8da70e 100644
--- a/pysteps/tests/test_datasets.py
+++ b/pysteps/tests/test_datasets.py
@@ -13,8 +13,6 @@
 )
 from pysteps.exceptions import DirectoryNotEmpty
 
-# TODO: Fix tests for xarray fields
-pytestmark = pytest.mark.skip("Needs migration to xarray")
 
 _datasets_opt_deps = dict(
     fmi=["pyproj"],
diff --git a/pysteps/tests/test_detcatscores.py b/pysteps/tests/test_detcatscores.py
index e1dc3de4b..ba14b23e6 100644
--- a/pysteps/tests/test_detcatscores.py
+++ b/pysteps/tests/test_detcatscores.py
@@ -6,8 +6,6 @@
 
 from pysteps.verification import det_cat_fct
 
-# TODO: Fix tests for xarray fields
-pytestmark = pytest.mark.skip("Needs migration to xarray")
 
 # CREATE A LARGE DATASET TO MATCH
 # EXAMPLES IN
diff --git a/pysteps/tests/test_downscaling_rainfarm.py b/pysteps/tests/test_downscaling_rainfarm.py
index f85a8cefe..c09c4009e 100644
--- a/pysteps/tests/test_downscaling_rainfarm.py
+++ b/pysteps/tests/test_downscaling_rainfarm.py
@@ -4,16 +4,13 @@
 
 from pysteps import downscaling
 from pysteps.tests.helpers import get_precipitation_fields
-from pysteps.utils import aggregate_fields_space, square_domain
-
 
 # load and preprocess input field
-precip, metadata = get_precipitation_fields(
-    num_prev_files=0, num_next_files=0, return_raw=False, metadata=True
+PRECIP = get_precipitation_fields(
+    source="bom",
+    filled=True,
+    convert_to="mm/h",
 )
-precip = precip.filled()
-precip, metadata = square_domain(precip, metadata, "crop")
-
 
 rainfarm_arg_names = ("alpha", "ds_factor", "threshold", "return_alpha")
 
@@ -23,18 +20,17 @@
 
 @pytest.mark.parametrize(rainfarm_arg_names, rainfarm_arg_values)
 def test_rainfarm_shape(alpha, ds_factor, threshold, return_alpha):
-    """Test that the output of rainfarm is consistent with the downscalnig factor."""
+    """Test that the output of rainfarm is consistent with the downscaling factor."""
 
-    window = metadata["xpixelsize"] * ds_factor
-    precip_lr, __ = aggregate_fields_space(precip, metadata, window)
+    precip_lr = PRECIP.coarsen(x=ds_factor, y=ds_factor).mean()
 
     rainfarm = downscaling.get_method("rainfarm")
 
     precip_hr = rainfarm(precip_lr, alpha, ds_factor, threshold, return_alpha)
 
-    assert precip_hr.ndim == precip.ndim
-    assert precip_hr.shape[0] == precip.shape[0]
-    assert precip_hr.shape[1] == precip.shape[1]
+    assert precip_hr.ndim == PRECIP.ndim
+    assert precip_hr.shape[0] == PRECIP.shape[0]
+    assert precip_hr.shape[1] == PRECIP.shape[1]
 
 
 rainfarm_arg_values = [(1.0, 2, 0, True), (None, 2, 0, True)]
@@ -44,9 +40,7 @@ def test_rainfarm_shape(alpha, ds_factor, threshold, return_alpha):
 def test_rainfarm_alpha(alpha, ds_factor, threshold, return_alpha):
     """Test that rainfarm computes and returns alpha."""
 
-    window = metadata["xpixelsize"] * ds_factor
-    precip_lr, __ = aggregate_fields_space(precip, metadata, window)
-
+    precip_lr = PRECIP.coarsen(x=ds_factor, y=ds_factor).mean()
     rainfarm = downscaling.get_method("rainfarm")
 
     precip_hr = rainfarm(precip_lr, alpha, ds_factor, threshold, return_alpha)
diff --git a/pysteps/tests/test_ensscores.py b/pysteps/tests/test_ensscores.py
index e1e7b89e8..f4bb3f582 100644
--- a/pysteps/tests/test_ensscores.py
+++ b/pysteps/tests/test_ensscores.py
@@ -7,13 +7,15 @@
 from pysteps.tests.helpers import get_precipitation_fields
 from pysteps.verification import ensscores
 
-precip = get_precipitation_fields(num_next_files=10, return_raw=True)
+PRECIP = get_precipitation_fields(
+    num_next_files=10,
+)
 np.random.seed(42)
 
 # rankhist
 test_data = [
-    (precip[:10], precip[-1], None, True, 11),
-    (precip[:10], precip[-1], None, False, 11),
+    (PRECIP.isel(t=slice(0, 10)), PRECIP.isel(t=-1), None, True, 11),
+    (PRECIP.isel(t=slice(0, 10)), PRECIP.isel(t=-1), None, False, 11),
 ]
 
 
@@ -28,14 +30,14 @@ def test_rankhist_size(X_f, X_o, X_min, normalize, expected):
 # ensemble_skill
 test_data = [
     (
-        precip[:10],
-        precip[-1],
+        PRECIP.isel(t=slice(0, 10)),
+        PRECIP.isel(t=-1),
         "RMSE",
         {"axis": None, "conditioning": "single"},
         0.26054151,
     ),
-    (precip[:10], precip[-1], "CSI", {"thr": 1.0, "axis": None}, 0.22017924),
-    (precip[:10], precip[-1], "FSS", {"thr": 1.0, "scale": 10}, 0.63239752),
+    (PRECIP.isel(t=slice(0, 10)), PRECIP.isel(t=-1), "CSI", {"thr": 1.0, "axis": None}, 0.22017924),
+    (PRECIP.isel(t=slice(0, 10)), PRECIP.isel(t=-1), "FSS", {"thr": 1.0, "scale": 10}, 0.63239752),
 ]
 
 
@@ -49,9 +51,9 @@ def test_ensemble_skill(X_f, X_o, metric, kwargs, expected):
 
 # ensemble_spread
 test_data = [
-    (precip, "RMSE", {"axis": None, "conditioning": "single"}, 0.22635757),
-    (precip, "CSI", {"thr": 1.0, "axis": None}, 0.25218158),
-    (precip, "FSS", {"thr": 1.0, "scale": 10}, 0.70235667),
+    (PRECIP, "RMSE", {"axis": None, "conditioning": "single"}, 0.22635757),
+    (PRECIP, "CSI", {"thr": 1.0, "axis": None}, 0.25218158),
+    (PRECIP, "FSS", {"thr": 1.0, "scale": 10}, 0.70235667),
 ]
 
 
diff --git a/pysteps/tests/test_exporters.py b/pysteps/tests/test_exporters.py
index 9515828b7..ca27c9533 100644
--- a/pysteps/tests/test_exporters.py
+++ b/pysteps/tests/test_exporters.py
@@ -5,7 +5,7 @@
 from datetime import datetime
 
 import numpy as np
-import pytest
+import pandas as pd
 from numpy.testing import assert_array_almost_equal
 
 from pysteps.io import import_netcdf_pysteps
@@ -15,11 +15,6 @@
 from pysteps.io.exporters import initialize_forecast_exporter_netcdf
 from pysteps.tests.helpers import get_precipitation_fields, get_invalid_mask
 
-# TODO: Fix tests for xarray fields
-pytestmark = pytest.mark.skip(
-    "Needs migration to xarray. Problem with cartesian_unit attribute."
-)
-
 
 def test_get_geotiff_filename():
     """Test the geotif name generator."""
@@ -43,21 +38,15 @@ def test_io_export_netcdf_one_member_one_time_step():
     """Test the export netcdf.
     Also, test that the exported file can be read by the importer."""
 
-    pytest.importorskip("pyproj")
-
-    precip, metadata = get_precipitation_fields(
-        return_raw=True, metadata=True, source="fmi"
-    )
-    precip = precip.squeeze()
-
+    precip = get_precipitation_fields(source="fmi")
     invalid_mask = get_invalid_mask(precip)
 
     # save it back to disk
     with tempfile.TemporaryDirectory() as outpath:
         outfnprefix = "test_netcdf_out"
         file_path = os.path.join(outpath, outfnprefix + ".nc")
-        startdate = metadata["timestamps"][0]
-        timestep = metadata["accutime"]
+        startdate = pd.to_datetime(precip.t.values).to_pydatetime()
+        timestep = precip.attrs.get("accutime")
         n_timesteps = 1
         shape = precip.shape
         exporter = initialize_forecast_exporter_netcdf(
@@ -67,10 +56,10 @@ def test_io_export_netcdf_one_member_one_time_step():
             timestep,
             n_timesteps,
             shape,
-            metadata,
+            precip.attrs,
             n_ens_members=1,
         )
-        export_forecast_dataset(precip[np.newaxis, :], exporter)
+        export_forecast_dataset(precip.expand_dims("t"), exporter)
         close_forecast_files(exporter)
 
         # assert if netcdf file was saved and file size is not zero
@@ -79,15 +68,15 @@ def test_io_export_netcdf_one_member_one_time_step():
         # Test that the file can be read by the nowcast_importer
         output_file_path = os.path.join(outpath, f"{outfnprefix}.nc")
 
-        precip_new, _ = import_netcdf_pysteps(output_file_path)
+        precip_new = import_netcdf_pysteps(output_file_path)
 
         assert_array_almost_equal(precip, precip_new.data)
         assert precip_new.dtype == "single"
 
-        precip_new, _ = import_netcdf_pysteps(output_file_path, dtype="double")
+        precip_new = import_netcdf_pysteps(output_file_path, dtype="double")
         assert_array_almost_equal(precip, precip_new.data)
         assert precip_new.dtype == "double"
 
-        precip_new, _ = import_netcdf_pysteps(output_file_path, fillna=-1000)
+        precip_new = import_netcdf_pysteps(output_file_path, fillna=-1000)
         new_invalid_mask = precip_new == -1000
         assert (new_invalid_mask == invalid_mask).all()
diff --git a/pysteps/tests/test_feature.py b/pysteps/tests/test_feature.py
index e3873efd9..24a366154 100644
--- a/pysteps/tests/test_feature.py
+++ b/pysteps/tests/test_feature.py
@@ -9,12 +9,12 @@
 
 @pytest.mark.parametrize(arg_names, arg_values)
 def test_feature(method, max_num_features):
-    input_field, _ = get_precipitation_fields(0, 0, True, True, None, "mch")
+    input_field = get_precipitation_fields(source="mch", convert_to="reflectivity")
 
     detector = feature.get_method(method)
 
     kwargs = {"max_num_features": max_num_features}
-    output = detector(input_field.squeeze(), **kwargs)
+    output = detector(input_field, **kwargs)
 
     assert isinstance(output, np.ndarray)
     assert output.ndim == 2
diff --git a/pysteps/tests/test_feature_tstorm.py b/pysteps/tests/test_feature_tstorm.py
index df7c3dc5a..3b58725d4 100644
--- a/pysteps/tests/test_feature_tstorm.py
+++ b/pysteps/tests/test_feature_tstorm.py
@@ -2,7 +2,6 @@
 import pytest
 
 from pysteps.feature.tstorm import detection
-from pysteps.utils import to_reflectivity
 from pysteps.tests.helpers import get_precipitation_fields
 
 try:
@@ -29,15 +28,13 @@ def test_feature_tstorm_detection(source, output_feat, dry_input, max_num_featur
     pytest.importorskip("pandas")
 
     if not dry_input:
-        input, metadata = get_precipitation_fields(0, 0, True, True, None, source)
-        input = input.squeeze()
-        input, __ = to_reflectivity(input, metadata)
+        input_image = get_precipitation_fields(source=source, convert_to="reflectivity")
     else:
-        input = np.zeros((50, 50))
+        input_image = np.zeros((50, 50))
 
     time = "000"
     output = detection(
-        input, time=time, output_feat=output_feat, max_num_features=max_num_features
+        input_image, time=time, output_feat=output_feat, max_num_features=max_num_features
     )
 
     if output_feat:
@@ -67,7 +64,7 @@ def test_feature_tstorm_detection(source, output_feat, dry_input, max_num_featur
         ]
         assert (output[0].time == time).all()
         assert output[1].ndim == 2
-        assert output[1].shape == input.shape
+        assert output[1].shape == input_image.shape
         if not dry_input:
             assert output[0].shape[0] > 0
             assert sorted(list(output[0].ID)) == sorted(list(np.unique(output[1]))[1:])
diff --git a/pysteps/tests/test_importer_decorator.py b/pysteps/tests/test_importer_decorator.py
index 85d09b946..ccd8d0292 100644
--- a/pysteps/tests/test_importer_decorator.py
+++ b/pysteps/tests/test_importer_decorator.py
@@ -6,40 +6,36 @@
 
 from pysteps.tests.helpers import get_precipitation_fields
 
-default_dtypes = dict(
-    fmi="double",
-    knmi="double",
-    mch="double",
-    opera="double",
-    saf="double",
-    mrms="single",
-)
-
-
-@pytest.mark.parametrize("source, default_dtype", default_dtypes.items())
-def test_postprocess_import_decorator(source, default_dtype):
+DEFAULT_DTYPE = "float32"
+
+def test_postprocess_import_decorator():
     """Test the postprocessing decorator for the importers."""
-    import_data = partial(get_precipitation_fields, return_raw=True, source=source)
+    import_data = partial(get_precipitation_fields, source="mch")
 
     precip = import_data()
     invalid_mask = ~np.isfinite(precip)
 
-    assert precip.dtype == default_dtype
+    assert precip.dtype == DEFAULT_DTYPE
 
-    if default_dtype == "single":
-        dtype = "double"
+    if DEFAULT_DTYPE == "float32":
+        dtype = "float64"
     else:
-        dtype = "single"
+        dtype = "float32"
 
-    precip = import_data(dtype=dtype)
+    precip = import_data(importer_kwargs=dict(dtype=dtype))
 
     assert precip.dtype == dtype
 
+    # Test that equivalent dtypes are handled correctly
+    for dtype1, dtype2 in zip(["float64", "float32",], ["double", "single"]):
+        precip = import_data(importer_kwargs=dict(dtype=dtype1))
+        assert precip.dtype == dtype2
+
     # Test that invalid types are handled correctly
     for dtype in ["int", "int64"]:
         with pytest.raises(ValueError):
-            _ = import_data(dtype=dtype)
+            import_data(importer_kwargs=dict(dtype=dtype))
 
-    precip = import_data(fillna=-1000)
+    precip = import_data(importer_kwargs=dict(fillna=-1000))
     new_invalid_mask = precip == -1000
     assert (new_invalid_mask == invalid_mask).all()
diff --git a/pysteps/tests/test_interfaces.py b/pysteps/tests/test_interfaces.py
index 49b6bb69b..e92b67ba6 100644
--- a/pysteps/tests/test_interfaces.py
+++ b/pysteps/tests/test_interfaces.py
@@ -248,7 +248,6 @@ def test_utils_interface():
     from pysteps.utils import arrays
     from pysteps.utils import cleansing
     from pysteps.utils import conversion
-    from pysteps.utils import dimension
     from pysteps.utils import images
     from pysteps.utils import interpolate
     from pysteps.utils import spectral
@@ -266,20 +265,16 @@ def test_utils_interface():
         ("raindepth", conversion.to_raindepth),
         ("dbz", conversion.to_reflectivity),
         ("reflectivity", conversion.to_reflectivity),
-        ("accumulate", dimension.aggregate_fields_time),
-        ("clip", dimension.clip_domain),
-        ("square", dimension.square_domain),
-        ("upscale", dimension.aggregate_fields_space),
         ("morph_opening", images.morph_opening),
         ("rbfinterp2d", interpolate.rbfinterp2d),
         ("rapsd", spectral.rapsd),
         ("rm_rdisc", spectral.remove_rain_norain_discontinuity),
         ("boxcox", transformation.boxcox_transform),
         ("box-cox", transformation.boxcox_transform),
-        ("db", transformation.dB_transform),
-        ("decibel", transformation.dB_transform),
+        ("db", transformation.db_transform),
+        ("decibel", transformation.db_transform),
         ("log", transformation.boxcox_transform),
-        ("nqt", transformation.NQ_transform),
+        ("nqt", transformation.nq_transform),
         ("sqrt", transformation.sqrt_transform),
     ]
 
diff --git a/pysteps/tests/test_io_bom_nwp.py b/pysteps/tests/test_io_bom_nwp.py
index d7deae63f..83368d273 100644
--- a/pysteps/tests/test_io_bom_nwp.py
+++ b/pysteps/tests/test_io_bom_nwp.py
@@ -11,14 +11,17 @@
 
 pytest.importorskip("netCDF4")
 
-root_path = pysteps.rcparams.data_sources["bom_nwp"]["root_path"]
-rel_path = os.path.join("2020", "10", "31")
-filename = os.path.join(root_path, rel_path, "20201031_0000_regrid_short.nc")
-data_array_xr = pysteps.io.import_bom_nwp_xr(filename)
+try:
+    root_path = pysteps.rcparams.data_sources["bom_nwp"]["root_path"]
+except KeyError:
+    pytestmark = pytest.mark.skip("all tests still WIP")
+else:
+    rel_path = os.path.join("2020", "10", "31")
+    filename = os.path.join(root_path, rel_path, "20201031_0000_regrid_short.nc")
+    data_array_xr = pysteps.io.import_bom_nwp_xr(filename)
 
 expected_proj = "+proj=aea  +lon_0=153.240 +lat_0=-27.718 +lat_1=-26.200 +lat_2=-29.300"
 
-
 def test_io_import_bom_nwp_xarray():
     """Test the importer Bom NWP."""
     assert isinstance(data_array_xr, xr.DataArray)
diff --git a/pysteps/tests/test_io_knmi_nwp.py b/pysteps/tests/test_io_knmi_nwp.py
index 95c332dcd..4d6a02a91 100644
--- a/pysteps/tests/test_io_knmi_nwp.py
+++ b/pysteps/tests/test_io_knmi_nwp.py
@@ -11,10 +11,14 @@
 
 pytest.importorskip("netCDF4")
 
-root_path = pysteps.rcparams.data_sources["knmi_nwp"]["root_path"]
-rel_path = os.path.join("2018", "09", "05")
-filename = os.path.join(root_path, rel_path, "20180905_0600_Pforecast_Harmonie.nc")
-data_array_xr = pysteps.io.import_knmi_nwp_xr(filename)
+try:
+    root_path = pysteps.rcparams.data_sources["knmi_nwp"]["root_path"]
+except KeyError:
+    pytestmark = pytest.mark.skip("all tests still WIP")
+else:
+    rel_path = os.path.join("2018", "09", "05")
+    filename = os.path.join(root_path, rel_path, "20180905_0600_Pforecast_Harmonie.nc")
+    data_array_xr = pysteps.io.import_knmi_nwp_xr(filename)
 
 expected_proj = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
 
diff --git a/pysteps/tests/test_io_mrms_grib.py b/pysteps/tests/test_io_mrms_grib.py
index afebb7b37..ea3740ab9 100644
--- a/pysteps/tests/test_io_mrms_grib.py
+++ b/pysteps/tests/test_io_mrms_grib.py
@@ -33,10 +33,7 @@ def test_io_import_mrms_grib():
         "transform": None,
         "zerovalue": 0,
         "threshold": 0.1,
-        "x1": -129.99999999999997,
-        "x2": -60.00000199999991,
-        "y1": 20.000001,
-        "y2": 55.00000000000001,
+        "bounding_box": (-129.99999999999997, -60.00000199999991, 20.000001, 55.00000000000001),
         "units": "degrees",
     }
     metadata = data_Array.attrs
@@ -48,8 +45,9 @@ def test_io_import_mrms_grib():
         else:
             assert metadata[key] == expected_metadata[key]
 
-    x = np.arange(metadata["x1"], metadata["x2"], metadata["xpixelsize"])
-    y = np.arange(metadata["y1"], metadata["y2"], metadata["ypixelsize"])
+    x1, x2, y1, y2 = metadata["bounding_box"]
+    x = np.arange(x1, x2, metadata["xpixelsize"])
+    y = np.arange(y1, y2, metadata["ypixelsize"])
     precip_full = data_Array.values
     assert y.size == precip_full.shape[0]
     assert x.size == precip_full.shape[1]
@@ -86,30 +84,28 @@ def test_io_import_mrms_grib():
     assert precip_single.dtype == "single"
     del precip_single
 
-    precip_donwscaled = pysteps.io.import_mrms_grib(
-        filename, dtype="single", fillna=0, window_size=2
-    )
-    assert precip_donwscaled.shape == (3500 / 2, 7000 / 2)
+    # TODO: add downscaling feauture to importer decorator
+    # precip_donwscaled = pysteps.io.import_mrms_grib(
+    #     filename, dtype="single", fillna=0, window_size=2
+    # )
+    # assert precip_donwscaled.shape == (3500 / 2, 7000 / 2)
 
-    precip_donwscaled = pysteps.io.import_mrms_grib(
-        filename, dtype="single", fillna=0, window_size=3
-    )
-    print(metadata)
-    expected_metadata.update(
-        xpixelsize=0.03,
-        ypixelsize=0.03,
-        x1=-130.00000000028575,
-        x2=-60.01000199942843,
-        y1=20.02000099914261,
-        y2=55.000000000285794,
-    )
+    # precip_donwscaled = pysteps.io.import_mrms_grib(
+    #     filename, dtype="single", fillna=0, window_size=3
+    # )
+    # print(metadata)
+    # expected_metadata.update(
+    #     xpixelsize=0.03,
+    #     ypixelsize=0.03,
+    #     bounding_box=(-130.00000000028575, 60.01000199942843, 20.02000099914261, 55.000000000285794),
+    # )
 
     # Remove the threshold keyword from the test when the window_size>1 is used.
     # The threshold is computed automatically from the minimum valid precipitation values.
     # The minimum value is affected by the the block_averaging (window_size=3 keyword)
     # of the precipitation fields. Hence, the "threshold" value will depend on the
     # precipitation pattern (i.e. the file being read).
-    expected_metadata.pop("threshold", None)
+    # expected_metadata.pop("threshold", None)
 
     # expected_metadata.update(
     #     xpixelsize=0.03,
diff --git a/pysteps/tests/test_io_rmi_nwp.py b/pysteps/tests/test_io_rmi_nwp.py
index e4ab07bed..9a5ca845f 100644
--- a/pysteps/tests/test_io_rmi_nwp.py
+++ b/pysteps/tests/test_io_rmi_nwp.py
@@ -14,10 +14,14 @@
 
 pytest.importorskip("netCDF4")
 
-root_path = pysteps.rcparams.data_sources["rmi_nwp"]["root_path"]
-rel_path = os.path.join("2021", "07", "04")
-filename = os.path.join(root_path, rel_path, "ao13_2021070412_native_5min.nc")
-data_array_xr = pysteps.io.import_rmi_nwp_xr(filename)
+try:
+    root_path = pysteps.rcparams.data_sources["rmi_nwp"]["root_path"]
+except KeyError:
+    pytestmark = pytest.mark.skip("all tests still WIP")
+else:
+    rel_path = os.path.join("2021", "07", "04")
+    filename = os.path.join(root_path, rel_path, "ao13_2021070412_native_5min.nc")
+    data_array_xr = pysteps.io.import_rmi_nwp_xr(filename)
 
 expected_proj = "+proj=lcc +lon_0=4.55 +lat_1=50.8 +lat_2=50.8 +a=6371229 +es=0 +lat_0=50.8 +x_0=365950 +y_0=-365950"
 
diff --git a/pysteps/tests/test_motion.py b/pysteps/tests/test_motion.py
index 1cae9ec98..ff66aeb65 100644
--- a/pysteps/tests/test_motion.py
+++ b/pysteps/tests/test_motion.py
@@ -225,8 +225,8 @@ def test_optflow_method_convergence(
 
         retrieved_motion = oflow_method(precip_obs, verbose=False)
 
-    precip_data, _ = stp.utils.dB_transform(precip_obs.max(axis=0), inverse=True)
-    precip_data.data[precip_data.mask] = 0
+    precip_data = precip_obs.pysteps.db_transform(precip_obs.max("t"), inverse=True)
+    precip_data = precip_data.fillna(0)
 
     precip_mask = (uniform_filter(precip_data, size=20) > 0.1) & ~precip_obs.mask.any(
         axis=0
diff --git a/pysteps/tests/test_motion_lk.py b/pysteps/tests/test_motion_lk.py
index 871dcd98b..6e1452bed 100644
--- a/pysteps/tests/test_motion_lk.py
+++ b/pysteps/tests/test_motion_lk.py
@@ -61,14 +61,12 @@ def test_lk(
         pytest.importorskip("pandas")
 
     # inputs
-    precip, metadata = get_precipitation_fields(
+    precip = get_precipitation_fields(
         num_prev_files=2,
         num_next_files=0,
-        return_raw=False,
-        metadata=True,
-        upscale=2000,
+        coarsen=2,
+        filled=True,
     )
-    precip = precip.filled()
 
     # Retrieve motion field
     oflow_method = motion.get_method("LK")
diff --git a/pysteps/tests/test_noise_fftgenerators.py b/pysteps/tests/test_noise_fftgenerators.py
index f0da5b759..ab636fc3e 100644
--- a/pysteps/tests/test_noise_fftgenerators.py
+++ b/pysteps/tests/test_noise_fftgenerators.py
@@ -8,10 +8,9 @@
     num_prev_files=0,
     num_next_files=0,
     return_raw=False,
-    metadata=False,
     upscale=2000,
 )
-PRECIP = PRECIP.filled()
+PRECIP = PRECIP.fillna(PRECIP.min())
 
 
 def test_noise_param_2d_fft_filter():
diff --git a/pysteps/tests/test_plt_cartopy.py b/pysteps/tests/test_plt_cartopy.py
index e20aaee8a..2aeb9ef71 100644
--- a/pysteps/tests/test_plt_cartopy.py
+++ b/pysteps/tests/test_plt_cartopy.py
@@ -3,7 +3,6 @@
 import pytest
 
 from pysteps.visualization import plot_precip_field
-from pysteps.utils import to_rainrate
 from pysteps.tests.helpers import get_precipitation_fields
 import matplotlib.pyplot as plt
 
@@ -27,7 +26,7 @@
 def test_visualization_plot_precip_field(source, map_kwargs, pass_geodata):
     field, metadata = get_precipitation_fields(0, 0, True, True, None, source)
     field = field.squeeze()
-    field, __ = to_rainrate(field, metadata)
+    field = field.pysteps.to_rainrate()
 
     if not pass_geodata:
         metadata = None
diff --git a/pysteps/tests/test_plt_precipfields.py b/pysteps/tests/test_plt_precipfields.py
index 3ac8dffae..5e53af478 100644
--- a/pysteps/tests/test_plt_precipfields.py
+++ b/pysteps/tests/test_plt_precipfields.py
@@ -43,18 +43,18 @@ def test_visualization_plot_precip_field(
 
         field, metadata = get_precipitation_fields(0, 0, True, True, None, source)
         field = field.squeeze()
-        field, metadata = conversion.to_rainrate(field, metadata)
+        field = field.pysteps.to_rainrate()
 
     elif plot_type == "depth":
 
         field, metadata = get_precipitation_fields(0, 0, True, True, None, source)
         field = field.squeeze()
-        field, metadata = conversion.to_raindepth(field, metadata)
+        field = field.pysteps.to_raindepth()
 
     elif plot_type == "prob":
 
         field, metadata = get_precipitation_fields(0, 10, True, True, None, source)
-        field, metadata = conversion.to_rainrate(field, metadata)
+        field = field.pysteps.to_rainrate()
         field = ensemblestats.excprob(field, probthr)
 
     field_orig = field.copy()
diff --git a/pysteps/tests/test_tracking_tdating.py b/pysteps/tests/test_tracking_tdating.py
index f8efe2d17..e62db1991 100644
--- a/pysteps/tests/test_tracking_tdating.py
+++ b/pysteps/tests/test_tracking_tdating.py
@@ -4,7 +4,6 @@
 import pytest
 
 from pysteps.tracking.tdating import dating
-from pysteps.utils import to_reflectivity
 from pysteps.tests.helpers import get_precipitation_fields
 
 arg_names = ("source", "dry_input")
@@ -23,7 +22,7 @@ def test_tracking_tdating_dating(source, dry_input):
 
     if not dry_input:
         input, metadata = get_precipitation_fields(0, 2, True, True, 4000, source)
-        input, __ = to_reflectivity(input, metadata)
+        input = input.pysteps.to_reflectivity()
     else:
         input = np.zeros((3, 50, 50))
         metadata = {"timestamps": ["00", "01", "02"]}
diff --git a/pysteps/tests/test_utils_conversion.py b/pysteps/tests/test_utils_conversion.py
index 169cdb50e..cd7251f34 100644
--- a/pysteps/tests/test_utils_conversion.py
+++ b/pysteps/tests/test_utils_conversion.py
@@ -1,336 +1,306 @@
 # -*- coding: utf-8 -*-
 
-import numpy as np
 import pytest
+import xarray as xr
 from numpy.testing import assert_array_almost_equal
 
-from pysteps.utils import conversion
-
 # to_rainrate
 test_data = [
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
-            "transform": None,
             "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([1]),
+        1.0,
     ),
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": None,
             "unit": "mm",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([12]),
+        12.0,
     ),
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": "dB",
+            "offset": 0.0,
             "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([1.25892541]),
+        1.2589,
     ),
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": "dB",
+            "offset": 0.0,
             "unit": "mm",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([15.10710494]),
+        15.1071,
     ),
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": "dB",
+            "offset": 0.0,
             "unit": "dBZ",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([0.04210719]),
+        0.0421,
     ),
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": "log",
+            "offset": 0.0,
             "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([2.71828183]),
+        2.7183,
     ),
     (
-        np.array([1.0]),
+        1.0,
         {
             "accutime": 5,
             "transform": "log",
+            "offset": 0.0,
             "unit": "mm",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([32.61938194]),
+        32.6194,
     ),
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": "sqrt",
             "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([1]),
+        1.0,
     ),
     (
-        np.array([1.0]),
+        1.0,
         {
             "accutime": 5,
             "transform": "sqrt",
             "unit": "mm",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([12.0]),
+        12.0,
     ),
 ]
 
 
-@pytest.mark.parametrize("R, metadata, expected", test_data)
-def test_to_rainrate(R, metadata, expected):
+@pytest.mark.parametrize("data, attrs, expected", test_data)
+def test_to_rainrate(data, attrs, expected):
     """Test the to_rainrate."""
-    assert_array_almost_equal(conversion.to_rainrate(R, metadata)[0], expected)
+    data_array = xr.DataArray([data], dims="x", attrs=attrs)
+    output = data_array.pysteps.to_rainrate()
+    assert output.attrs.get("unit") == "mm/h"
+    assert_array_almost_equal(output.values, [expected], decimal=4)
 
 
 # to_raindepth
 test_data = [
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": None,
             "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([0.08333333]),
+        0.0833,
     ),
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": None,
             "unit": "mm",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([1]),
+        1.0,
     ),
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": "dB",
+            "offset": 0.0,
             "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([0.10491045]),
+        0.1049,
     ),
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": "dB",
+            "offset": 0.0,
             "unit": "mm",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([1.25892541]),
+        1.2589,
     ),
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": "dB",
+            "offset": 0.0,
             "unit": "dBZ",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([0.00350893]),
+        0.00353,
     ),
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": "log",
+            "offset": 0.0,
             "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([0.22652349]),
+        0.2265,
     ),
     (
-        np.array([1.0]),
+        1.0,
         {
             "accutime": 5,
             "transform": "log",
+            "offset": 0.0,
             "unit": "mm",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([2.71828183]),
+        2.7183,
     ),
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": "sqrt",
             "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([0.08333333]),
+        0.0833,
     ),
     (
-        np.array([1.0]),
+        1.0,
         {
             "accutime": 5,
             "transform": "sqrt",
             "unit": "mm",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([1.0]),
+        1.0,
     ),
 ]
 
 
-@pytest.mark.parametrize("R, metadata, expected", test_data)
-def test_to_raindepth(R, metadata, expected):
+@pytest.mark.parametrize("data, attrs, expected", test_data)
+def test_to_raindepth(data, attrs, expected):
     """Test the to_raindepth."""
-    assert_array_almost_equal(conversion.to_raindepth(R, metadata)[0], expected)
+    data_array = xr.DataArray([data], dims="x", attrs=attrs)
+    output = data_array.pysteps.to_raindepth()
+    assert output.attrs.get("unit") == "mm"
+    assert_array_almost_equal(output.values, [expected], decimal=4)
 
 
 # to_reflectivity
 test_data = [
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": None,
             "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([23.01029996]),
+        23.0103,
     ),
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": None,
             "unit": "mm",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([40.27719989]),
+        40.2772,
     ),
     (
-        np.array([1]),
+        1,
         {
             "accutime": 5,
             "transform": "dB",
+            "offset": 0.0,
             "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([24.61029996]),
+        24.6103,
     ),
     (
-        np.array([1]),
+        1,
         {
             "accutime": 5,
             "transform": "dB",
+            "offset": 0.0,
             "unit": "mm",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([41.87719989]),
+        41.8772,
     ),
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": "dB",
+            "offset": 0.0,
             "unit": "dBZ",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([1]),
+        1.0,
     ),
     (
-        np.array([1]),
+        1.0,
         {
             "accutime": 5,
             "transform": "log",
+            "offset": 0.0,
             "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([29.95901167]),
+        29.9590,
     ),
     (
-        np.array([1.0]),
+        1.0,
         {
             "accutime": 5,
             "transform": "log",
+            "offset": 0.0,
             "unit": "mm",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([47.2259116]),
+        47.2259,
     ),
     (
-        np.array([1]),
+        1,
         {
             "accutime": 5,
             "transform": "sqrt",
             "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([23.01029996]),
+        23.0103,
     ),
     (
-        np.array([1.0]),
+        1.0,
         {
             "accutime": 5,
             "transform": "sqrt",
             "unit": "mm",
-            "threshold": 0,
-            "zerovalue": 0,
         },
-        np.array([40.27719989]),
+        40.2772,
     ),
 ]
 
 
-@pytest.mark.parametrize("R, metadata, expected", test_data)
-def test_to_reflectivity(R, metadata, expected):
+@pytest.mark.parametrize("data, attrs, expected", test_data)
+def test_to_reflectivity(data, attrs, expected):
     """Test the to_reflectivity."""
-    assert_array_almost_equal(conversion.to_reflectivity(R, metadata)[0], expected)
+    data_array = xr.DataArray([data], dims="x", attrs=attrs)
+    output = data_array.pysteps.to_reflectivity(offset=0.0, to_decibels=True)
+    assert output.attrs.get("unit") == "dBZ"
+    assert_array_almost_equal(output.values, [expected], decimal=4)
+    output = data_array.pysteps.to_reflectivity(offset=0.0, to_decibels=False)
+    assert output.attrs.get("unit") == "Z"
+    assert_array_almost_equal(output.values, [10.0 ** (expected / 10.0)], decimal=0)
diff --git a/pysteps/tests/test_utils_dimension.py b/pysteps/tests/test_utils_dimension.py
deleted file mode 100644
index e35940d0b..000000000
--- a/pysteps/tests/test_utils_dimension.py
+++ /dev/null
@@ -1,286 +0,0 @@
-# -*- coding: utf-8 -*-
-
-import datetime as dt
-
-import numpy as np
-import pytest
-from numpy.testing import assert_array_equal
-from pytest import raises
-
-from pysteps.utils import dimension
-
-test_data_not_trim = (
-    # "data, window_size, axis, method, expected"
-    (np.arange(6), 2, 0, "mean", np.array([0.5, 2.5, 4.5])),
-    (
-        np.arange(4 * 6).reshape(4, 6),
-        (2, 3),
-        (0, 1),
-        "sum",
-        np.array([[24, 42], [96, 114]]),
-    ),
-    (
-        np.arange(4 * 6).reshape(4, 6),
-        (2, 2),
-        (0, 1),
-        "sum",
-        np.array([[14, 22, 30], [62, 70, 78]]),
-    ),
-    (
-        np.arange(4 * 6).reshape(4, 6),
-        2,
-        (0, 1),
-        "sum",
-        np.array([[14, 22, 30], [62, 70, 78]]),
-    ),
-    (
-        np.arange(4 * 6).reshape(4, 6),
-        (2, 3),
-        (0, 1),
-        "mean",
-        np.array([[4.0, 7.0], [16.0, 19.0]]),
-    ),
-    (
-        np.arange(4 * 6).reshape(4, 6),
-        (2, 2),
-        (0, 1),
-        "mean",
-        np.array([[3.5, 5.5, 7.5], [15.5, 17.5, 19.5]]),
-    ),
-    (
-        np.arange(4 * 6).reshape(4, 6),
-        2,
-        (0, 1),
-        "mean",
-        np.array([[3.5, 5.5, 7.5], [15.5, 17.5, 19.5]]),
-    ),
-)
-
-
-@pytest.mark.parametrize(
-    "data, window_size, axis, method, expected", test_data_not_trim
-)
-def test_aggregate_fields(data, window_size, axis, method, expected):
-    """
-    Test the aggregate_fields function.
-    The windows size must divide exactly the data dimensions.
-    Internally, additional test are generated for situations where the
-    windows size does not divide the data dimensions.
-    The length of each dimension should be larger than 2.
-    """
-
-    assert_array_equal(
-        dimension.aggregate_fields(data, window_size, axis=axis, method=method),
-        expected,
-    )
-
-    # Test the trimming capabilities.
-    data = np.pad(data, (0, 1))
-    assert_array_equal(
-        dimension.aggregate_fields(
-            data, window_size, axis=axis, method=method, trim=True
-        ),
-        expected,
-    )
-
-    with raises(ValueError):
-        dimension.aggregate_fields(data, window_size, axis=axis, method=method)
-
-
-def test_aggregate_fields_errors():
-    """
-    Test that the errors are correctly captured in the aggregate_fields
-    function.
-    """
-    data = np.arange(4 * 6).reshape(4, 6)
-
-    with raises(ValueError):
-        dimension.aggregate_fields(data, -1, axis=0)
-    with raises(ValueError):
-        dimension.aggregate_fields(data, 0, axis=0)
-    with raises(ValueError):
-        dimension.aggregate_fields(data, 1, method="invalid")
-
-    with raises(TypeError):
-        dimension.aggregate_fields(data, (1, 1), axis=0)
-
-
-# aggregate_fields_time
-timestamps = [dt.datetime.now() + dt.timedelta(minutes=t) for t in range(10)]
-test_data = [
-    (
-        np.ones((10, 1, 1)),
-        {"unit": "mm/h", "timestamps": timestamps},
-        2,
-        False,
-        np.ones((5, 1, 1)),
-    ),
-    (
-        np.ones((10, 1, 1)),
-        {"unit": "mm", "timestamps": timestamps},
-        2,
-        False,
-        2 * np.ones((5, 1, 1)),
-    ),
-]
-
-
-@pytest.mark.parametrize(
-    "R, metadata, time_window_min, ignore_nan, expected", test_data
-)
-def test_aggregate_fields_time(R, metadata, time_window_min, ignore_nan, expected):
-    """Test the aggregate_fields_time."""
-    assert_array_equal(
-        dimension.aggregate_fields_time(R, metadata, time_window_min, ignore_nan)[0],
-        expected,
-    )
-
-
-# aggregate_fields_space
-test_data = [
-    (
-        np.ones((1, 10, 10)),
-        {"unit": "mm/h", "xpixelsize": 1, "ypixelsize": 1},
-        2,
-        False,
-        np.ones((1, 5, 5)),
-    ),
-    (
-        np.ones((1, 10, 10)),
-        {"unit": "mm", "xpixelsize": 1, "ypixelsize": 1},
-        2,
-        False,
-        4 * np.ones((1, 5, 5)),
-    ),
-]
-
-
-@pytest.mark.parametrize("R, metadata, space_window, ignore_nan, expected", test_data)
-def test_aggregate_fields_space(R, metadata, space_window, ignore_nan, expected):
-    """Test the aggregate_fields_space."""
-    assert_array_equal(
-        dimension.aggregate_fields_space(R, metadata, space_window, ignore_nan)[0],
-        expected,
-    )
-
-
-# clip_domain
-R = np.zeros((4, 4))
-R[:2, :] = 1
-test_data = [
-    (
-        R,
-        {
-            "x1": 0,
-            "x2": 4,
-            "y1": 0,
-            "y2": 4,
-            "xpixelsize": 1,
-            "ypixelsize": 1,
-            "zerovalue": 0,
-            "yorigin": "upper",
-        },
-        None,
-        R,
-    ),
-    (
-        R,
-        {
-            "x1": 0,
-            "x2": 4,
-            "y1": 0,
-            "y2": 4,
-            "xpixelsize": 1,
-            "ypixelsize": 1,
-            "zerovalue": 0,
-            "yorigin": "lower",
-        },
-        (2, 4, 2, 4),
-        np.zeros((2, 2)),
-    ),
-    (
-        R,
-        {
-            "x1": 0,
-            "x2": 4,
-            "y1": 0,
-            "y2": 4,
-            "xpixelsize": 1,
-            "ypixelsize": 1,
-            "zerovalue": 0,
-            "yorigin": "upper",
-        },
-        (2, 4, 2, 4),
-        np.ones((2, 2)),
-    ),
-]
-
-
-@pytest.mark.parametrize("R, metadata, extent, expected", test_data)
-def test_clip_domain(R, metadata, extent, expected):
-    """Test the clip_domain."""
-    assert_array_equal(dimension.clip_domain(R, metadata, extent)[0], expected)
-
-
-# square_domain
-R = np.zeros((4, 2))
-test_data = [
-    # square by padding
-    (
-        R,
-        {"x1": 0, "x2": 2, "y1": 0, "y2": 4, "xpixelsize": 1, "ypixelsize": 1},
-        "pad",
-        False,
-        np.zeros((4, 4)),
-    ),
-    # square by cropping
-    (
-        R,
-        {"x1": 0, "x2": 2, "y1": 0, "y2": 4, "xpixelsize": 1, "ypixelsize": 1},
-        "crop",
-        False,
-        np.zeros((2, 2)),
-    ),
-    # inverse square by padding
-    (
-        np.zeros((4, 4)),
-        {
-            "x1": -1,
-            "x2": 3,
-            "y1": 0,
-            "y2": 4,
-            "xpixelsize": 1,
-            "ypixelsize": 1,
-            "orig_domain": (4, 2),
-            "square_method": "pad",
-        },
-        "pad",
-        True,
-        R,
-    ),
-    # inverse square by cropping
-    (
-        np.zeros((2, 2)),
-        {
-            "x1": 0,
-            "x2": 2,
-            "y1": 1,
-            "y2": 3,
-            "xpixelsize": 1,
-            "ypixelsize": 1,
-            "orig_domain": (4, 2),
-            "square_method": "crop",
-        },
-        "crop",
-        True,
-        R,
-    ),
-]
-
-
-@pytest.mark.parametrize("R, metadata, method, inverse, expected", test_data)
-def test_square_domain(R, metadata, method, inverse, expected):
-    """Test the square_domain."""
-    assert_array_equal(
-        dimension.square_domain(R, metadata, method, inverse)[0], expected
-    )
diff --git a/pysteps/tests/test_utils_images.py b/pysteps/tests/test_utils_images.py
new file mode 100644
index 000000000..d6c18ac30
--- /dev/null
+++ b/pysteps/tests/test_utils_images.py
@@ -0,0 +1,41 @@
+import numpy as np
+import xarray as xr
+
+from pysteps.utils import images
+
+
+def test_morph_opening_nans():
+    """"""
+    input_array = np.zeros((20, 20))
+    input_array[12:14, 12:14] = np.nan
+    input_array = xr.DataArray(input_array, dims=("y", "x"))
+    output_array = images.morph_opening(input_array, 0, 3)
+    assert isinstance(output_array, xr.DataArray)
+    assert output_array.sizes == input_array.sizes
+    assert np.isnan(input_array).sum() == 4
+    assert np.isnan(output_array).sum() == 4
+
+def test_morph_opening_2d():
+    """"""
+    input_array = np.zeros((20, 20))
+    input_array[1:11, 1:11] = 1
+    input_array[12:14, 12:14] = 1
+    input_array = xr.DataArray(input_array, dims=("y", "x"))
+    output_array = images.morph_opening(input_array, 0, np.ones((3, 3)))
+    assert isinstance(output_array, xr.DataArray)
+    assert output_array.sizes == input_array.sizes
+    assert int(input_array.sum()) == 104
+    assert int(output_array.sum()) == 100
+
+def test_morph_opening_3d():
+    """"""
+    input_array = np.zeros((3, 20, 20))
+    input_array[0, 1:11, 1:11] = 1
+    input_array[:, 12:14, 12:14] = 1
+    input_array = xr.DataArray(input_array, dims=("t", "y", "x"))
+    output_array = images.morph_opening(input_array, 0, np.ones((3, 3)))
+    assert isinstance(output_array, xr.DataArray)
+    assert output_array.sizes == input_array.sizes
+    assert int(input_array.sum()) == 112
+    assert int(output_array.sum()) == 100
+
diff --git a/pysteps/tests/test_utils_interpolate.py b/pysteps/tests/test_utils_interpolate.py
index eee0e6b2c..64d88c7c6 100644
--- a/pysteps/tests/test_utils_interpolate.py
+++ b/pysteps/tests/test_utils_interpolate.py
@@ -1,127 +1,175 @@
 # -*- coding: utf-8 -*-
+from pkg_resources import parse_version
+
+import xarray as xr
 import numpy as np
 import pytest
+from scipy import __version__ as scipy_version
 
+from pysteps.exceptions import MissingOptionalDependency
 from pysteps.utils import get_method
 
-
-interp_methods = (
-    "idwinterp2d",
-    "rbfinterp2d",
-)
+RBF_IMPORTED = parse_version(scipy_version) >= parse_version("1.7.0")
+
+interp_methods = ["idwinterp2d"]
+if RBF_IMPORTED:
+    interp_methods.append("rbfinterp2d")
+
+
+def make_dummy_input(n_vars, input_array=None):
+    xgrid, ygrid = np.linspace(0, 1, 10), np.linspace(0, 1, 8)
+    coord = np.random.rand(2, 10)
+    coord[-1, -1] = np.nan  # robust to missing coordinate ?
+    if input_array is None:
+        input_array = np.random.rand(10, n_vars)
+        input_array[-1, -1] = np.nan  # robust to missing values ?
+    sparse_data = xr.DataArray(
+        input_array.astype(np.float32),
+        dims=("sample", "variable"),
+        coords={
+            "x": ("sample", coord[0]),
+            "y": ("sample", coord[1]),
+            "variable": ("variable", list(map(lambda x: f"var{x+1}", range(n_vars)))),
+        },
+        attrs=dict(units="units"),
+    )
+    if n_vars == 1:
+        sparse_data = sparse_data.squeeze("variable")
+    return sparse_data, xgrid, ygrid
 
 
 @pytest.mark.parametrize("interp_method", interp_methods)
 def test_interp_univariate(interp_method):
-    coord = np.random.rand(10, 2)
-    input_array = np.random.rand(10)
-    xgrid, ygrid = np.linspace(0, 1, 10), np.linspace(0, 1, 10)
-
     interp = get_method(interp_method)
-    output = interp(coord, input_array, xgrid, ygrid)
+    sparse_data, xgrid, ygrid = make_dummy_input(1)
 
-    assert isinstance(output, np.ndarray)
+    # input a DataArray
+    output = interp(sparse_data, xgrid, ygrid)
+    assert isinstance(output, xr.DataArray)
     assert output.ndim == 2
     assert output.shape == (ygrid.size, xgrid.size)
     assert np.isfinite(output).all()
+    assert output.dtype == sparse_data.dtype
+    assert output.attrs == sparse_data.attrs
+    assert np.allclose(output.x, xgrid)
+    assert np.allclose(output.y, ygrid)
+
+    # input a Dataset
+    sparse_data = sparse_data.to_dataset(name="var1")
+    output = interp(sparse_data, xgrid, ygrid)
+    assert isinstance(output, xr.Dataset)
+    assert output.var1.ndim == 2
+    assert output.var1.shape == (ygrid.size, xgrid.size)
+    assert output.var1.dtype == sparse_data.var1.dtype
+    assert np.isfinite(output.var1).all()
+    assert output.attrs == sparse_data.attrs
+    assert np.allclose(output.x, xgrid)
+    assert np.allclose(output.y, ygrid)
 
 
 @pytest.mark.parametrize("interp_method", interp_methods)
 def test_interp_multivariate(interp_method):
-    coord = np.random.rand(10, 2)
-    input_array = np.random.rand(10, 2)
-    xgrid, ygrid = np.linspace(0, 1, 10), np.linspace(0, 1, 10)
-
     interp = get_method(interp_method)
-    output = interp(coord, input_array, xgrid, ygrid)
+    sparse_data, xgrid, ygrid = make_dummy_input(2)
 
-    assert isinstance(output, np.ndarray)
+    # input a DataArray
+    output = interp(sparse_data, xgrid, ygrid)
+    assert isinstance(output, xr.DataArray)
     assert output.ndim == 3
-    assert output.shape[0] == 2
-    assert output.shape[1:] == (ygrid.size, xgrid.size)
+    assert output.sizes["variable"] == 2
+    assert output.isel(variable=0).shape == (ygrid.size, xgrid.size)
     assert np.isfinite(output).all()
+    assert output.dtype == sparse_data.dtype
+    assert output.attrs == sparse_data.attrs
+    assert np.allclose(output.x, xgrid)
+    assert np.allclose(output.y, ygrid)
+
+    # input a Dataset
+    sparse_data = sparse_data.to_dataset(dim="variable")
+    output = interp(sparse_data, xgrid, ygrid)
+    assert isinstance(output, xr.Dataset)
+    assert "var1" in sparse_data.data_vars
+    assert "var2" in sparse_data.data_vars
+    assert output.var1.ndim == 2
+    assert output.var1.shape == (ygrid.size, xgrid.size)
+    assert output.var1.dtype == sparse_data.var1.dtype
+    assert np.isfinite(output.var1).all()
+    assert output.attrs == sparse_data.attrs
+    assert np.allclose(output.x, xgrid)
+    assert np.allclose(output.y, ygrid)
 
 
 @pytest.mark.parametrize("interp_method", interp_methods)
 def test_wrong_inputs(interp_method):
-    coord = np.random.rand(10, 2)
-    input_array = np.random.rand(10, 2)
-    xgrid, ygrid = np.linspace(0, 1, 10), np.linspace(0, 1, 10)
-
     interp = get_method(interp_method)
+    sparse_data, xgrid, ygrid = make_dummy_input(2)
 
-    # nan in the input values
+    # not an xarray object
     with pytest.raises(ValueError):
-        input_with_nans = input_array.copy()
-        input_with_nans[0, 0] = np.nan
-        interp(coord, input_with_nans, xgrid, ygrid)
+        interp(sparse_data.values, xgrid, ygrid)
 
-    # nan in the input coordinates
+    # missing 'sample' dimension
     with pytest.raises(ValueError):
-        coord_with_nans = coord.copy()
-        coord_with_nans[0, 0] = np.nan
-        interp(coord_with_nans, input_array, xgrid, ygrid)
+        interp(sparse_data.rename({"sample": "dummy"}), xgrid, ygrid)
 
     # too many dimensions in the input values
     with pytest.raises(ValueError):
-        interp(coord, np.random.rand(10, 2, 1), xgrid, ygrid)
-
-    # wrong dimension size in the input coordinates
-    with pytest.raises(ValueError):
-        interp(np.random.rand(10, 1), input_array, xgrid, ygrid)
+        interp(sparse_data.expand_dims("dummy"), xgrid, ygrid)
 
-    # wrong number of dimensions in the input coordinates
+    # missing x or y coordinates
     with pytest.raises(ValueError):
-        interp(np.random.rand(10, 2, 1), input_array, xgrid, ygrid)
+        interp(sparse_data.drop_vars("x"), xgrid, ygrid)
+        interp(sparse_data.drop_vars("y"), xgrid, ygrid)
 
-    # wrong number of coordinates
+    # missing values in target coordinates
     with pytest.raises(ValueError):
-        interp(np.random.rand(9, 2), input_array, xgrid, ygrid)
+        xgridn, ygridn = xgrid.copy(), ygrid.copy()
+        xgridn[0] = np.nan
+        ygridn[0] = np.nan
+        interp(sparse_data, xgridn, ygrid)
+        interp(sparse_data, xgrid, ygridn)
 
 
 @pytest.mark.parametrize("interp_method", interp_methods)
 def test_one_sample_input(interp_method):
-    coord = np.random.rand(1, 2)
-    input_array = np.array([1, 2])[None, :]
-    xgrid, ygrid = np.linspace(0, 1, 10), np.linspace(0, 1, 10)
-
     interp = get_method(interp_method)
+    sparse_data, xgrid, ygrid = make_dummy_input(2)
 
     # one sample returns uniform grids
-    output = interp(coord, input_array, xgrid, ygrid)
+    output = interp(sparse_data.isel(sample=slice(0, 1)), xgrid, ygrid)
+    assert isinstance(output, xr.DataArray)
+    assert output.ndim == 3
+    assert output.sizes["variable"] == 2
+    assert output.isel(variable=0).shape == (ygrid.size, xgrid.size)
     assert np.isfinite(output).all()
-    assert output[0, ...].max() == output[0, ...].min() == 1
-    assert output[1, ...].max() == output[1, ...].min() == 2
+    assert output.dtype == sparse_data.dtype
+    assert output.attrs == sparse_data.attrs
+    assert np.allclose(output.x, xgrid)
+    assert np.allclose(output.y, ygrid)
+    assert output[0, ...].max() == output[0, ...].min()
+    assert output[1, ...].max() == output[1, ...].min()
 
 
 @pytest.mark.parametrize("interp_method", interp_methods)
 def test_uniform_input(interp_method):
-    coord = np.random.rand(10, 2)
-    xgrid, ygrid = np.linspace(0, 1, 10), np.linspace(0, 1, 10)
-
     interp = get_method(interp_method)
 
     # same value across all variables
     input_array = np.ones((10, 2))
-    output = interp(coord, input_array, xgrid, ygrid)
+    sparse_data, xgrid, ygrid = make_dummy_input(2, input_array)
+    output = interp(sparse_data, xgrid, ygrid)
     assert np.isfinite(output).all()
     assert output.max() == output.min() == input_array.ravel()[0]
 
-    # # same value in one variable only
-    # input_array = np.vstack((np.ones(10), np.random.rand(10))).T
-    # output = interp(coord, input_array, xgrid, ygrid)
-    # assert output[0,].max() == output[0,].min() == input_array[0,0]
+    # TODO: test same value in one variable only
 
 
 def test_idwinterp2d_k1():
-    coord = np.random.rand(10, 2)
-    input_array = np.random.rand(10, 2)
-    xgrid, ygrid = np.linspace(0, 1, 10), np.linspace(0, 1, 10)
-
     interp = get_method("idwinterp2d")
-    output = interp(coord, input_array, xgrid, ygrid, k=1)
+    sparse_data, xgrid, ygrid = make_dummy_input(2)
 
-    assert isinstance(output, np.ndarray)
+    output = interp(sparse_data, xgrid, ygrid, k=1)
+    assert isinstance(output, xr.DataArray)
     assert output.ndim == 3
     assert output.shape[0] == 2
     assert output.shape[1:] == (ygrid.size, xgrid.size)
@@ -129,15 +177,21 @@ def test_idwinterp2d_k1():
 
 
 def test_idwinterp2d_kNone():
-    coord = np.random.rand(10, 2)
-    input_array = np.random.rand(10, 2)
-    xgrid, ygrid = np.linspace(0, 1, 10), np.linspace(0, 1, 10)
-
     interp = get_method("idwinterp2d")
-    output = interp(coord, input_array, xgrid, ygrid, k=None)
+    sparse_data, xgrid, ygrid = make_dummy_input(2)
 
-    assert isinstance(output, np.ndarray)
+    output = interp(sparse_data, xgrid, ygrid, k=None)
+    assert isinstance(output, xr.DataArray)
     assert output.ndim == 3
     assert output.shape[0] == 2
     assert output.shape[1:] == (ygrid.size, xgrid.size)
     assert np.isfinite(output).all()
+
+
+@pytest.mark.skipif(RBF_IMPORTED, reason="RBFInterpolator is available")
+def test_rbf_missingoptionaldependency():
+    """For scipy<1.7, RBFInterpolator is not available"""
+    interp = get_method("rbfinterp2d")
+    sparse_data, xgrid, ygrid = make_dummy_input(2)
+    with pytest.raises(MissingOptionalDependency):
+        interp(sparse_data, xgrid, ygrid)
diff --git a/pysteps/tests/test_utils_transformation.py b/pysteps/tests/test_utils_transformation.py
index 101e6b9d5..bc56a9d0a 100644
--- a/pysteps/tests/test_utils_transformation.py
+++ b/pysteps/tests/test_utils_transformation.py
@@ -2,189 +2,72 @@
 
 import numpy as np
 import pytest
-from numpy.testing import assert_array_almost_equal
-
-from pysteps.utils import transformation
-
-# boxcox_transform
-test_data = [
-    (
-        np.array([1]),
-        {
-            "accutime": 5,
-            "transform": None,
-            "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
-        },
-        None,
-        None,
-        None,
-        False,
-        np.array([0]),
-    ),
-    (
-        np.array([1]),
-        {
-            "accutime": 5,
-            "transform": "BoxCox",
-            "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
-        },
-        None,
-        None,
-        None,
-        True,
-        np.array([np.exp(1)]),
-    ),
-    (
-        np.array([1]),
-        {
-            "accutime": 5,
-            "transform": None,
-            "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
-        },
-        1.0,
-        None,
-        None,
-        False,
-        np.array([0]),
-    ),
-    (
-        np.array([1]),
-        {
-            "accutime": 5,
-            "transform": "BoxCox",
-            "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
-        },
-        1.0,
-        None,
-        None,
-        True,
-        np.array([2.0]),
-    ),
-]
-
-
-@pytest.mark.parametrize(
-    "R, metadata, Lambda, threshold, zerovalue, inverse, expected", test_data
-)
-def test_boxcox_transform(R, metadata, Lambda, threshold, zerovalue, inverse, expected):
-    """Test the boxcox_transform."""
-    assert_array_almost_equal(
-        transformation.boxcox_transform(
-            R, metadata, Lambda, threshold, zerovalue, inverse
-        )[0],
-        expected,
-    )
-
-
-# dB_transform
-test_data = [
-    (
-        np.array([1]),
-        {
-            "accutime": 5,
-            "transform": None,
-            "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
-        },
-        None,
-        None,
-        False,
-        np.array([0]),
-    ),
-    (
-        np.array([1]),
-        {
-            "accutime": 5,
-            "transform": "dB",
-            "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
-        },
-        None,
-        None,
-        True,
-        np.array([1.25892541]),
-    ),
-]
-
-
-@pytest.mark.parametrize(
-    "R, metadata, threshold, zerovalue, inverse, expected", test_data
+import xarray as xr
+
+
+xr.set_options(keep_attrs=True)
+
+
+data_array = xr.DataArray(
+    np.random.randn(10, 10),
+    coords={"x": np.arange(10), "y": np.arange(10)},
+    dims=("y", "x"),
+    attrs={"unit": "dummy"},
 )
-def test_dB_transform(R, metadata, threshold, zerovalue, inverse, expected):
-    """Test the dB_transform."""
-    assert_array_almost_equal(
-        transformation.dB_transform(R, metadata, threshold, zerovalue, inverse)[0],
-        expected,
-    )
-
-
-# NQ_transform
-test_data = [
-    (
-        np.array([1, 2]),
-        {
-            "accutime": 5,
-            "transform": None,
-            "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
-        },
-        False,
-        np.array([-0.4307273, 0.4307273]),
-    )
-]
-
-
-@pytest.mark.parametrize("R, metadata, inverse, expected", test_data)
-def test_NQ_transform(R, metadata, inverse, expected):
-    """Test the NQ_transform."""
-    assert_array_almost_equal(
-        transformation.NQ_transform(R, metadata, inverse)[0], expected
-    )
-
-
-# sqrt_transform
-test_data = [
-    (
-        np.array([1]),
-        {
-            "accutime": 5,
-            "transform": None,
-            "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
-        },
-        False,
-        np.array([1]),
-    ),
-    (
-        np.array([1]),
-        {
-            "accutime": 5,
-            "transform": "sqrt",
-            "unit": "mm/h",
-            "threshold": 0,
-            "zerovalue": 0,
-        },
-        True,
-        np.array([1]),
-    ),
-]
-
-
-@pytest.mark.parametrize("R, metadata, inverse, expected", test_data)
-def test_sqrt_transform(R, metadata, inverse, expected):
-    """Test the sqrt_transform."""
-    assert_array_almost_equal(
-        transformation.sqrt_transform(R, metadata, inverse)[0], expected
-    )
+data_array = data_array.where(data_array > -1)
+data_array = xr.where(data_array < 0, 0, data_array)
+
+
+test_args = [(1.0, 0.01), (0.0, 0.01)]
+
+
+@pytest.mark.parametrize("boxcox_lambda, offset", test_args)
+def test_boxcox_transform(boxcox_lambda, offset):
+
+    data_transformed = data_array.pysteps.boxcox_transform(boxcox_lambda, offset)
+    assert data_transformed.attrs.get("boxcox_lambda") == boxcox_lambda
+    assert data_transformed.attrs.get("offset") == offset
+    assert data_transformed.attrs.get("transform") == "BoxCox"
+    assert data_transformed.attrs.get("unit") == "dummy"
+
+    data_back = data_transformed.pysteps.boxcox_transform(inverse=True)
+    assert data_back.attrs.get("transform") is None
+    assert data_back.attrs.get("unit") == "dummy"
+    assert "offset" not in data_back.attrs
+    xr.testing.assert_allclose(data_back, data_array)
+
+
+def test_db_transform():
+    offset = 0.01
+    data_transformed = data_array.pysteps.db_transform(offset)
+    assert data_transformed.attrs.get("offset") == offset
+    assert data_transformed.attrs.get("transform") == "dB"
+
+    data_back = data_transformed.pysteps.db_transform(inverse=True)
+    assert data_back.attrs.get("transform") is None
+    assert data_back.attrs.get("unit") == "dummy"
+    assert "offset" not in data_back.attrs
+    xr.testing.assert_allclose(data_back, data_array)
+
+
+def test_nq_transform():
+    nq_a = 0.0
+    data_transformed = data_array.pysteps.nq_transform(nq_a)
+    assert data_transformed.attrs.get("transform") == "NQ"
+
+    data_back = data_transformed.pysteps.nq_transform(template=data_array, inverse=True)
+    assert data_back.attrs.get("transform") is None
+    assert data_back.attrs.get("unit") == "dummy"
+    xr.testing.assert_equal(np.isnan(data_back), np.isnan(data_array))
+    # we do not expect an exact match since the extremes are lost
+    xr.testing.assert_allclose(data_back, data_array, rtol=1, atol=0.1)
+
+
+def test_sqrt_transform():
+    data_transformed = data_array.pysteps.sqrt_transform()
+    assert data_transformed.attrs.get("transform") == "sqrt"
+
+    data_back = data_transformed.pysteps.sqrt_transform(inverse=True)
+    assert data_back.attrs.get("transform") is None
+    assert data_back.attrs.get("unit") == "dummy"
+    xr.testing.assert_allclose(data_back, data_array)