Repository: pySTEPS/pysteps (language: python)

Issue:
Implement importer for NWP (netcdf) data
Implement function to read in CF-compliant NetCDF data, not (yet) including interpolation to the radar / nowcasting grid. This might be simplified by moving to xarray first. Required for #212 .


Discussion / hints:
@ladc working on that. I have just added one example with radar data and NWP rainfall forecasts into https://github.com/pySTEPS/pysteps-data
I just saw that the nwp importers do not give the coordinates of the middle of the cell, but of the left cell edge instead. The "_to_xarray" function in the decorators is skipped with the nwp data, because it is read as an xarray DataArray and therefore skipped by:
`if not isinstance(data_array, xr.DataArray):`

Should we directly set everything to the cell centers in the nwp importers, or should we adjust the decorators in such a way that it also knows how to work with this data? 

Fix patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index c80571951..f35ad636b 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,6 +1,6 @@
 repos:
 -   repo: https://github.com/psf/black
-    rev: 21.6b0
+    rev: 21.7b0
     hooks:
     - id: black
       language_version: python3
\ No newline at end of file
diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst
index c8fbcd8d6..473297b14 100644
--- a/CONTRIBUTING.rst
+++ b/CONTRIBUTING.rst
@@ -231,7 +231,8 @@ Since PEP8 is a set of recommendations, these are the most important good coding
           </tr>
 
         </table>
-   (source: `Section 3.16.4, Google's python style guide <http://google.github.io/styleguide/pyguide.html>`__)
+
+(source: `Section 3.16.4, Google's python style guide <http://google.github.io/styleguide/pyguide.html>`__)
 
 - If you need to ignore part of the variables returned by a function,
   use "_" (single underscore) or __ (double underscore)::
diff --git a/doc/requirements.txt b/doc/requirements.txt
index 05e2450cd..2448201db 100644
--- a/doc/requirements.txt
+++ b/doc/requirements.txt
@@ -1,8 +1,8 @@
-# Additional requeriments related to the documentation build only
+# Additional requirements related to the documentation build only
 sphinx
 sphinxcontrib.bibtex
 sphinx-book-theme
 sphinx_gallery
 scikit-image
 pandas
-
+git+https://github.com/pySTEPS/pysteps-nwp-importers.git@main#egg=pysteps_nwp_importers
\ No newline at end of file
diff --git a/doc/source/pysteps_reference/blending.rst b/doc/source/pysteps_reference/blending.rst
new file mode 100644
index 000000000..35fe47c01
--- /dev/null
+++ b/doc/source/pysteps_reference/blending.rst
@@ -0,0 +1,12 @@
+================
+pysteps.blending
+================
+
+Implementation of blending methods for blending (ensemble) nowcasts with Numerical Weather Prediction (NWP) models.
+
+.. automodule:: pysteps.blending.interface
+.. automodule:: pysteps.blending.clim
+.. automodule:: pysteps.blending.linear_blending
+.. automodule:: pysteps.blending.skill_scores
+.. automodule:: pysteps.blending.steps
+.. automodule:: pysteps.blending.utils
diff --git a/doc/source/pysteps_reference/index.rst b/doc/source/pysteps_reference/index.rst
index f7b049578..a8c4c83df 100644
--- a/doc/source/pysteps_reference/index.rst
+++ b/doc/source/pysteps_reference/index.rst
@@ -14,6 +14,7 @@ available in pysteps.
     :caption: API Reference
 
     pysteps
+    blending
     cascade
     decorators
     extrapolation
diff --git a/doc/source/pysteps_reference/utils.rst b/doc/source/pysteps_reference/utils.rst
index 7e5474145..0622e6bbe 100644
--- a/doc/source/pysteps_reference/utils.rst
+++ b/doc/source/pysteps_reference/utils.rst
@@ -16,3 +16,4 @@ Implementation of miscellaneous utility functions.
 .. automodule:: pysteps.utils.spectral
 .. automodule:: pysteps.utils.tapering
 .. automodule:: pysteps.utils.transformation
+.. automodule:: pysteps.utils.reprojection
diff --git a/doc/source/references.bib b/doc/source/references.bib
index 69fafd3c7..af96f3d3e 100644
--- a/doc/source/references.bib
+++ b/doc/source/references.bib
@@ -1,4 +1,14 @@
 
+@TECHREPORT{BPS2004,
+  AUTHOR = "N. E. Bowler and C. E. Pierce and A. W. Seed",
+  TITLE = "{STEPS}: A probabilistic precipitation forecasting scheme which merges an extrapolation nowcast with downscaled {NWP}",
+  INSTITUTION = "UK Met Office",
+  TYPE = "Forecasting Research Technical Report",
+  NUMBER = 433,
+  ADDRESS = "Wallingford, United Kingdom",
+  YEAR = 2004,
+}
+
 @ARTICLE{BPS2006,
   AUTHOR = "N. E. Bowler and C. E. Pierce and A. W. Seed",
   TITLE = "{STEPS}: A probabilistic precipitation forecasting scheme which merges an extrapolation nowcast with downscaled {NWP}",
diff --git a/doc/source/user_guide/install_pysteps.rst b/doc/source/user_guide/install_pysteps.rst
index 46de6d598..f32f14fbd 100644
--- a/doc/source/user_guide/install_pysteps.rst
+++ b/doc/source/user_guide/install_pysteps.rst
@@ -36,7 +36,9 @@ Other optional dependencies include:
 * `pywavelets <https://pywavelets.readthedocs.io/en/latest/>`_
   (for intensity-scale verification)
 * `pandas <https://pandas.pydata.org/>`_ and
-  `scikit-image <https://scikit-image.org/>`_ (for the DATing and LINDA nowcast methods)
+  `scikit-image <https://scikit-image.org/>`_ (for advanced feature detection methods)
+* `rasterio <https://rasterio.readthedocs.io/en/latest/>`_ (for the reprojection module)
+
 
 **Important**: If you only want to use pysteps, you can continue reading below.
 But, if you want to contribute to pysteps or edit the package, you need to install
diff --git a/environment_dev.yml b/environment_dev.yml
index 378611d98..0625e73f4 100644
--- a/environment_dev.yml
+++ b/environment_dev.yml
@@ -30,3 +30,4 @@ dependencies:
   - cartopy>=0.18
   - scikit-image
   - pandas
+  - rasterio
diff --git a/examples/anvil_nowcast.py b/examples/anvil_nowcast.py
index 57802ad4e..b922bdca4 100644
--- a/examples/anvil_nowcast.py
+++ b/examples/anvil_nowcast.py
@@ -118,9 +118,7 @@
     date, root_path, path_fmt, fn_pattern, fn_ext, timestep=5, num_next_files=3
 )
 
-refobs_field, quality, metadata = io.read_timeseries(
-    filenames, importer, **importer_kwargs
-)
+refobs_field, _, metadata = io.read_timeseries(filenames, importer, **importer_kwargs)
 
 refobs_field, metadata = utils.to_rainrate(refobs_field[-1], metadata)
 refobs_field[refobs_field < 0.5] = 0.0
diff --git a/examples/blended_forecast.py b/examples/blended_forecast.py
new file mode 100644
index 000000000..0745b3ca0
--- /dev/null
+++ b/examples/blended_forecast.py
@@ -0,0 +1,230 @@
+# -*- coding: utf-8 -*-
+"""
+Blended forecast
+====================
+
+This tutorial shows how to construct a blended forecast from an ensemble nowcast
+using the STEPS approach and a Numerical Weather Prediction (NWP) rainfall
+forecast. The used datasets are from the Bureau of Meteorology, Australia.
+
+"""
+
+import os
+from datetime import datetime
+
+import numpy as np
+from matplotlib import pyplot as plt
+
+import pysteps
+from pysteps import io, rcparams, blending
+from pysteps.visualization import plot_precip_field
+
+
+################################################################################
+# Read the radar images and the NWP forecast
+# ------------------------------------------
+#
+# First, we import a sequence of 3 images of 10-minute radar composites
+# and the corresponding NWP rainfall forecast that was available at that time.
+#
+# You need the pysteps-data archive downloaded and the pystepsrc file
+# configured with the data_source paths pointing to data folders.
+# Additionally, the pysteps-nwp-importers plugin needs to be installed, see
+# https://github.com/pySTEPS/pysteps-nwp-importers.
+
+# Selected case
+date_radar = datetime.strptime("202010310400", "%Y%m%d%H%M")
+# The last NWP forecast was issued at 00:00
+date_nwp = datetime.strptime("202010310000", "%Y%m%d%H%M")
+radar_data_source = rcparams.data_sources["bom"]
+nwp_data_source = rcparams.data_sources["bom_nwp"]
+
+###############################################################################
+# Load the data from the archive
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+root_path = radar_data_source["root_path"]
+path_fmt = "prcp-c10/66/%Y/%m/%d"
+fn_pattern = "66_%Y%m%d_%H%M00.prcp-c10"
+fn_ext = radar_data_source["fn_ext"]
+importer_name = radar_data_source["importer"]
+importer_kwargs = radar_data_source["importer_kwargs"]
+timestep = 10.0
+
+# Find the radar files in the archive
+fns = io.find_by_date(
+    date_radar, root_path, path_fmt, fn_pattern, fn_ext, timestep, num_prev_files=2
+)
+
+# Read the radar composites
+importer = io.get_method(importer_name, "importer")
+radar_precip, _, radar_metadata = io.read_timeseries(fns, importer, **importer_kwargs)
+
+# Import the NWP data
+filename = os.path.join(
+    nwp_data_source["root_path"],
+    datetime.strftime(date_nwp, nwp_data_source["path_fmt"]),
+    datetime.strftime(date_nwp, nwp_data_source["fn_pattern"])
+    + "."
+    + nwp_data_source["fn_ext"],
+)
+
+nwp_importer = io.get_method("bom_nwp", "importer")
+nwp_precip, _, nwp_metadata = nwp_importer(filename)
+
+# Only keep the NWP forecasts from the last radar observation time (2020-10-31 04:00)
+# onwards
+
+nwp_precip = nwp_precip[24:43, :, :]
+
+
+################################################################################
+# Pre-processing steps
+# --------------------
+
+# Make sure the units are in mm/h
+converter = pysteps.utils.get_method("mm/h")
+radar_precip, radar_metadata = converter(radar_precip, radar_metadata)
+nwp_precip, nwp_metadata = converter(nwp_precip, nwp_metadata)
+
+# Threshold the data
+radar_precip[radar_precip < 0.1] = 0.0
+nwp_precip[nwp_precip < 0.1] = 0.0
+
+# Plot the radar rainfall field and the first time step of the NWP forecast.
+date_str = datetime.strftime(date_radar, "%Y-%m-%d %H:%M")
+plt.figure(figsize=(10, 5))
+plt.subplot(121)
+plot_precip_field(
+    radar_precip[-1, :, :],
+    geodata=radar_metadata,
+    title=f"Radar observation at {date_str}",
+)
+plt.subplot(122)
+plot_precip_field(
+    nwp_precip[0, :, :], geodata=nwp_metadata, title=f"NWP forecast at {date_str}"
+)
+plt.tight_layout()
+plt.show()
+
+# transform the data to dB
+transformer = pysteps.utils.get_method("dB")
+radar_precip, radar_metadata = transformer(radar_precip, radar_metadata, threshold=0.1)
+nwp_precip, nwp_metadata = transformer(nwp_precip, nwp_metadata, threshold=0.1)
+
+# r_nwp has to be four dimentional (n_models, time, y, x).
+# If we only use one model:
+if nwp_precip.ndim == 3:
+    nwp_precip = nwp_precip[None, :]
+
+###############################################################################
+# For the initial time step (t=0), the NWP rainfall forecast is not that different
+# from the observed radar rainfall, but it misses some of the locations and
+# shapes of the observed rainfall fields. Therefore, the NWP rainfall forecast will
+# initially get a low weight in the blending process.
+#
+# Determine the velocity fields
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+oflow_method = pysteps.motion.get_method("lucaskanade")
+
+# First for the radar images
+velocity_radar = oflow_method(radar_precip)
+
+# Then for the NWP forecast
+velocity_nwp = []
+# Loop through the models
+for n_model in range(nwp_precip.shape[0]):
+    # Loop through the timesteps. We need two images to construct a motion
+    # field, so we can start from timestep 1. Timestep 0 will be the same
+    # as timestep 1.
+    _v_nwp_ = []
+    for t in range(1, nwp_precip.shape[1]):
+        v_nwp_ = oflow_method(nwp_precip[n_model, t - 1 : t + 1, :])
+        _v_nwp_.append(v_nwp_)
+        v_nwp_ = None
+    # Add the velocity field at time step 1 to time step 0.
+    _v_nwp_ = np.insert(_v_nwp_, 0, _v_nwp_[0], axis=0)
+    velocity_nwp.append(_v_nwp_)
+velocity_nwp = np.stack(velocity_nwp)
+
+
+################################################################################
+# The blended forecast
+# --------------------
+
+precip_forecast = blending.steps.forecast(
+    precip=radar_precip,
+    precip_models=nwp_precip,
+    velocity=velocity_radar,
+    velocity_models=velocity_nwp,
+    timesteps=18,
+    timestep=timestep,
+    issuetime=date_radar,
+    n_ens_members=1,
+    precip_thr=radar_metadata["threshold"],
+    kmperpixel=radar_metadata["xpixelsize"] / 1000.0,
+    noise_stddev_adj="auto",
+    vel_pert_method=None,
+)
+
+# Transform the data back into mm/h
+precip_forecast, _ = converter(precip_forecast, radar_metadata)
+radar_precip, _ = converter(radar_precip, radar_metadata)
+nwp_precip, _ = converter(nwp_precip, nwp_metadata)
+
+
+################################################################################
+# Visualize the output
+# ~~~~~~~~~~~~~~~~~~~~
+#
+# The NWP rainfall forecast has a lower weight than the radar-based extrapolation
+# forecast at the issue time of the forecast (+0 min). Therefore, the first time
+# steps consist mostly of the extrapolation.
+# However, near the end of the forecast (+180 min), the NWP share in the blended
+# forecast has become more important and the forecast starts to resemble the
+# NWP forecast more.
+
+fig = plt.figure(figsize=(5, 12))
+
+leadtimes_min = [30, 60, 90, 120, 150, 180]
+n_leadtimes = len(leadtimes_min)
+for n, leadtime in enumerate(leadtimes_min):
+
+    # Nowcast with blending into NWP
+    plt.subplot(n_leadtimes, 2, n * 2 + 1)
+    plot_precip_field(
+        precip_forecast[0, int(leadtime / timestep) - 1, :, :],
+        geodata=radar_metadata,
+        title=f"Nowcast +{leadtime} min",
+        axis="off",
+        colorbar=False,
+    )
+
+    # Raw NWP forecast
+    plt.subplot(n_leadtimes, 2, n * 2 + 2)
+    plot_precip_field(
+        nwp_precip[0, int(leadtime / timestep) - 1, :, :],
+        geodata=nwp_metadata,
+        title=f"NWP +{leadtime} min",
+        axis="off",
+        colorbar=False,
+    )
+
+
+################################################################################
+# References
+# ~~~~~~~~~~
+#
+# Bowler, N. E., and C. E. Pierce, and A. W. Seed. 2004. "STEPS: A probabilistic
+# precipitation forecasting scheme which merges an extrapolation nowcast with
+# downscaled NWP." Forecasting Research Technical Report No. 433. Wallingford, UK.
+#
+# Bowler, N. E., and C. E. Pierce, and A. W. Seed. 2006. "STEPS: A probabilistic
+# precipitation forecasting scheme which merges an extrapolation nowcast with
+# downscaled NWP." Quarterly Journal of the Royal Meteorological Society 132(16):
+# 2127-2155. https://doi.org/10.1256/qj.04.100
+#
+# Seed, A. W., and C. E. Pierce, and K. Norman. 2013. "Formulation and evaluation
+# of a scale decomposition-based stochastic precipitation nowcast scheme." Water
+# Resources Research 49(10): 6624-664. https://doi.org/10.1002/wrcr.20536
diff --git a/examples/plot_linear_blending.py b/examples/plot_linear_blending.py
new file mode 100644
index 000000000..b5319ddf1
--- /dev/null
+++ b/examples/plot_linear_blending.py
@@ -0,0 +1,189 @@
+# -*- coding: utf-8 -*-
+
+"""
+Linear blending
+===============
+
+This tutorial shows how to construct a simple linear blending between a STEPS
+ensemble nowcast and a Numerical Weather Prediction (NWP) rainfall forecast. The
+used datasets are from the Bureau of Meteorology, Australia.
+"""
+
+import os
+from datetime import datetime
+
+from matplotlib import pyplot as plt
+
+import pysteps
+from pysteps import io, rcparams, blending
+from pysteps.visualization import plot_precip_field
+
+
+################################################################################
+# Read the radar images and the NWP forecast
+# ------------------------------------------
+#
+# First, we import a sequence of 3 images of 10-minute radar composites
+# and the corresponding NWP rainfall forecast that was available at that time.
+#
+# You need the pysteps-data archive downloaded and the pystepsrc file
+# configured with the data_source paths pointing to data folders.
+# Additionally, the pysteps-nwp-importers plugin needs to be installed, see
+# https://github.com/pySTEPS/pysteps-nwp-importers.
+
+# Selected case
+date_radar = datetime.strptime("202010310400", "%Y%m%d%H%M")
+# The last NWP forecast was issued at 00:00
+date_nwp = datetime.strptime("202010310000", "%Y%m%d%H%M")
+radar_data_source = rcparams.data_sources["bom"]
+nwp_data_source = rcparams.data_sources["bom_nwp"]
+
+###############################################################################
+# Load the data from the archive
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+root_path = radar_data_source["root_path"]
+path_fmt = "prcp-c10/66/%Y/%m/%d"
+fn_pattern = "66_%Y%m%d_%H%M00.prcp-c10"
+fn_ext = radar_data_source["fn_ext"]
+importer_name = radar_data_source["importer"]
+importer_kwargs = radar_data_source["importer_kwargs"]
+timestep = 10.0
+
+# Find the radar files in the archive
+fns = io.find_by_date(
+    date_radar, root_path, path_fmt, fn_pattern, fn_ext, timestep, num_prev_files=2
+)
+
+# Read the radar composites
+importer = io.get_method(importer_name, "importer")
+radar_precip, _, radar_metadata = io.read_timeseries(fns, importer, **importer_kwargs)
+
+# Import the NWP data
+filename = os.path.join(
+    nwp_data_source["root_path"],
+    datetime.strftime(date_nwp, nwp_data_source["path_fmt"]),
+    datetime.strftime(date_nwp, nwp_data_source["fn_pattern"])
+    + "."
+    + nwp_data_source["fn_ext"],
+)
+
+nwp_importer = io.get_method("bom_nwp", "importer")
+nwp_precip, _, nwp_metadata = nwp_importer(filename)
+
+# Only keep the NWP forecasts from the last radar observation time (2020-10-31 04:00)
+# End of the forecast is 18 time steps (+3 hours) in advance.
+precip_nwp = nwp_precip[24:43, :, :]
+
+
+################################################################################
+# Pre-processing steps
+# --------------------
+
+# Make sure the units are in mm/h
+converter = pysteps.utils.get_method("mm/h")
+radar_precip, radar_metadata = converter(radar_precip, radar_metadata)
+precip_nwp, nwp_metadata = converter(precip_nwp, nwp_metadata)
+
+# Threshold the data
+radar_precip[radar_precip < 0.1] = 0.0
+precip_nwp[precip_nwp < 0.1] = 0.0
+
+# Plot the radar rainfall field and the first time step of the NWP forecast.
+# For the initial time step (t=0), the NWP rainfall forecast is not that different
+# from the observed radar rainfall, but it misses some of the locations and
+# shapes of the observed rainfall fields. Therefore, the NWP rainfall forecast will
+# initially get a low weight in the blending process.
+date_str = datetime.strftime(date_radar, "%Y-%m-%d %H:%M")
+plt.figure(figsize=(10, 5))
+plt.subplot(121)
+plot_precip_field(
+    radar_precip[-1, :, :],
+    geodata=radar_metadata,
+    title=f"Radar observation at {date_str}",
+)
+plt.subplot(122)
+plot_precip_field(
+    precip_nwp[0, :, :], geodata=nwp_metadata, title=f"NWP forecast at {date_str}"
+)
+plt.tight_layout()
+plt.show()
+
+# Only keep the NWP forecasts from 2020-10-31 04:05 onwards, because the first
+# forecast lead time starts at 04:05.
+precip_nwp = precip_nwp[1:]
+
+# Transform the radar data to dB - this transformation is useful for the motion
+# field estimation and the subsequent nowcasts. The NWP forecast is not
+# transformed, because the linear blending code sets everything back in mm/h
+# after the nowcast.
+transformer = pysteps.utils.get_method("dB")
+radar_precip, radar_metadata = transformer(radar_precip, radar_metadata, threshold=0.1)
+
+
+################################################################################
+# Determine the velocity field for the radar rainfall nowcast
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+oflow_method = pysteps.motion.get_method("lucaskanade")
+velocity_radar = oflow_method(radar_precip)
+
+
+################################################################################
+# The linear blending of nowcast and NWP rainfall forecast
+# --------------------------------------------------------
+
+# Calculate the blended precipitation field
+precip_blended = blending.linear_blending.forecast(
+    precip=radar_precip[-1, :, :],
+    precip_metadata=radar_metadata,
+    velocity=velocity_radar,
+    timesteps=18,
+    timestep=10,
+    nowcast_method="extrapolation",  # simple advection nowcast
+    precip_nwp=precip_nwp,
+    precip_nwp_metadata=nwp_metadata,
+    start_blending=60,  # in minutes (this is an arbritrary choice)
+    end_blending=120,  # in minutes (this is an arbritrary choice)
+)
+
+
+################################################################################
+# Visualize the output
+# ~~~~~~~~~~~~~~~~~~~~
+#
+# The linear blending starts at 60 min, so during the first 60 minutes the
+# blended forecast only consists of the extrapolation forecast (consisting of an
+# extrapolation nowcast). Between 60 and 120 min, the NWP forecast gradually gets more
+# weight, whereas the extrapolation forecasts gradually gets less weight.
+# After 120 min, the blended forecast entirely consists of the NWP rainfall
+# forecast.
+
+fig = plt.figure(figsize=(4, 8))
+
+leadtimes_min = [30, 60, 90, 120]
+n_leadtimes = len(leadtimes_min)
+for n, leadtime in enumerate(leadtimes_min):
+
+    # Nowcast with blending into NWP
+    plt.subplot(n_leadtimes, 2, n * 2 + 1)
+    plot_precip_field(
+        precip_blended[int(leadtime / timestep) - 1, :, :],
+        geodata=radar_metadata,
+        title=f"Nowcast +{leadtime} min",
+        axis="off",
+        colorbar=False,
+    )
+
+    # Raw NWP forecast
+    plt.subplot(n_leadtimes, 2, n * 2 + 2)
+    plot_precip_field(
+        precip_nwp[int(leadtime / timestep) - 1, :, :],
+        geodata=nwp_metadata,
+        title=f"NWP +{leadtime} min",
+        axis="off",
+        colorbar=False,
+    )
+
+plt.tight_layout()
+plt.show()
diff --git a/pysteps/__init__.py b/pysteps/__init__.py
index 22bbac4aa..c58b9f17c 100644
--- a/pysteps/__init__.py
+++ b/pysteps/__init__.py
@@ -8,6 +8,7 @@
 from jsonschema import Draft4Validator
 
 # import subpackages
+from . import blending
 from . import cascade
 from . import datasets
 from . import decorators
diff --git a/pysteps/blending/__init__.py b/pysteps/blending/__init__.py
new file mode 100644
index 000000000..640d81d01
--- /dev/null
+++ b/pysteps/blending/__init__.py
@@ -0,0 +1,7 @@
+# -*- coding: utf-8 -*-
+"""Methods for blending NWP model(s) with nowcasts."""
+
+from pysteps.blending.interface import get_method
+from .clim import *
+from .skill_scores import *
+from .utils import *
diff --git a/pysteps/blending/clim.py b/pysteps/blending/clim.py
new file mode 100644
index 000000000..808e30a96
--- /dev/null
+++ b/pysteps/blending/clim.py
@@ -0,0 +1,218 @@
+"""
+pysteps.blending.clim
+=====================
+
+Module with methods to read, write and compute past and climatological NWP model
+skill scores. The module stores the average daily skill score for the past t days
+and updates it every day. The resulting average climatological skill score is
+the skill the NWP model skill regresses to during the blended forecast. If no
+climatological values are present, the default skill from :cite:`BPS2006` is used.
+
+.. autosummary::
+    :toctree: ../generated/
+
+    get_default_skill
+    save_skill
+    calc_clim_skill
+"""
+
+import pickle
+from pathlib import Path
+
+import numpy as np
+
+
+def get_default_skill(n_cascade_levels=8, n_models=1):
+    """
+    Get the default climatological skill values as given in :cite:`BPS2006`.
+    Take subset of n_cascade_levels or add entries with small values (1e-4) if
+    n_cascade_levels differs from 8.
+
+    Parameters
+    ----------
+    n_cascade_levels: int, optional
+      Number of cascade levels. Defaults to 8.
+    n_models: int, optional
+      Number of NWP models. Defaults to 1.
+
+    Returns
+    -------
+    default_skill: array-like
+      Array of shape [model, scale_level] containing the climatological skill
+      values.
+
+    """
+
+    default_skill = np.array(
+        [0.848, 0.537, 0.237, 0.065, 0.020, 0.0044, 0.0052, 0.0040]
+    )
+    n_skill = default_skill.shape[0]
+    if n_cascade_levels < n_skill:
+        default_skill = default_skill[0:n_cascade_levels]
+    elif n_cascade_levels > n_skill:
+        default_skill = np.append(
+            default_skill, np.repeat(1e-4, n_cascade_levels - n_skill)
+        )
+    return np.resize(default_skill, (n_models, n_cascade_levels))
+
+
+def save_skill(
+    current_skill,
+    validtime,
+    outdir_path,
+    n_models=1,
+    window_length=30,
+):
+    """
+    Add the current NWP skill to update today's daily average skill. If the day
+    is over, update the list of daily average skill covering a rolling window.
+
+    Parameters
+    ----------
+    current_skill: array-like
+      Array of shape [model, scale_level, ...]
+      containing the current skill of the different NWP models per cascade
+      level.
+    validtime: datetime
+      Datetime object containing the date and time for which the current
+      skill are valid.
+    outdir_path: string
+      Path to folder where the historical skill are stored. Defaults to
+      path_workdir from rcparams.
+    n_models: int, optional
+      Number of NWP models. Defaults to 1.
+    window_length: int, optional
+      Length of window (in days) of daily skill that should be retained.
+      Defaults to 30.
+
+    Returns
+    -------
+    None
+
+    """
+
+    n_cascade_levels = current_skill.shape[1]
+
+    # Load skill_today, a dictionary containing {mean_skill, n, last_validtime}
+    new_skill_today_file = False
+
+    skill_today_file = Path(outdir_path) / "NWP_skill_today.pkl"
+    if skill_today_file.is_file():
+        with open(skill_today_file, "rb") as f:
+            skill_today = pickle.load(f)
+        if skill_today["mean_skill"].shape != current_skill.shape:
+            new_skill_today_file = True
+    else:
+        new_skill_today_file = True
+
+    if new_skill_today_file:
+        skill_today = {
+            "mean_skill": np.copy(current_skill),
+            "n": 0,
+            "last_validtime": validtime,
+        }
+
+    # Load the past skill which is an array with dimensions day x model x scale_level
+    past_skill_file = Path(outdir_path) / "NWP_skill_window.npy"
+    past_skill = None
+    if past_skill_file.is_file():
+        past_skill = np.load(past_skill_file)
+    # First check if we have started a new day wrt the last written skill, in which
+    # case we should update the daily skill file and reset daily statistics.
+    if skill_today["last_validtime"].date() < validtime.date():
+        # Append skill to the list of the past X daily averages.
+        if (
+            past_skill is not None
+            and past_skill.shape[2] == n_cascade_levels
+            and past_skill.shape[1] == n_models
+        ):
+            past_skill = np.append(past_skill, [skill_today["mean_skill"]], axis=0)
+        else:
+            past_skill = np.array([skill_today["mean_skill"]])
+
+        # Remove oldest if the number of entries exceeds the window length.
+        if past_skill.shape[0] > window_length:
+            past_skill = np.delete(past_skill, 0, axis=0)
+        # FIXME also write out last_validtime.date() in this file?
+        # In that case it will need pickling or netcdf.
+        # Write out the past skill within the rolling window.
+        np.save(past_skill_file, past_skill)
+        # Reset statistics for today.
+        skill_today["n"] = 0
+        skill_today["mean_skill"] = 0
+
+    # Reset today's skill if needed and/or compute online average from the
+    # current skill using numerically stable algorithm
+    skill_today["n"] += 1
+    skill_today["mean_skill"] += (
+        current_skill - skill_today["mean_skill"]
+    ) / skill_today["n"]
+    skill_today["last_validtime"] = validtime
+    # Make path if path does not exist
+    skill_today_file.parent.mkdir(exist_ok=True, parents=True)
+    # Open and write to skill file
+    with open(skill_today_file, "wb") as f:
+        pickle.dump(skill_today, f)
+
+    return None
+
+
+def calc_clim_skill(
+    outdir_path,
+    n_cascade_levels=8,
+    n_models=1,
+    window_length=30,
+):
+    """
+    Return the climatological skill based on the daily average skill in the
+    rolling window. This is done using a geometric mean.
+
+    Parameters
+    ----------
+    n_cascade_levels: int, optional
+      Number of cascade levels.
+    outdir_path: string
+      Path to folder where the historical skill are stored. Defaults to
+      path_workdir from rcparams.
+    n_models: int, optional
+      Number of NWP models. Defaults to 1.
+    window_length: int, optional
+      Length of window (in days) over which to compute the climatological
+      skill. Defaults to 30.
+
+    Returns
+    -------
+    climatological_mean_skill: array-like
+      Array of shape [model, scale_level, ...] containing the climatological
+      (geometric) mean skill.
+
+    """
+    past_skill_file = Path(outdir_path) / "NWP_skill_window.npy"
+    # past_skill has dimensions date x model x scale_level  x ....
+    if past_skill_file.is_file():
+        past_skill = np.load(past_skill_file)
+    else:
+        past_skill = np.array(None)
+    # check if there is enough data to compute the climatological skill
+    if not past_skill.any():
+        return get_default_skill(n_cascade_levels, n_models)
+    elif past_skill.shape[0] < window_length:
+        return get_default_skill(n_cascade_levels, n_models)
+    # reduce window if necessary
+    else:
+        past_skill = past_skill[-window_length:]
+
+    # Make sure past_skill cannot be lower than 10e-5
+    past_skill = np.where(past_skill < 10e-5, 10e-5, past_skill)
+
+    # Calculate climatological skill from the past_skill using the
+    # geometric mean.
+    geomean_skill = np.exp(np.log(past_skill).mean(axis=0))
+
+    # Make sure skill is always a positive value and a finite value
+    geomean_skill = np.where(geomean_skill < 10e-5, 10e-5, geomean_skill)
+    geomean_skill = np.nan_to_num(
+        geomean_skill, copy=True, nan=10e-5, posinf=10e-5, neginf=10e-5
+    )
+
+    return geomean_skill
diff --git a/pysteps/blending/interface.py b/pysteps/blending/interface.py
new file mode 100644
index 000000000..ed8046f1d
--- /dev/null
+++ b/pysteps/blending/interface.py
@@ -0,0 +1,57 @@
+# -*- coding: utf-8 -*-
+"""
+pysteps.blending.interface
+==========================
+Interface for the blending module. It returns a callable function for computing
+blended nowcasts with NWP models.
+
+.. autosummary::
+    :toctree: ../generated/
+
+    get_method
+"""
+
+from pysteps.blending import linear_blending
+from pysteps.blending import steps
+
+_blending_methods = dict()
+_blending_methods["linear_blending"] = linear_blending.forecast
+_blending_methods["steps"] = steps.forecast
+
+
+def get_method(name):
+    """
+    Return a callable function for computing nowcasts blending into an NWP
+    forecast.
+
+    Implemented methods:
+
+    +------------------+------------------------------------------------------+
+    |     Name         |              Description                             |
+    +==================+======================================================+
+    | linear_blending  | the linear blending of a nowcast method with other   |
+    |                  | data (e.g. NWP data).                                |
+    +------------------+------------------------------------------------------+
+    | steps            | the STEPS stochastic nowcasting blending method      |
+    |                  | described in :cite:`Seed2003`, :cite:`BPS2006` and   |
+    |                  | :cite:`SPN2013`. The blending weights approach       |
+    |                  | currently follows :cite:`BPS2006`.                   |
+    +------------------+------------------------------------------------------+
+    """
+    if isinstance(name, str):
+        name = name.lower()
+    else:
+        raise TypeError(
+            "Only strings supported for the method's names.\n"
+            + "Available names:"
+            + str(list(_blending_methods.keys()))
+        ) from None
+
+    try:
+        return _blending_methods[name]
+    except KeyError:
+        raise ValueError(
+            f"Unknown blending method {name}."
+            "The available methods are: "
+            f"{*list(_blending_methods.keys()),}"
+        ) from None
diff --git a/pysteps/blending/linear_blending.py b/pysteps/blending/linear_blending.py
new file mode 100644
index 000000000..d2835d70c
--- /dev/null
+++ b/pysteps/blending/linear_blending.py
@@ -0,0 +1,210 @@
+# -*- coding: utf-8 -*-
+"""
+pysteps.nowcasts.linear_blending
+================================
+
+Linear blending method to blend (ensemble) nowcasts and (ensemble) Numerical
+Weather Prediction (NWP) forecasts. The method uses a predefined start and end
+time for the blending. Before this start time, the resulting blended forecasts only
+consists of the nowcast(s). In between the start and end time, the nowcast(s)
+weight decreases and NWP forecasts weight increases linearly from 1(0) to
+0(1). After the end time, the blended forecast entirely consists of the NWP
+forecasts.
+
+Implementation of the linear blending between nowcast and NWP data.
+
+.. autosummary::
+    :toctree: ../generated/
+
+    forecast
+"""
+
+import numpy as np
+from pysteps import nowcasts
+from pysteps.utils import conversion
+
+
+def forecast(
+    precip,
+    precip_metadata,
+    velocity,
+    timesteps,
+    timestep,
+    nowcast_method,
+    precip_nwp=None,
+    precip_nwp_metadata=None,
+    start_blending=120,
+    end_blending=240,
+    fill_nwp=True,
+    nowcast_kwargs=None,
+):
+
+    """Generate a forecast by linearly blending nowcasts with NWP data
+
+    Parameters
+    ----------
+    precip: array_like
+      Array containing the input precipitation field(s) ordered by timestamp
+      from oldest to newest. The time steps between the inputs are assumed
+      to be regular.
+    precip_metadata: dict
+        Metadata dictionary containing (at least) the transform, unit and threshold
+        attributes as described in the documentation of :py:mod:`pysteps.io.importers`.
+    velocity; array_like
+      Array of shape (2, m, n) containing the x- and y-components of the advection
+      field. The velocities are assumed to represent one time step between the
+      inputs. All values are required to be finite.
+    timesteps: int
+      Number of time steps to forecast.
+    timestep: int or float
+      The time difference (in minutes) between consecutive forecast fields.
+    nowcast_method: str
+      Name of the nowcasting method. See :py:mod:`pysteps.nowcasts.interface`
+      for the list of available methods.
+    precip_nwp: array_like or NoneType, optional
+      Array of shape (timesteps, m, n) in the case of no ensemble or
+      of shape (n_ens_members, timesteps, m, n) in the case of an ensemble
+      containing the NWP precipitation fields ordered by timestamp from oldest
+      to newest. The time steps between the inputs are assumed to be regular
+      (and identical to the time step between the nowcasts). If no NWP
+      data is given the value of precip_nwp is None and no blending will be performed.
+    precip_nwp_metadata: dict or NoneType, optional
+        NWP metadata dictionary containing (at least) the transform, unit and threshold
+        attributes as described in the documentation of :py:mod:`pysteps.io.importers`.
+    start_blending: int, optional
+      Time stamp (in minutes) after which the blending should start. Before this
+      only the nowcast data is used.
+    end_blending: int, optional
+      Time stamp (in minutes) after which the blending should end. Between
+      start_blending and end_blending the nowcasts and NWP data are linearly
+      merged with each other. After end_blending only the NWP data is used.
+    fill_nwp: bool, optional
+      Standard value is True. If True, the NWP data will be used to fill in the
+      no data mask of the nowcast.
+    nowcast_kwargs: dict, optional
+      Dictionary containing keyword arguments for the nowcast method.
+
+    Returns
+    -------
+    R_blended: ndarray
+      Array of shape (timesteps, m, n) in the case of no ensemble or
+      of shape (n_ens_members, timesteps, m, n) in the case of an ensemble
+      containing the precipation forecast generated by linearly blending
+      the nowcasts and the NWP data. n_ens_members equals the maximum no. of
+      ensemble members in either the nowcast or nwp model(s).
+    """
+
+    if nowcast_kwargs is None:
+        nowcast_kwargs = dict()
+
+    # Calculate the nowcasts
+    nowcast_method_func = nowcasts.get_method(nowcast_method)
+    precip_nowcast = nowcast_method_func(
+        precip,
+        velocity,
+        timesteps,
+        **nowcast_kwargs,
+    )
+
+    # Make sure that precip_nowcast and precip_nwp are in mm/h
+    precip_nowcast, _ = conversion.to_rainrate(precip_nowcast, metadata=precip_metadata)
+
+    # Check if NWP data is given as input
+    if precip_nwp is not None:
+        precip_nwp, _ = conversion.to_rainrate(precip_nwp, metadata=precip_nwp_metadata)
+
+        if len(precip_nowcast.shape) == 4:
+            n_ens_members_nowcast = precip_nowcast.shape[0]
+            if n_ens_members_nowcast == 1:
+                precip_nowcast = np.squeeze(precip_nowcast)
+        else:
+            n_ens_members_nowcast = 1
+
+        if len(precip_nwp.shape) == 4:
+            n_ens_members_nwp = precip_nwp.shape[0]
+            if n_ens_members_nwp == 1:
+                precip_nwp = np.squeeze(precip_nwp)
+        else:
+            n_ens_members_nwp = 1
+
+        # Now, repeat the nowcast ensemble members or the nwp models/members until
+        # it has the same amount of members as n_ens_members_max. For instance, if
+        # you have 10 ensemble nowcasts members and 3 NWP members, the output will
+        # be an ensemble of 10 members. Hence, the three NWP members are blended
+        # with the first three members of the nowcast (member one with member one,
+        # two with two, etc.), subsequently, the same NWP members are blended with
+        # the next three members (NWP member one with member 4, NWP member 2 with
+        # member 5, etc.), until 10 is reached.
+        n_ens_members_max = max(n_ens_members_nowcast, n_ens_members_nwp)
+        n_ens_members_min = min(n_ens_members_nowcast, n_ens_members_nwp)
+
+        if n_ens_members_min != n_ens_members_max:
+            if n_ens_members_nwp == 1:
+                precip_nwp = np.repeat(
+                    precip_nwp[np.newaxis, :, :], n_ens_members_max, axis=0
+                )
+            elif n_ens_members_nowcast == 1:
+                precip_nowcast = np.repeat(
+                    precip_nowcast[np.newaxis, :, :], n_ens_members_max, axis=0
+                )
+            else:
+                repeats = [
+                    (n_ens_members_max + i) // n_ens_members_min
+                    for i in range(n_ens_members_min)
+                ]
+
+                if n_ens_members_nwp == n_ens_members_min:
+                    precip_nwp = np.repeat(precip_nwp, repeats, axis=0)
+                elif n_ens_members_nowcast == n_ens_members_min:
+                    precip_nowcast = np.repeat(precip_nowcast, repeats, axis=0)
+
+        # Check if dimensions are correct
+        assert (
+            precip_nwp.shape == precip_nowcast.shape
+        ), "The dimensions of precip_nowcast and precip_nwp need to be identical: dimension of precip_nwp = {} and dimension of precip_nowcast = {}".format(
+            precip_nwp.shape, precip_nowcast.shape
+        )
+
+        # Initialise output
+        R_blended = np.zeros_like(precip_nowcast)
+
+        # Calculate the weights
+        for i in range(timesteps):
+            # Calculate what time we are at
+            t = (i + 1) * timestep
+
+            # Calculate the weight with a linear relation (weight_nwp at start_blending = 0.0)
+            # and (weight_nwp at end_blending = 1.0)
+            weight_nwp = (t - start_blending) / (end_blending - start_blending)
+
+            # Set weights at times before start_blending and after end_blending
+            if weight_nwp < 0.0:
+                weight_nwp = 0.0
+            elif weight_nwp > 1.0:
+                weight_nwp = 1.0
+
+            # Calculate weight_nowcast
+            weight_nowcast = 1.0 - weight_nwp
+
+            # Calculate output by combining precip_nwp and precip_nowcast,
+            # while distinguishing between ensemble and non-ensemble methods
+            if n_ens_members_max == 1:
+                R_blended[i, :, :] = (
+                    weight_nwp * precip_nwp[i, :, :]
+                    + weight_nowcast * precip_nowcast[i, :, :]
+                )
+            else:
+                R_blended[:, i, :, :] = (
+                    weight_nwp * precip_nwp[:, i, :, :]
+                    + weight_nowcast * precip_nowcast[:, i, :, :]
+                )
+
+            # Find where the NaN values are and replace them with NWP data
+            if fill_nwp:
+                nan_indices = np.isnan(R_blended)
+                R_blended[nan_indices] = precip_nwp[nan_indices]
+    else:
+        # If no NWP data is given, the blended field is simply equal to the nowcast field
+        R_blended = precip_nowcast
+
+    return R_blended
diff --git a/pysteps/blending/skill_scores.py b/pysteps/blending/skill_scores.py
new file mode 100644
index 000000000..51447ca11
--- /dev/null
+++ b/pysteps/blending/skill_scores.py
@@ -0,0 +1,289 @@
+# -*- coding: utf-8 -*-
+"""
+pysteps.blending.skill_scores
+==============================
+
+Methods for computing skill scores, needed for the blending weights, of two-
+dimensional model fields with the latest observation field.
+
+.. autosummary::
+    :toctree: ../generated/
+
+    spatial_correlation
+    lt_dependent_cor_nwp
+    lt_dependent_cor_extrapolation
+    clim_regr_values
+"""
+
+import numpy as np
+from pysteps.blending import clim
+
+
+def spatial_correlation(obs, mod, domain_mask):
+    """Determine the spatial correlation between the cascade of the latest
+    available observed (radar) rainfall field and a time-synchronous cascade
+    derived from a model (generally NWP) field. Both fields are assumed to use
+    the same grid.
+
+
+    Parameters
+    ----------
+    obs : array-like
+        Array of shape [cascade_level, y, x] with per cascade_level the
+        normalized cascade of the observed (radar) rainfall field.
+    mod : array-like
+        Array of shape [cascade_level, y, x] with per cascade_level the
+        normalized cascade of the model field.
+    domain_mask : array-like
+        Boolean array of shape [y, x] indicating which cells fall outside the
+        radar domain.
+
+    Returns
+    -------
+    rho : array-like
+        Array of shape [n_cascade_levels] containing per cascade_level the
+        correlation between the normalized cascade of the observed (radar)
+        rainfall field and the normalized cascade of the model field.
+
+    References
+    ----------
+    :cite:`BPS2006`
+    :cite:`SPN2013`
+
+    """
+    rho = []
+    # Fill rho per cascade level, so loop through the cascade levels
+    for cascade_level in range(0, obs.shape[0]):
+        # Only calculate the skill for the pixels within the radar domain
+        # (as that is where there are observations)
+        obs_cascade_level = obs[cascade_level, :, :]
+        mod_cascade_level = mod[cascade_level, :, :]
+        obs_cascade_level[domain_mask] = np.nan
+        mod_cascade_level[domain_mask] = np.nan
+
+        # Flatten both arrays
+        obs_1d = obs_cascade_level.flatten()
+        mod_1d = mod_cascade_level.flatten()
+        # Calculate the correlation between the two
+        cov = np.nansum(
+            (mod_1d - np.nanmean(mod_1d)) * (obs_1d - np.nanmean(obs_1d))
+        )  # Without 1/n, as this cancels out (same for stdx and -y)
+        std_obs = np.sqrt(np.nansum((obs_1d - np.nanmean(obs_1d)) ** 2.0))
+        std_mod = np.sqrt(np.nansum((mod_1d - np.nanmean(mod_1d)) ** 2.0))
+        rho.append(cov / (std_mod * std_obs))
+
+    # Make sure rho is always a (finite) number
+    rho = np.nan_to_num(rho, copy=True, nan=10e-5, posinf=10e-5, neginf=10e-5)
+
+    return rho
+
+
+def lt_dependent_cor_nwp(lt, correlations, outdir_path, n_model=0, skill_kwargs=None):
+    """Determine the correlation of a model field for lead time lt and
+    cascade k, by assuming that the correlation determined at t=0 regresses
+    towards the climatological values.
+
+
+    Parameters
+    ----------
+    lt : int
+        The lead time of the forecast in minutes.
+    correlations : array-like
+        Array of shape [n_cascade_levels] containing per cascade_level the
+        correlation between the normalized cascade of the observed (radar)
+        rainfall field and the normalized cascade of the model field.
+    outdir_path: string
+      Path to folder where the historical skill are stored. Defaults to
+      path_workdir from rcparams.
+    n_model: int, optional
+        The index number of the (NWP) model when the climatological skill of
+        multiple (NWP) models is stored. For calculations with one model, or
+        when n_model is not provided, n_model = 0.
+    skill_kwargs : dict, optional
+        Dictionary containing e.g. the outdir_path, nmodels and window_length
+        parameters.
+
+    Returns
+    -------
+    rho : array-like
+        Array of shape [n_cascade_levels] containing, for lead time lt, per
+        cascade_level the correlation between the normalized cascade of the
+        observed (radar) rainfall field and the normalized cascade of the
+        model field.
+
+    References
+    ----------
+    :cite:`BPS2004`
+    :cite:`BPS2006`
+    """
+
+    if skill_kwargs is None:
+        skill_kwargs = dict()
+
+    # Obtain the climatological values towards which the correlations will
+    # regress
+    clim_cor_values, regr_pars = clim_regr_values(
+        n_cascade_levels=len(correlations),
+        outdir_path=outdir_path,
+        n_model=n_model,
+        skill_kwargs=skill_kwargs,
+    )
+    # Determine the speed of the regression (eq. 24 in BPS2004)
+    qm = np.exp(-lt / regr_pars[0, :]) * (2 - np.exp(-lt / regr_pars[1, :]))
+    # Determine the correlation for lead time lt
+    rho = qm * correlations + (1 - qm) * clim_cor_values
+
+    return rho
+
+
+def lt_dependent_cor_extrapolation(PHI, correlations=None, correlations_prev=None):
+    """Determine the correlation of the extrapolation (nowcast) component for
+    lead time lt and cascade k, by assuming that the correlation determined at
+    t=0 regresses towards the climatological values.
+
+
+    Parameters
+    ----------
+    PHI : array-like
+        Array of shape [n_cascade_levels, ar_order + 1] containing per
+        cascade level the autoregression parameters.
+    correlations : array-like, optional
+        Array of shape [n_cascade_levels] containing per cascade_level the
+        latest available correlation from the extrapolation component that can
+        be found from the AR-2 model.
+    correlations_prev : array-like, optional
+        Similar to correlations, but from the timestep before that.
+
+    Returns
+    -------
+    rho : array-like
+        Array of shape [n_cascade_levels] containing, for lead time lt, per
+        cascade_level the correlation of the extrapolation component.
+
+    References
+    ----------
+    :cite:`BPS2004`
+    :cite:`BPS2006`
+
+    """
+    # Check if correlations_prev exists, if not, we set it to 1.0
+    if correlations_prev is None:
+        correlations_prev = np.repeat(1.0, PHI.shape[0])
+    # Same for correlations at first time step, we set it to
+    # phi1 / (1 - phi2), see BPS2004
+    if correlations is None:
+        correlations = PHI[:, 0] / (1.0 - PHI[:, 1])
+
+    # Calculate the correlation for lead time lt
+    rho = PHI[:, 0] * correlations + PHI[:, 1] * correlations_prev
+
+    # Finally, set the current correlations array as the previous one for the
+    # next time step
+    rho_prev = correlations
+
+    return rho, rho_prev
+
+
+def clim_regr_values(n_cascade_levels, outdir_path, n_model=0, skill_kwargs=None):
+    """Obtains the climatological correlation values and regression parameters
+    from a file called NWP_weights_window.bin in the outdir_path. If this file
+    is not present yet, the values from :cite:`BPS2004` are used.
+
+
+    Parameters
+    ----------
+    n_cascade_levels : int
+        The number of cascade levels to use.
+    outdir_path: string
+      Path to folder where the historical skill are stored. Defaults to
+      path_workdir from rcparams.
+    n_model: int, optional
+        The index number of the (NWP) model when the climatological skill of
+        multiple (NWP) models is stored. For calculations with one model, or
+        when n_model is not provided, n_model = 0.
+    skill_kwargs : dict, optional
+        Dictionary containing e.g. the outdir_path, nmodels and window_length
+        parameters.
+
+    Returns
+    -------
+    clim_cor_values : array-like
+        Array of shape [n_cascade_levels] containing the
+        climatological values of the lag 1 and lag 2 auto-correlation
+        coefficients, obtained by calling a method implemented in
+        pysteps.blending.skill_scores.get_clim_skill_scores.
+    regr_pars : array-like
+        Array of shape [2, n_cascade_levels] containing the regression
+        parameters. These are fixed values that should be hard-coded in this
+        function. Unless changed by the user, the standard values from
+        `BPS2004` are used.
+
+    Notes
+    -----
+    The literature climatological values assume 8 cascade levels. In case
+    less or more cascade levels are used, the script will handle this by taking
+    the first n values or extending the array with a small value. This is not
+    ideal, but will be fixed once the clim_regr_file is made. Hence, this
+    requires a warm-up period of the model.
+    In addition, the regression parameter values (eq. 24 in BPS2004) are hard-
+    coded and can only be optimized by the user after (re)fitting of the
+    equation.
+
+    """
+
+    if skill_kwargs is None:
+        skill_kwargs = dict()
+
+    # First, obtain climatological skill values
+    try:
+        clim_cor_values = clim.calc_clim_skill(
+            outdir_path=outdir_path, n_cascade_levels=n_cascade_levels, **skill_kwargs
+        )
+    except FileNotFoundError:
+        # The climatological skill values file does not exist yet, so we'll
+        # use the default values from BPS2004.
+        clim_cor_values = np.array(
+            [[0.848, 0.537, 0.237, 0.065, 0.020, 0.0044, 0.0052, 0.0040]]
+        )
+
+    clim_cor_values = clim_cor_values[n_model, :]
+
+    # Check if clim_cor_values has only one model, otherwise it has
+    # returned the skill values for multiple models
+    if clim_cor_values.ndim != 1:
+        raise IndexError(
+            "The climatological cor. values of multiple models were returned, but only one model should be specified. Make sure to pass the argument nmodels in the function"
+        )
+
+    # Also check whether the number of cascade_levels is correct
+    if clim_cor_values.shape[0] > n_cascade_levels:
+        clim_cor_values = clim_cor_values[0:n_cascade_levels]
+    elif clim_cor_values.shape[0] < n_cascade_levels:
+        # Get the number of cascade levels that is missing
+        n_extra_lev = n_cascade_levels - clim_cor_values.shape[0]
+        # Append the array with correlation values of 10e-4
+        clim_cor_values = np.append(clim_cor_values, np.repeat(1e-4, n_extra_lev))
+
+    # Get the regression parameters (L in eq. 24 in BPS2004) - Hard coded,
+    # change to own values when present.
+    regr_pars = np.array(
+        [
+            [130.0, 165.0, 120.0, 55.0, 50.0, 15.0, 15.0, 10.0],
+            [155.0, 220.0, 200.0, 75.0, 10e4, 10e4, 10e4, 10e4],
+        ]
+    )
+
+    # Check whether the number of cascade_levels is correct
+    if regr_pars.shape[1] > n_cascade_levels:
+        regr_pars = regr_pars[:, 0:n_cascade_levels]
+    elif regr_pars.shape[1] < n_cascade_levels:
+        # Get the number of cascade levels that is missing
+        n_extra_lev = n_cascade_levels - regr_pars.shape[1]
+        # Append the array with correlation values of 10e-4
+        regr_pars = np.append(
+            regr_pars,
+            [np.repeat(10.0, n_extra_lev), np.repeat(10e4, n_extra_lev)],
+            axis=1,
+        )
+
+    return clim_cor_values, regr_pars
diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py
new file mode 100644
index 000000000..3cef7257e
--- /dev/null
+++ b/pysteps/blending/steps.py
@@ -0,0 +1,2039 @@
+# -*- coding: utf-8 -*-
+"""
+pysteps.blending.steps
+======================
+
+Implementation of the STEPS stochastic blending method as described in
+:cite:`BPS2004`, :cite:`BPS2006` and :cite:`SPN2013`. The STEPS blending method
+consists of the following main steps:
+
+    #. Set the radar rainfall fields in a Lagrangian space.
+    #. Initialize the noise method.
+    #. Perform the cascade decomposition for the input radar rainfall fields.
+       The method assumes that the cascade decomposition of the NWP model fields is
+       already done prior to calling the function, as the NWP model fields are
+       generally not updated with the same frequency (which is more efficient). A
+       method to decompose and store the NWP model fields whenever a new NWP model
+       field is present, is present in pysteps.blending.utils.decompose_NWP.
+    #. Estimate AR parameters for the extrapolation nowcast and noise cascade.
+    #. Before starting the forecast loop, determine which NWP models will be
+       combined with which nowcast ensemble members. The number of output ensemble
+       members equals the maximum number of (ensemble) members in the input, which
+       can be either the defined number of (nowcast) ensemble members or the number
+       of NWP models/members.
+    #. Initialize all the random generators.
+    #. Calculate the initial skill of the NWP model forecasts at t=0.
+    #. Start the forecasting loop:
+        #. Determine the skill and weights of the forecasting components
+           (extrapolation, NWP and noise) for that lead time.
+        #. Regress the extrapolation and noise cascades separately to the subsequent
+           time step.
+        #. Extrapolate the extrapolation and noise cascades to the current time step.
+        #. Blend the cascades.
+        #. Recompose the cascade to a rainfall field.
+        #. Post-processing steps (masking and probability matching, which are
+           different from the original blended STEPS implementation).
+
+.. autosummary::
+    :toctree: ../generated/
+
+    forecast
+    calculate_ratios
+    calculate_weights_bps
+    calculate_weights_spn
+    blend_means_sigmas
+"""
+
+import time
+
+import numpy as np
+import scipy.ndimage
+
+from pysteps import cascade
+from pysteps import extrapolation
+from pysteps import noise
+from pysteps import utils
+from pysteps.nowcasts import utils as nowcast_utils
+from pysteps.postprocessing import probmatching
+from pysteps.timeseries import autoregression, correlation
+from pysteps import blending
+
+try:
+    import dask
+
+    DASK_IMPORTED = True
+except ImportError:
+    DASK_IMPORTED = False
+
+
+def forecast(
+    precip,
+    precip_models,
+    velocity,
+    velocity_models,
+    timesteps,
+    timestep,
+    issuetime,
+    n_ens_members,
+    n_cascade_levels=8,
+    blend_nwp_members=False,
+    precip_thr=None,
+    kmperpixel=None,
+    extrap_method="semilagrangian",
+    decomp_method="fft",
+    bandpass_filter_method="gaussian",
+    noise_method="nonparametric",
+    noise_stddev_adj=None,
+    ar_order=2,
+    vel_pert_method="bps",
+    weights_method="bps",
+    conditional=False,
+    probmatching_method="cdf",
+    mask_method="incremental",
+    callback=None,
+    return_output=True,
+    seed=None,
+    num_workers=1,
+    fft_method="numpy",
+    domain="spatial",
+    outdir_path_skill="./tmp/",
+    extrap_kwargs=None,
+    filter_kwargs=None,
+    noise_kwargs=None,
+    vel_pert_kwargs=None,
+    clim_kwargs=None,
+    mask_kwargs=None,
+    measure_time=False,
+):
+    """
+    Generate a blended nowcast ensemble by using the Short-Term Ensemble
+    Prediction System (STEPS) method.
+
+    Parameters
+    ----------
+    precip: array-like
+      Array of shape (ar_order+1,m,n) containing the input precipitation fields
+      ordered by timestamp from oldest to newest. The time steps between the
+      inputs are assumed to be regular.
+    precip_models: array-like
+      Array of shape (n_models,timesteps+1) containing, per timestep (t=0 to
+      lead time here) and per (NWP) model or model ensemble member, a
+      dictionary with a list of cascades obtained by calling a method
+      implemented in :py:mod:`pysteps.cascade.decomposition`. In case of one
+      (deterministic) model as input, add an extra dimension to make sure
+      precip_models is five dimensional prior to calling this function.
+      It is also possible to supply the original (NWP) model forecasts containing
+      rainfall fields as an array of shape (n_models,timestep+1,m,n), which will
+      then be decomposed in this function. Note that for an operational application
+      or for testing with multiple model runs, it is recommended to decompose
+      the model forecasts outside beforehand, as this reduces calculation times.
+      This is possible with :py:func:`pysteps.blending.utils.decompose_NWP`,
+      :py:func:`pysteps.blending.utils.compute_store_nwp_motion`, and
+      :py:func:`pysteps.blending.utils.load_NWP`.
+    velocity: array-like
+      Array of shape (2,m,n) containing the x- and y-components of the advection
+      field. The velocities are assumed to represent one time step between the
+      inputs. All values are required to be finite.
+    velocity_models: array-like
+      Array of shape (n_models,timestep,2,m,n) containing the x- and y-components
+      of the advection field for the (NWP) model field per forecast lead time.
+      All values are required to be finite.
+    timesteps: int or list of floats
+      Number of time steps to forecast or a list of time steps for which the
+      forecasts are computed (relative to the input time step). The elements of
+      the list are required to be in ascending order.
+    timestep: float
+      Time step of the motion vectors (minutes). Required if vel_pert_method is
+      not None or mask_method is 'incremental'.
+    issuetime: datetime
+      Datetime object containing the date and time for which the forecast
+      is issued.
+    n_ens_members: int
+      The number of ensemble members to generate.
+    n_cascade_levels: int, optional
+      The number of cascade levels to use. Default set to 8 due to default
+      climatological skill values on 8 levels.
+    blend_nwp_members: bool
+      Check if NWP models/members should be used individually, or if all of
+      them are blended together per nowcast ensemble member. Standard set to
+      false.
+    precip_thr: float, optional
+      Specifies the threshold value for minimum observable precipitation
+      intensity. Required if mask_method is not None or conditional is True.
+    kmperpixel: float, optional
+      Spatial resolution of the input data (kilometers/pixel). Required if
+      vel_pert_method is not None or mask_method is 'incremental'.
+    extrap_method: str, optional
+      Name of the extrapolation method to use. See the documentation of
+      :py:mod:`pysteps.extrapolation.interface`.
+    decomp_method: {'fft'}, optional
+      Name of the cascade decomposition method to use. See the documentation
+      of :py:mod:`pysteps.cascade.interface`.
+    bandpass_filter_method: {'gaussian', 'uniform'}, optional
+      Name of the bandpass filter method to use with the cascade decomposition.
+      See the documentation of :py:mod:`pysteps.cascade.interface`.
+    noise_method: {'parametric','nonparametric','ssft','nested',None}, optional
+      Name of the noise generator to use for perturbating the precipitation
+      field. See the documentation of :py:mod:`pysteps.noise.interface`. If set to None,
+      no noise is generated.
+    noise_stddev_adj: {'auto','fixed',None}, optional
+      Optional adjustment for the standard deviations of the noise fields added
+      to each cascade level. This is done to compensate incorrect std. dev.
+      estimates of casace levels due to presence of no-rain areas. 'auto'=use
+      the method implemented in :py:func:`pysteps.noise.utils.compute_noise_stddev_adjs`.
+      'fixed'= use the formula given in :cite:`BPS2006` (eq. 6), None=disable
+      noise std. dev adjustment.
+    ar_order: int, optional
+      The order of the autoregressive model to use. Must be >= 1.
+    vel_pert_method: {'bps',None}, optional
+      Name of the noise generator to use for perturbing the advection field. See
+      the documentation of :py:mod:`pysteps.noise.interface`. If set to None, the advection
+      field is not perturbed.
+    weights_method: {'bps','spn'}, optional
+      The calculation method of the blending weights. Options are the method
+      by :cite:`BPS2006` and the covariance-based method by :cite:`SPN2013`.
+      Defaults to bps.
+    conditional: bool, optional
+      If set to True, compute the statistics of the precipitation field
+      conditionally by excluding pixels where the values are below the threshold
+      precip_thr.
+    mask_method: {'obs','incremental',None}, optional
+      The method to use for masking no precipitation areas in the forecast field.
+      The masked pixels are set to the minimum value of the observations.
+      'obs' = apply precip_thr to the most recently observed precipitation intensity
+      field, 'incremental' = iteratively buffer the mask with a certain rate
+      (currently it is 1 km/min), None=no masking.
+    probmatching_method: {'cdf','mean',None}, optional
+      Method for matching the statistics of the forecast field with those of
+      the most recently observed one. 'cdf'=map the forecast CDF to the observed
+      one, 'mean'=adjust only the conditional mean value of the forecast field
+      in precipitation areas, None=no matching applied. Using 'mean' requires
+      that mask_method is not None.
+    callback: function, optional
+      Optional function that is called after computation of each time step of
+      the nowcast. The function takes one argument: a three-dimensional array
+      of shape (n_ens_members,h,w), where h and w are the height and width
+      of the input field precip, respectively. This can be used, for instance,
+      writing the outputs into files.
+    return_output: bool, optional
+      Set to False to disable returning the outputs as numpy arrays. This can
+      save memory if the intermediate results are written to output files using
+      the callback function.
+    seed: int, optional
+      Optional seed number for the random generators.
+    num_workers: int, optional
+      The number of workers to use for parallel computation. Applicable if dask
+      is enabled or pyFFTW is used for computing the FFT. When num_workers>1, it
+      is advisable to disable OpenMP by setting the environment variable
+      OMP_NUM_THREADS to 1. This avoids slowdown caused by too many simultaneous
+      threads.
+    fft_method: str, optional
+      A string defining the FFT method to use (see FFT methods in
+      :py:func:`pysteps.utils.interface.get_method`).
+      Defaults to 'numpy' for compatibility reasons. If pyFFTW is installed,
+      the recommended method is 'pyfftw'.
+    domain: {"spatial", "spectral"}
+      If "spatial", all computations are done in the spatial domain (the
+      classical STEPS model). If "spectral", the AR(2) models and stochastic
+      perturbations are applied directly in the spectral domain to reduce
+      memory footprint and improve performance :cite:`PCH2019b`.
+    outdir_path_skill: string, optional
+      Path to folder where the historical skill are stored. Defaults to
+      path_workdir from rcparams. If no path is given, './tmp' will be used.
+    extrap_kwargs: dict, optional
+      Optional dictionary containing keyword arguments for the extrapolation
+      method. See the documentation of :py:func:`pysteps.extrapolation.interface`.
+    filter_kwargs: dict, optional
+      Optional dictionary containing keyword arguments for the filter method.
+      See the documentation of :py:mod:`pysteps.cascade.bandpass_filters`.
+    noise_kwargs: dict, optional
+      Optional dictionary containing keyword arguments for the initializer of
+      the noise generator. See the documentation of :py:mod:`pysteps.noise.fftgenerators`.
+    vel_pert_kwargs: dict, optional
+      Optional dictionary containing keyword arguments 'p_par' and 'p_perp' for
+      the initializer of the velocity perturbator. The choice of the optimal
+      parameters depends on the domain and the used optical flow method.
+
+      Default parameters from :cite:`BPS2006`:
+      p_par  = [10.88, 0.23, -7.68]
+      p_perp = [5.76, 0.31, -2.72]
+
+      Parameters fitted to the data (optical flow/domain):
+
+      darts/fmi:
+      p_par  = [13.71259667, 0.15658963, -16.24368207]
+      p_perp = [8.26550355, 0.17820458, -9.54107834]
+
+      darts/mch:
+      p_par  = [24.27562298, 0.11297186, -27.30087471]
+      p_perp = [-7.80797846e+01, -3.38641048e-02, 7.56715304e+01]
+
+      darts/fmi+mch:
+      p_par  = [16.55447057, 0.14160448, -19.24613059]
+      p_perp = [14.75343395, 0.11785398, -16.26151612]
+
+      lucaskanade/fmi:
+      p_par  = [2.20837526, 0.33887032, -2.48995355]
+      p_perp = [2.21722634, 0.32359621, -2.57402761]
+
+      lucaskanade/mch:
+      p_par  = [2.56338484, 0.3330941, -2.99714349]
+      p_perp = [1.31204508, 0.3578426, -1.02499891]
+
+      lucaskanade/fmi+mch:
+      p_par  = [2.31970635, 0.33734287, -2.64972861]
+      p_perp = [1.90769947, 0.33446594, -2.06603662]
+
+      vet/fmi:
+      p_par  = [0.25337388, 0.67542291, 11.04895538]
+      p_perp = [0.02432118, 0.99613295, 7.40146505]
+
+      vet/mch:
+      p_par  = [0.5075159, 0.53895212, 7.90331791]
+      p_perp = [0.68025501, 0.41761289, 4.73793581]
+
+      vet/fmi+mch:
+      p_par  = [0.29495222, 0.62429207, 8.6804131 ]
+      p_perp = [0.23127377, 0.59010281, 5.98180004]
+
+      fmi=Finland, mch=Switzerland, fmi+mch=both pooled into the same data set
+
+      The above parameters have been fitten by using run_vel_pert_analysis.py
+      and fit_vel_pert_params.py located in the scripts directory.
+
+      See :py:mod:`pysteps.noise.motion` for additional documentation.
+    clim_kwargs: dict, optional
+      Optional dictionary containing keyword arguments for the climatological
+      skill file. Arguments can consist of: 'outdir_path', 'n_models'
+      (the number of NWP models) and 'window_length' (the minimum number of
+      days the clim file should have, otherwise the default is used).
+    mask_kwargs: dict
+      Optional dictionary containing mask keyword arguments 'mask_f' and
+      'mask_rim', the factor defining the the mask increment and the rim size,
+      respectively.
+      The mask increment is defined as mask_f*timestep/kmperpixel.
+    measure_time: bool
+      If set to True, measure, print and return the computation time.
+
+    Returns
+    -------
+    out: ndarray
+      If return_output is True, a four-dimensional array of shape
+      (n_ens_members,num_timesteps,m,n) containing a time series of forecast
+      precipitation fields for each ensemble member. Otherwise, a None value
+      is returned. The time series starts from t0+timestep, where timestep is
+      taken from the input precipitation fields precip. If measure_time is True, the
+      return value is a three-element tuple containing the nowcast array, the
+      initialization time of the nowcast generator and the time used in the
+      main loop (seconds).
+
+    See also
+    --------
+    :py:mod:`pysteps.extrapolation.interface`, :py:mod:`pysteps.cascade.interface`,
+    :py:mod:`pysteps.noise.interface`, :py:func:`pysteps.noise.utils.compute_noise_stddev_adjs`
+
+    References
+    ----------
+    :cite:`Seed2003`, :cite:`BPS2004`, :cite:`BPS2006`, :cite:`SPN2013`, :cite:`PCH2019b`
+
+    Notes
+    -----
+    1. The blending currently does not blend the beta-parameters in the parametric
+    noise method. It is recommended to use the non-parameteric noise method.
+
+    2. If blend_nwp_members is True, the BPS2006 method for the weights is
+    suboptimal. It is recommended to use the SPN2013 method instead.
+
+    3. Not yet implemented (and neither in the steps nowcasting module): The regression
+    of the lag-1 and lag-2 parameters to their climatological values. See also eq.
+    12 - 19 in :cite: `BPS2004`. By doing so, the Phi parameters change over time,
+    which enhances the AR process. This can become a future development if this
+    turns out to be a warranted functionality.
+    """
+
+    # 0.1 Start with some checks
+    _check_inputs(precip, precip_models, velocity, velocity_models, timesteps, ar_order)
+
+    if extrap_kwargs is None:
+        extrap_kwargs = dict()
+
+    if filter_kwargs is None:
+        filter_kwargs = dict()
+
+    if noise_kwargs is None:
+        noise_kwargs = dict()
+
+    if vel_pert_kwargs is None:
+        vel_pert_kwargs = dict()
+
+    if clim_kwargs is None:
+        clim_kwargs = dict()
+
+    if mask_kwargs is None:
+        mask_kwargs = dict()
+
+    if np.any(~np.isfinite(velocity)):
+        raise ValueError("velocity contains non-finite values")
+
+    if mask_method not in ["obs", "incremental", None]:
+        raise ValueError(
+            "unknown mask method %s: must be 'obs', 'incremental' or None" % mask_method
+        )
+
+    if conditional and precip_thr is None:
+        raise ValueError("conditional=True but precip_thr is not set")
+
+    if mask_method is not None and precip_thr is None:
+        raise ValueError("mask_method!=None but precip_thr=None")
+
+    if noise_stddev_adj not in ["auto", "fixed", None]:
+        raise ValueError(
+            "unknown noise_std_dev_adj method %s: must be 'auto', 'fixed', or None"
+            % noise_stddev_adj
+        )
+
+    if kmperpixel is None:
+        if vel_pert_method is not None:
+            raise ValueError("vel_pert_method is set but kmperpixel=None")
+        if mask_method == "incremental":
+            raise ValueError("mask_method='incremental' but kmperpixel=None")
+
+    if timestep is None:
+        if vel_pert_method is not None:
+            raise ValueError("vel_pert_method is set but timestep=None")
+        if mask_method == "incremental":
+            raise ValueError("mask_method='incremental' but timestep=None")
+
+    # 0.2 Log some settings
+    print("STEPS blending")
+    print("==============")
+    print("")
+
+    print("Inputs")
+    print("------")
+    print(f"forecast issue time:         {issuetime.isoformat()}")
+    print(f"input dimensions:            {precip.shape[1]}x{precip.shape[2]}")
+    if kmperpixel is not None:
+        print(f"km/pixel:                    {kmperpixel}")
+    if timestep is not None:
+        print(f"time step:                   {timestep} minutes")
+    print("")
+
+    print("NWP and blending inputs")
+    print("-----------------------")
+    print(f"number of (NWP) models:      {precip_models.shape[0]}")
+    print(f"blend (NWP) model members:   {blend_nwp_members}")
+    print(f"decompose (NWP) models:      {'yes' if precip_models.ndim == 4 else 'no'}")
+    print("")
+
+    print("Methods")
+    print("-------")
+    print(f"extrapolation:               {extrap_method}")
+    print(f"bandpass filter:             {bandpass_filter_method}")
+    print(f"decomposition:               {decomp_method}")
+    print(f"noise generator:             {noise_method}")
+    print(f"noise adjustment:            {'yes' if noise_stddev_adj else 'no'}")
+    print(f"velocity perturbator:        {vel_pert_method}")
+    print(f"blending weights method:     {weights_method}")
+    print(f"conditional statistics:      {'yes' if conditional else 'no'}")
+    print(f"precip. mask method:         {mask_method}")
+    print(f"probability matching:        {probmatching_method}")
+    print(f"FFT method:                  {fft_method}")
+    print(f"domain:                      {domain}")
+    print("")
+
+    print("Parameters")
+    print("----------")
+    if isinstance(timesteps, int):
+        print(f"number of time steps:        {timesteps}")
+    else:
+        print(f"time steps:                  {timesteps}")
+    print(f"ensemble size:               {n_ens_members}")
+    print(f"parallel threads:            {num_workers}")
+    print(f"number of cascade levels:    {n_cascade_levels}")
+    print(f"order of the AR(p) model:    {ar_order}")
+    if vel_pert_method == "bps":
+        vp_par = vel_pert_kwargs.get("p_par", noise.motion.get_default_params_bps_par())
+        vp_perp = vel_pert_kwargs.get(
+            "p_perp", noise.motion.get_default_params_bps_perp()
+        )
+        print(f"vel. pert., parallel:        {vp_par[0]},{vp_par[1]},{vp_par[2]}")
+        print(f"vel. pert., perpendicular:   {vp_perp[0]},{vp_perp[1]},{vp_perp[2]}")
+    else:
+        vp_par, vp_perp = None, None
+
+    if conditional or mask_method is not None:
+        print(f"precip. intensity threshold: {precip_thr}")
+    print("")
+
+    # 0.3 Get the methods that will be used
+    num_ensemble_workers = n_ens_members if num_workers > n_ens_members else num_workers
+
+    if measure_time:
+        starttime_init = time.time()
+
+    fft = utils.get_method(fft_method, shape=precip.shape[1:], n_threads=num_workers)
+
+    precip_shape = precip.shape[1:]
+
+    # initialize the band-pass filter
+    filter_method = cascade.get_method(bandpass_filter_method)
+    bp_filter = filter_method(precip_shape, n_cascade_levels, **filter_kwargs)
+
+    decompositor, recompositor = cascade.get_method(decomp_method)
+
+    extrapolator = extrapolation.get_method(extrap_method)
+
+    x_values, y_values = np.meshgrid(
+        np.arange(precip.shape[2]), np.arange(precip.shape[1])
+    )
+
+    xy_coords = np.stack([x_values, y_values])
+
+    precip = precip[-(ar_order + 1) :, :, :].copy()
+
+    # determine the domain mask from non-finite values
+    domain_mask = np.logical_or.reduce(
+        [~np.isfinite(precip[i, :]) for i in range(precip.shape[0])]
+    )
+
+    # determine the precipitation threshold mask
+    if conditional:
+        MASK_thr = np.logical_and.reduce(
+            [precip[i, :, :] >= precip_thr for i in range(precip.shape[0])]
+        )
+    else:
+        MASK_thr = None
+
+    # 1. Start with the radar rainfall fields. We want the fields in a
+    # Lagrangian space
+    precip = _transform_to_lagrangian(
+        precip, velocity, ar_order, xy_coords, extrapolator, extrap_kwargs, num_workers
+    )
+
+    # 2. Initialize the noise method
+    pp, generate_noise, noise_std_coeffs = _init_noise(
+        precip,
+        precip_thr,
+        n_cascade_levels,
+        bp_filter,
+        decompositor,
+        fft,
+        noise_method,
+        noise_kwargs,
+        noise_stddev_adj,
+        measure_time,
+        num_workers,
+    )
+
+    # 3. Perform the cascade decomposition for the input precip fields and
+    # and, if necessary, for the (NWP) model fields
+
+    # 3.1 Compute the cascade decompositions of the input precipitation fields
+    (
+        precip_cascade,
+        mu_extrapolation,
+        sigma_extrapolation,
+    ) = _compute_cascade_decomposition_radar(
+        precip,
+        ar_order,
+        n_cascade_levels,
+        n_ens_members,
+        MASK_thr,
+        domain,
+        bp_filter,
+        decompositor,
+        fft,
+    )
+
+    # 3.2 If necessary, decompose (NWP) model forecasts and stack cascades
+    (
+        precip_models_cascade,
+        mu_models,
+        sigma_models,
+        precip_models_pm,
+    ) = _compute_cascade_decomposition_nwp(
+        precip_models, bp_filter, decompositor, recompositor, fft, domain
+    )
+
+    # 4. Estimate AR parameters for the radar rainfall field
+    PHI = _estimate_ar_parameters_radar(
+        precip_cascade, ar_order, n_cascade_levels, MASK_thr
+    )
+
+    # 5. Before calling the worker for the forecast loop, determine which (NWP)
+    # models will be combined with which nowcast ensemble members. With the
+    # way it is implemented at this moment: n_ens_members of the output equals
+    # the maximum number of (ensemble) members in the input (either the nowcasts or NWP).
+    (
+        precip_cascade,
+        precip_models_cascade,
+        precip_models_pm,
+        velocity_models,
+        mu_models,
+        sigma_models,
+        n_ens_members,
+        n_model_indices,
+    ) = _find_nwp_combination(
+        precip_cascade,
+        precip_models_cascade,
+        precip_models_pm,
+        velocity_models,
+        mu_models,
+        sigma_models,
+        n_ens_members,
+        ar_order,
+        n_cascade_levels,
+        blend_nwp_members,
+    )
+
+    # Also initialize the cascade of temporally correlated noise, which has the
+    # same shape as R_c, but starts with value zero.
+    noise_cascade = np.zeros(precip_cascade.shape)
+
+    # 6. Initialize all the random generators and prepare for the forecast loop
+    randgen_prec, vps, generate_vel_noise = _init_random_generators(
+        velocity,
+        noise_method,
+        vel_pert_method,
+        vp_par,
+        vp_perp,
+        seed,
+        n_ens_members,
+        kmperpixel,
+        timestep,
+    )
+    D, D_Yn, D_pb, R_f, R_m, mask_rim, struct, fft_objs = _prepare_forecast_loop(
+        precip_cascade,
+        noise_method,
+        fft_method,
+        n_cascade_levels,
+        n_ens_members,
+        mask_method,
+        mask_kwargs,
+        timestep,
+        kmperpixel,
+    )
+
+    precip = precip[-1, :, :]
+
+    # 7. Calculate the initial skill of the (NWP) model forecasts at t=0
+    rho_nwp_models = _compute_initial_nwp_skill(
+        precip_cascade,
+        precip_models_cascade,
+        domain_mask,
+        issuetime,
+        outdir_path_skill,
+        clim_kwargs,
+    )
+
+    # Also initizalize the current and previous extrapolation forecast scale
+    # for the nowcasting component
+    rho_extr_prev = np.repeat(1.0, PHI.shape[0])
+    rho_extr = PHI[:, 0] / (1.0 - PHI[:, 1])  # phi1 / (1 - phi2), see BPS2004
+
+    if measure_time:
+        init_time = time.time() - starttime_init
+
+    ###
+    # 8. Start the forecasting loop
+    ###
+    print("Starting blended nowcast computation.")
+
+    if measure_time:
+        starttime_mainloop = time.time()
+
+    if isinstance(timesteps, int):
+        timesteps = range(timesteps + 1)
+        timestep_type = "int"
+    else:
+        original_timesteps = [0] + list(timesteps)
+        timesteps = nowcast_utils.binned_timesteps(original_timesteps)
+        timestep_type = "list"
+
+    extrap_kwargs["return_displacement"] = True
+    forecast_prev = precip_cascade
+    noise_prev = noise_cascade
+    t_prev = [0.0 for j in range(n_ens_members)]
+    t_total = [0.0 for j in range(n_ens_members)]
+
+    # iterate each time step
+    for t, subtimestep_idx in enumerate(timesteps):
+        if timestep_type == "list":
+            subtimesteps = [original_timesteps[t_] for t_ in subtimestep_idx]
+        else:
+            subtimesteps = [t]
+
+        if (timestep_type == "list" and subtimesteps) or (
+            timestep_type == "int" and t > 0
+        ):
+            is_nowcast_time_step = True
+        else:
+            is_nowcast_time_step = False
+
+        if is_nowcast_time_step:
+            print(
+                f"Computing nowcast for time step {t}... ",
+                end="",
+                flush=True,
+            )
+
+        if measure_time:
+            starttime = time.time()
+
+        # 8.1.1 Determine the skill of the components for lead time (t0 + t)
+        # First for the extrapolation component. Only calculate it when t > 0.
+        if t > 0:
+            (
+                rho_extr,
+                rho_extr_prev,
+            ) = blending.skill_scores.lt_dependent_cor_extrapolation(
+                PHI=PHI, correlations=rho_extr, correlations_prev=rho_extr_prev
+            )
+
+        # the nowcast iteration for each ensemble member
+        def worker(j):
+            # 8.1.2 Determine the skill of the nwp components for lead time (t0 + t)
+            # Then for the model components
+            if blend_nwp_members:
+                rho_nwp_fc = [
+                    blending.skill_scores.lt_dependent_cor_nwp(
+                        lt=(t * int(timestep)),
+                        correlations=rho_nwp_models[n_model],
+                        outdir_path=outdir_path_skill,
+                        n_model=n_model,
+                        skill_kwargs=clim_kwargs,
+                    )
+                    for n_model in range(rho_nwp_models.shape[0])
+                ]
+                rho_nwp_fc = np.stack(rho_nwp_fc)
+                # Concatenate rho_extr and rho_nwp
+                rho_fc = np.concatenate((rho_extr[None, :], rho_nwp_fc), axis=0)
+            else:
+                rho_nwp_fc = blending.skill_scores.lt_dependent_cor_nwp(
+                    lt=(t * int(timestep)),
+                    correlations=rho_nwp_models[j],
+                    outdir_path=outdir_path_skill,
+                    n_model=n_model_indices[j],
+                    skill_kwargs=clim_kwargs,
+                )
+                # Concatenate rho_extr and rho_nwp
+                rho_fc = np.concatenate(
+                    (rho_extr[None, :], rho_nwp_fc[None, :]), axis=0
+                )
+
+            # 8.2 Determine the weights per component
+
+            # Weights following the bps method. These are needed for the velocity
+            # weights prior to the advection step. If weights method spn is
+            # selected, weights will be overwritten with those weights prior to
+            # blending step.
+            # weight = [(extr_field, n_model_fields, noise), n_cascade_levels, ...]
+            weights = calculate_weights_bps(rho_fc)
+
+            # The model only weights
+            if weights_method == "bps":
+                # Determine the weights of the components without the extrapolation
+                # cascade, in case this is no data or outside the mask.
+                weights_model_only = calculate_weights_bps(rho_fc[1:, :])
+            elif weights_method == "spn":
+                # Only the weights of the components without the extrapolation
+                # cascade will be determined here. The full set of weights are
+                # determined after the extrapolation step in this method.
+                if blend_nwp_members and precip_models_cascade.shape[0] > 1:
+                    weights_model_only = np.zeros(
+                        (precip_models_cascade.shape[0] + 1, n_cascade_levels)
+                    )
+                    for i in range(n_cascade_levels):
+                        # Determine the normalized covariance matrix (containing)
+                        # the cross-correlations between the models
+                        cov = np.corrcoef(
+                            np.stack(
+                                [
+                                    precip_models_cascade[n_model, t, i, :, :].flatten()
+                                    for n_model in range(precip_models_cascade.shape[0])
+                                ]
+                            )
+                        )
+                        # Determine the weights for this cascade level
+                        weights_model_only[:, i] = calculate_weights_spn(
+                            correlations=rho_fc[1:, i], cov=cov
+                        )
+                else:
+                    # Same as correlation and noise is 1 - correlation
+                    weights_model_only = calculate_weights_bps(rho_fc[1:, :])
+            else:
+                raise ValueError(
+                    "Unknown weights method %s: must be 'bps' or 'spn'" % weights_method
+                )
+
+            # 8.3 Determine the noise cascade and regress this to the subsequent
+            # time step + regress the extrapolation component to the subsequent
+            # time step
+
+            # 8.3.1 Determine the epsilon, a cascade of temporally independent
+            # but spatially correlated noise
+            if noise_method is not None:
+                # generate noise field
+                EPS = generate_noise(
+                    pp, randstate=randgen_prec[j], fft_method=fft_objs[j], domain=domain
+                )
+
+                # decompose the noise field into a cascade
+                EPS = decompositor(
+                    EPS,
+                    bp_filter,
+                    fft_method=fft_objs[j],
+                    input_domain=domain,
+                    output_domain=domain,
+                    compute_stats=True,
+                    normalize=True,
+                    compact_output=True,
+                )
+            else:
+                EPS = None
+
+            # 8.3.2 regress the extrapolation component to the subsequent time
+            # step
+            # iterate the AR(p) model for each cascade level
+            for i in range(n_cascade_levels):
+                # apply AR(p) process to extrapolation cascade level
+                if EPS is not None or vel_pert_method is not None:
+                    precip_cascade[j][i] = autoregression.iterate_ar_model(
+                        precip_cascade[j][i], PHI[i, :]
+                    )
+
+                else:
+                    # use the deterministic AR(p) model computed above if
+                    # perturbations are disabled
+                    precip_cascade[j][i] = R_m[i]
+
+            # 8.3.3 regress the noise component to the subsequent time step
+            # iterate the AR(p) model for each cascade level
+            for i in range(n_cascade_levels):
+                # normalize the noise cascade
+                if EPS is not None:
+                    EPS_ = EPS["cascade_levels"][i]
+                    EPS_ *= noise_std_coeffs[i]
+                else:
+                    EPS_ = None
+                # apply AR(p) process to noise cascade level
+                # (Returns zero noise if EPS is None)
+                noise_cascade[j][i] = autoregression.iterate_ar_model(
+                    noise_cascade[j][i], PHI[i, :], eps=EPS_
+                )
+
+            EPS = None
+            EPS_ = None
+
+            # 8.4 Perturb and blend the advection fields + advect the
+            # extrapolation and noise cascade to the current time step
+            # (or subtimesteps if non-integer time steps are given)
+
+            # Settings and initialize the output
+            extrap_kwargs_ = extrap_kwargs.copy()
+            extrap_kwargs_noise = extrap_kwargs.copy()
+            extrap_kwargs_pb = extrap_kwargs.copy()
+            velocity_pert = velocity
+            R_f_ep_out = []
+            Yn_ep_out = []
+            R_pm_ep = []
+
+            # Extrapolate per sub time step
+            for t_sub in subtimesteps:
+                if t_sub > 0:
+                    t_diff_prev_int = t_sub - int(t_sub)
+                    if t_diff_prev_int > 0.0:
+                        R_f_ip = [
+                            (1.0 - t_diff_prev_int) * forecast_prev[j][i][-1, :]
+                            + t_diff_prev_int * precip_cascade[j][i][-1, :]
+                            for i in range(n_cascade_levels)
+                        ]
+                        Yn_ip = [
+                            (1.0 - t_diff_prev_int) * noise_prev[j][i][-1, :]
+                            + t_diff_prev_int * noise_cascade[j][i][-1, :]
+                            for i in range(n_cascade_levels)
+                        ]
+
+                    else:
+                        R_f_ip = [
+                            forecast_prev[j][i][-1, :] for i in range(n_cascade_levels)
+                        ]
+                        Yn_ip = [
+                            noise_prev[j][i][-1, :] for i in range(n_cascade_levels)
+                        ]
+
+                    R_f_ip = np.stack(R_f_ip)
+                    Yn_ip = np.stack(Yn_ip)
+
+                    t_diff_prev = t_sub - t_prev[j]
+                    t_total[j] += t_diff_prev
+
+                    # compute the perturbed motion field - include the NWP
+                    # velocities and the weights. Note that we only perturb
+                    # the extrapolation velocity field, as the NWP velocity
+                    # field is present per time step
+                    if vel_pert_method is not None:
+                        velocity_pert = velocity + generate_vel_noise(
+                            vps[j], t_total[j] * timestep
+                        )
+
+                    # Stack the perturbed extrapolation and the NWP velocities
+                    if blend_nwp_members:
+                        V_stack = np.concatenate(
+                            (
+                                velocity_pert[None, :, :, :],
+                                velocity_models[:, t, :, :, :],
+                            ),
+                            axis=0,
+                        )
+                    else:
+                        V_model_ = velocity_models[j, t, :, :, :]
+                        V_stack = np.concatenate(
+                            (velocity_pert[None, :, :, :], V_model_[None, :, :, :]),
+                            axis=0,
+                        )
+                        V_model_ = None
+
+                    # Obtain a blended optical flow, using the weights of the
+                    # second cascade following eq. 24 in BPS2006
+                    velocity_blended = blending.utils.blend_optical_flows(
+                        flows=V_stack,
+                        weights=weights[
+                            :-1, 1
+                        ],  # [(extr_field, n_model_fields), cascade_level=2]
+                    )
+
+                    # Extrapolate both cascades to the next time step
+                    R_f_ep = np.zeros(R_f_ip.shape)
+                    Yn_ep = np.zeros(Yn_ip.shape)
+
+                    for i in range(n_cascade_levels):
+                        extrap_kwargs_["displacement_prev"] = D[j][i]
+                        extrap_kwargs_noise["displacement_prev"] = D_Yn[j][i]
+                        extrap_kwargs_noise["map_coordinates_mode"] = "wrap"
+
+                        # First, extrapolate the extrapolation component
+                        # Apply the domain mask to the extrapolation component
+                        R_f_ip[i][domain_mask] = np.NaN
+                        R_f_ep_, D[j][i] = extrapolator(
+                            R_f_ip[i],
+                            velocity_blended,
+                            [t_diff_prev],
+                            allow_nonfinite_values=True,
+                            **extrap_kwargs_,
+                        )
+                        R_f_ep[i] = R_f_ep_[0]
+
+                        # Then, extrapolate the noise component
+                        Yn_ep_, D_Yn[j][i] = extrapolator(
+                            Yn_ip[i],
+                            velocity_blended,
+                            [t_diff_prev],
+                            allow_nonfinite_values=True,
+                            **extrap_kwargs_noise,
+                        )
+                        Yn_ep[i] = Yn_ep_[0]
+
+                    # Append the results to the output lists
+                    R_f_ep_out.append(R_f_ep)
+                    Yn_ep_out.append(Yn_ep)
+                    R_f_ep_ = None
+                    Yn_ep_ = None
+
+                    # Finally, also extrapolate the initial radar rainfall
+                    # field. This will be blended with the rainfall field(s)
+                    # of the (NWP) model(s) for Lagrangian blended prob. matching
+                    # min_R = np.min(precip)
+                    extrap_kwargs_pb["displacement_prev"] = D_pb[j]
+                    # Apply the domain mask to the extrapolation component
+                    R_ = precip.copy()
+                    R_[domain_mask] = np.NaN
+                    R_pm_ep_, D_pb[j] = extrapolator(
+                        R_,
+                        velocity_blended,
+                        [t_diff_prev],
+                        allow_nonfinite_values=True,
+                        **extrap_kwargs_pb,
+                    )
+                    R_pm_ep.append(R_pm_ep_[0])
+
+                    t_prev[j] = t_sub
+
+                forecast_prev[j] = precip_cascade[j]
+
+            if len(R_f_ep_out) > 0:
+                R_f_ep_out = np.stack(R_f_ep_out)
+                Yn_ep_out = np.stack(Yn_ep_out)
+                R_pm_ep = np.stack(R_pm_ep)
+
+            # advect the forecast field by one time step if no subtimesteps in the
+            # current interval were found
+            if not subtimesteps:
+                t_diff_prev = t + 1 - t_prev[j]
+                t_total[j] += t_diff_prev
+
+                # compute the perturbed motion field - include the NWP
+                # velocities and the weights
+                if vel_pert_method is not None:
+                    velocity_pert = velocity + generate_vel_noise(
+                        vps[j], t_total[j] * timestep
+                    )
+
+                # Stack the perturbed extrapolation and the NWP velocities
+                if blend_nwp_members:
+                    V_stack = np.concatenate(
+                        (velocity_pert[None, :, :, :], velocity_models[:, t, :, :, :]),
+                        axis=0,
+                    )
+                else:
+                    V_model_ = velocity_models[j, t, :, :, :]
+                    V_stack = np.concatenate(
+                        (velocity_pert[None, :, :, :], V_model_[None, :, :, :]), axis=0
+                    )
+                    V_model_ = None
+
+                # Obtain a blended optical flow, using the weights of the
+                # second cascade following eq. 24 in BPS2006
+                velocity_blended = blending.utils.blend_optical_flows(
+                    flows=V_stack,
+                    weights=weights[
+                        :-1, 1
+                    ],  # [(extr_field, n_model_fields), cascade_level=2]
+                )
+
+                # Extrapolate the extrapolation and noise cascade
+                for i in range(n_cascade_levels):
+                    extrap_kwargs_["displacement_prev"] = D[j][i]
+                    extrap_kwargs_noise["displacement_prev"] = D_Yn[j][i]
+                    extrap_kwargs_noise["map_coordinates_mode"] = "wrap"
+
+                    _, D[j][i] = extrapolator(
+                        None,
+                        velocity_blended,
+                        [t_diff_prev],
+                        allow_nonfinite_values=True,
+                        **extrap_kwargs_,
+                    )
+
+                    _, D_Yn[j][i] = extrapolator(
+                        None,
+                        velocity_blended,
+                        [t_diff_prev],
+                        allow_nonfinite_values=True,
+                        **extrap_kwargs_noise,
+                    )
+
+                # Also extrapolate the radar observation, used for the probability
+                # matching and post-processing steps
+                extrap_kwargs_pb["displacement_prev"] = D_pb[j]
+                _, D_pb[j] = extrapolator(
+                    None,
+                    velocity_blended,
+                    [t_diff_prev],
+                    allow_nonfinite_values=True,
+                    **extrap_kwargs_pb,
+                )
+
+                t_prev[j] = t + 1
+
+            # 8.5 Blend the cascades
+            R_f_out = []
+
+            for t_sub in subtimesteps:
+                # TODO: does it make sense to use sub time steps - check if it works?
+                if t_sub > 0:
+                    t_index = np.where(np.array(subtimesteps) == t_sub)[0][0]
+                    # First concatenate the cascades and the means and sigmas
+                    # precip_models = [n_models,timesteps,n_cascade_levels,m,n]
+                    if blend_nwp_members:
+                        cascades_stacked = np.concatenate(
+                            (
+                                R_f_ep_out[None, t_index],
+                                precip_models_cascade[:, t],
+                                Yn_ep_out[None, t_index],
+                            ),
+                            axis=0,
+                        )  # [(extr_field, n_model_fields, noise), n_cascade_levels, ...]
+                        means_stacked = np.concatenate(
+                            (mu_extrapolation[None, :], mu_models[:, t]), axis=0
+                        )
+                        sigmas_stacked = np.concatenate(
+                            (sigma_extrapolation[None, :], sigma_models[:, t]),
+                            axis=0,
+                        )
+                    else:
+                        cascades_stacked = np.concatenate(
+                            (
+                                R_f_ep_out[None, t_index],
+                                precip_models_cascade[None, j, t],
+                                Yn_ep_out[None, t_index],
+                            ),
+                            axis=0,
+                        )  # [(extr_field, n_model_fields, noise), n_cascade_levels, ...]
+                        means_stacked = np.concatenate(
+                            (mu_extrapolation[None, :], mu_models[None, j, t]), axis=0
+                        )
+                        sigmas_stacked = np.concatenate(
+                            (sigma_extrapolation[None, :], sigma_models[None, j, t]),
+                            axis=0,
+                        )
+
+                    # First determine the blending weights if method is spn. The
+                    # weights for method bps have already been determined.
+                    if weights_method == "spn":
+                        weights = np.zeros(
+                            (cascades_stacked.shape[0], n_cascade_levels)
+                        )
+                        for i in range(n_cascade_levels):
+                            # Determine the normalized covariance matrix (containing)
+                            # the cross-correlations between the models
+                            cascades_stacked_ = np.stack(
+                                [
+                                    cascades_stacked[n_model, i, :, :].flatten()
+                                    for n_model in range(cascades_stacked.shape[0] - 1)
+                                ]
+                            )  # -1 to exclude the noise component
+                            cov = np.ma.corrcoef(
+                                np.ma.masked_invalid(cascades_stacked_)
+                            )
+                            # Determine the weights for this cascade level
+                            weights[:, i] = calculate_weights_spn(
+                                correlations=rho_fc[:, i], cov=cov
+                            )
+
+                    # Blend the extrapolation, (NWP) model(s) and noise cascades
+                    R_f_blended = blending.utils.blend_cascades(
+                        cascades_norm=cascades_stacked, weights=weights
+                    )
+
+                    # Also blend the cascade without the extrapolation component
+                    R_f_blended_mod_only = blending.utils.blend_cascades(
+                        cascades_norm=cascades_stacked[1:, :],
+                        weights=weights_model_only,
+                    )
+
+                    # Blend the means and standard deviations
+                    # Input is array of shape [number_components, scale_level, ...]
+                    means_blended, sigmas_blended = blend_means_sigmas(
+                        means=means_stacked, sigmas=sigmas_stacked, weights=weights
+                    )
+                    # Also blend the means and sigmas for the cascade without extrapolation
+                    (
+                        means_blended_mod_only,
+                        sigmas_blended_mod_only,
+                    ) = blend_means_sigmas(
+                        means=means_stacked[1:, :],
+                        sigmas=sigmas_stacked[1:, :],
+                        weights=weights_model_only,
+                    )
+
+                    # 8.6 Recompose the cascade to a precipitation field
+                    # (The function first normalizes the blended cascade, R_f_blended
+                    # again)
+                    R_f_new = blending.utils.recompose_cascade(
+                        combined_cascade=R_f_blended,
+                        combined_mean=means_blended,
+                        combined_sigma=sigmas_blended,
+                    )
+                    # The recomposed cascade without the extrapolation (for NaN filling
+                    # outside the radar domain)
+                    R_f_new_mod_only = blending.utils.recompose_cascade(
+                        combined_cascade=R_f_blended_mod_only,
+                        combined_mean=means_blended_mod_only,
+                        combined_sigma=sigmas_blended_mod_only,
+                    )
+                    if domain == "spectral":
+                        # TODO: Check this! (Only tested with domain == 'spatial')
+                        R_f_new = fft_objs[j].irfft2(R_f_new)
+                        R_f_new_mod_only = fft_objs[j].irfft2(R_f_new_mod_only)
+
+                    # 8.7 Post-processing steps - use the mask and fill no data with
+                    # the blended NWP forecast. Probability matching following
+                    # Lagrangian blended probability matching which uses the
+                    # latest extrapolated radar rainfall field blended with the
+                    # nwp model(s) rainfall forecast fields as 'benchmark'.
+
+                    # TODO: Check probability matching method
+                    # 8.7.1 first blend the extrapolated rainfall field (the field
+                    # that is only used for post-processing steps) with the NWP
+                    # rainfall forecast for this time step using the weights
+                    # at scale level 2.
+                    weights_pm = weights[:-1, 1]  # Weights without noise, level 2
+                    weights_pm_normalized = weights_pm / np.sum(weights_pm)
+                    # And the weights for outside the radar domain
+                    weights_pm_mod_only = weights_model_only[
+                        :-1, 1
+                    ]  # Weights without noise, level 2
+                    weights_pm_normalized_mod_only = weights_pm_mod_only / np.sum(
+                        weights_pm_mod_only
+                    )
+                    # Stack the fields
+                    if blend_nwp_members:
+                        R_pm_stacked = np.concatenate(
+                            (
+                                R_pm_ep[None, t_index],
+                                precip_models_pm[:, t],
+                            ),
+                            axis=0,
+                        )
+                    else:
+                        R_pm_stacked = np.concatenate(
+                            (
+                                R_pm_ep[None, t_index],
+                                precip_models_pm[None, j, t],
+                            ),
+                            axis=0,
+                        )
+                    # Blend it
+                    R_pm_blended = np.sum(
+                        weights_pm_normalized.reshape(
+                            weights_pm_normalized.shape[0], 1, 1
+                        )
+                        * R_pm_stacked,
+                        axis=0,
+                    )
+                    if blend_nwp_members:
+                        R_pm_blended_mod_only = np.sum(
+                            weights_pm_normalized_mod_only.reshape(
+                                weights_pm_normalized_mod_only.shape[0], 1, 1
+                            )
+                            * precip_models_pm[:, t],
+                            axis=0,
+                        )
+                    else:
+                        R_pm_blended_mod_only = precip_models_pm[j, t]
+
+                    # The extrapolation components are NaN outside the advected
+                    # radar domain. This results in NaN values in the blended
+                    # forecast outside the radar domain. Therefore, fill these
+                    # areas with the "..._mod_only" blended forecasts, consisting
+                    # of the NWP and noise components.
+                    nan_indices = np.isnan(R_f_new)
+                    R_f_new[nan_indices] = R_f_new_mod_only[nan_indices]
+                    nan_indices = np.isnan(R_pm_blended)
+                    R_pm_blended[nan_indices] = R_pm_blended_mod_only[nan_indices]
+                    # Finally, fill the remaining nan values, if present, with
+                    # the minimum value in the forecast
+                    nan_indices = np.isnan(R_f_new)
+                    R_f_new[nan_indices] = np.nanmin(R_f_new)
+                    nan_indices = np.isnan(R_pm_blended)
+                    R_pm_blended[nan_indices] = np.nanmin(R_pm_blended)
+
+                    # 8.7.2. Apply the masking and prob. matching
+                    if mask_method is not None:
+                        # apply the precipitation mask to prevent generation of new
+                        # precipitation into areas where it was not originally
+                        # observed
+                        R_cmin = R_f_new.min()
+                        if mask_method == "incremental":
+                            # The incremental mask is slightly different from
+                            # the implementation in the non-blended steps.py, as
+                            # it is not based on the last forecast, but instead
+                            # on R_pm_blended. Therefore, the buffer does not
+                            # increase over time.
+                            # Get the mask for this forecast
+                            MASK_prec = R_pm_blended >= precip_thr
+                            # Buffer the mask
+                            MASK_prec = _compute_incremental_mask(
+                                MASK_prec, struct, mask_rim
+                            )
+                            # Get the final mask
+                            R_f_new = R_cmin + (R_f_new - R_cmin) * MASK_prec
+                            MASK_prec_ = R_f_new > R_cmin
+                        elif mask_method == "obs":
+                            # The mask equals the most recent benchmark
+                            # rainfall field
+                            MASK_prec_ = R_pm_blended >= precip_thr
+
+                        # Set to min value outside of mask
+                        R_f_new[~MASK_prec_] = R_cmin
+
+                    if probmatching_method == "cdf":
+                        # adjust the CDF of the forecast to match the most recent
+                        # benchmark rainfall field (R_pm_blended)
+                        R_f_new = probmatching.nonparam_match_empirical_cdf(
+                            R_f_new, R_pm_blended
+                        )
+                    elif probmatching_method == "mean":
+                        # Use R_pm_blended as benchmark field and
+                        mu_0 = np.mean(R_pm_blended[R_pm_blended >= precip_thr])
+                        MASK = R_f_new >= precip_thr
+                        mu_fct = np.mean(R_f_new[MASK])
+                        R_f_new[MASK] = R_f_new[MASK] - mu_fct + mu_0
+
+                    R_f_out.append(R_f_new)
+
+            return R_f_out
+
+        res = []
+        for j in range(n_ens_members):
+            if not DASK_IMPORTED or n_ens_members == 1:
+                res.append(worker(j))
+            else:
+                res.append(dask.delayed(worker)(j))
+
+        R_f_ = (
+            dask.compute(*res, num_workers=num_ensemble_workers)
+            if DASK_IMPORTED and n_ens_members > 1
+            else res
+        )
+        res = None
+
+        if is_nowcast_time_step:
+            if measure_time:
+                print(f"{time.time() - starttime:.2f} seconds.")
+            else:
+                print("done.")
+
+        if callback is not None:
+            R_f_stacked = np.stack(R_f_)
+            if R_f_stacked.shape[1] > 0:
+                callback(R_f_stacked.squeeze())
+
+        if return_output:
+            for j in range(n_ens_members):
+                R_f[j].extend(R_f_[j])
+
+        R_f_ = None
+
+    if measure_time:
+        mainloop_time = time.time() - starttime_mainloop
+
+    if return_output:
+        outarr = np.stack([np.stack(R_f[j]) for j in range(n_ens_members)])
+        if measure_time:
+            return outarr, init_time, mainloop_time
+        else:
+            return outarr
+    else:
+        return None
+
+
+def calculate_ratios(correlations):
+    """Calculate explained variance ratios from correlation.
+
+    Parameters
+    ----------
+    Array of shape [component, scale_level, ...]
+      containing correlation (skills) for each component (NWP and nowcast),
+      scale level, and optionally along [y, x] dimensions.
+
+    Returns
+    -------
+    out : numpy array
+      An array containing the ratios of explain variance for each
+      component, scale level, ...
+    """
+    # correlations: [component, scale, ...]
+    square_corrs = np.square(correlations)
+    # Calculate the ratio of the explained variance to the unexplained
+    # variance of the nowcast and NWP model components
+    out = square_corrs / (1 - square_corrs)
+    # out: [component, scale, ...]
+    return out
+
+
+def calculate_weights_bps(correlations):
+    """Calculate BPS blending weights for STEPS blending from correlation.
+
+    Parameters
+    ----------
+    correlations : array-like
+      Array of shape [component, scale_level, ...]
+      containing correlation (skills) for each component (NWP and nowcast),
+      scale level, and optionally along [y, x] dimensions.
+
+    Returns
+    -------
+    weights : array-like
+      Array of shape [component+1, scale_level, ...]
+      containing the weights to be used in STEPS blending for
+      each original component plus an addtional noise component, scale level,
+      and optionally along [y, x] dimensions.
+
+    References
+    ----------
+    :cite:`BPS2006`
+
+    Notes
+    -----
+    The weights in the BPS method can sum op to more than 1.0.
+    """
+    # correlations: [component, scale, ...]
+    # Check if the correlations are positive, otherwise rho = 10e-5
+    correlations = np.where(correlations < 10e-5, 10e-5, correlations)
+
+    # If we merge more than one component with the noise cascade, we follow
+    # the weights impolementation in either :cite:`BPS2006` or :cite:`SPN2013`.
+    if correlations.shape[0] > 1:
+        # Calculate weights for each source
+        ratios = calculate_ratios(correlations)
+        # ratios: [component, scale, ...]
+        total_ratios = np.sum(ratios, axis=0)
+        # total_ratios: [scale, ...] - the denominator of eq. 11 & 12 in BPS2006
+        weights = correlations * np.sqrt(ratios / total_ratios)
+        # weights: [component, scale, ...]
+        # Calculate the weight of the noise component.
+        # Original BPS2006 method in the following two lines (eq. 13)
+        total_square_weights = np.sum(np.square(weights), axis=0)
+        noise_weight = np.sqrt(1.0 - total_square_weights)
+        # Finally, add the noise_weights to the weights variable.
+        weights = np.concatenate((weights, noise_weight[None, ...]), axis=0)
+
+    # Otherwise, the weight equals the correlation on that scale level and
+    # the noise component weight equals 1 - this weight. This only occurs for
+    # the weights calculation outside the radar domain where in the case of 1
+    # NWP model or ensemble member, no blending of multiple models has to take
+    # place
+    else:
+        noise_weight = 1.0 - correlations
+        weights = np.concatenate((correlations, noise_weight), axis=0)
+
+    return weights
+
+
+def calculate_weights_spn(correlations, cov):
+    """Calculate SPN blending weights for STEPS blending from correlation.
+
+    Parameters
+    ----------
+    correlations : array-like
+      Array of shape [n_components]
+      containing correlation (skills) for each component (NWP models and nowcast).
+    cov : array-like
+        Array of shape [n_components, n_components] containing the covariance
+        matrix of the models that will be blended. If cov is set to None and
+        correlations only contains one model, the weight equals the correlation
+        on that scale level and the noise component weight equals 1 - this weight.
+
+    Returns
+    -------
+    weights : array-like
+      Array of shape [component+1]
+      containing the weights to be used in STEPS blending for each original
+      component plus an addtional noise component.
+
+    References
+    ----------
+    :cite:`SPN2013`
+    """
+    # Check if the correlations are positive, otherwise rho = 10e-5
+    correlations = np.where(correlations < 10e-5, 10e-5, correlations)
+
+    if correlations.shape[0] > 1:
+        if isinstance(cov, type(None)):
+            raise ValueError("cov must contain a covariance matrix")
+        else:
+            # Make a numpy matrix out of cov and get the inverse
+            cov_matrix = np.asmatrix(cov)
+            cov_matrix_inv = cov_matrix.getI()
+            # The component weights are the dot product between cov_matrix_inv
+            # and cor_vec
+            weights = cov_matrix_inv.dot(correlations)
+            # If the dot product of the weights with the correlations is
+            # larger than 1.0, we assign a weight of 0.0 to the noise (to make
+            # it numerically stable)
+            if weights.dot(correlations) > 1.0:
+                noise_weight = np.array([0])
+            # Calculate the noise weight
+            else:
+                noise_weight = np.asarray(np.sqrt(1 - weights.dot(correlations)))[0]
+            # Make sure the weights are positive, otherwise weight = 0.0
+            weights = np.where(weights < 0.0, 0.0, weights)[0]
+            # Finally, add the noise_weights to the weights variable.
+            weights = np.concatenate((weights, noise_weight), axis=0)
+
+    # Otherwise, the weight equals the correlation on that scale level and
+    # the noise component weight equals 1 - this weight. This only occurs for
+    # the weights calculation outside the radar domain where in the case of 1
+    # NWP model or ensemble member, no blending of multiple models has to take
+    # place
+    else:
+        noise_weight = 1.0 - correlations
+        weights = np.concatenate((correlations, noise_weight), axis=0)
+
+    return weights
+
+
+def blend_means_sigmas(means, sigmas, weights):
+    """Calculate the blended means and sigmas, the normalization parameters
+    needed to recompose the cascade. This procedure uses the weights of the
+    blending of the normalized cascades and follows eq. 32 and 33 in BPS2004.
+
+    Parameters
+    ----------
+    means : array-like
+      Array of shape [number_components, scale_level, ...]
+      with the mean for each component (NWP, nowcasts, noise).
+    sigmas : array-like
+      Array of shape [number_components, scale_level, ...]
+      with the standard deviation for each component.
+    weights : array-like
+      An array of shape [number_components + 1, scale_level, ...]
+      containing the weights to be used in this routine
+      for each component plus noise, scale level, and optionally [y, x]
+      dimensions, obtained by calling either
+      :py:func:`pysteps.blending.steps.calculate_weights_bps` or
+      :py:func:`pysteps.blending.steps.calculate_weights_spn`.
+
+    Returns
+    -------
+    combined_means : array-like
+      An array of shape [scale_level, ...]
+      containing per scale level (cascade) the weighted combination of
+      means from multiple components (NWP, nowcasts and noise).
+    combined_sigmas : array-like
+      An array of shape [scale_level, ...]
+      similar to combined_means, but containing the standard deviations.
+
+    """
+    # Check if the dimensions are the same
+    diff_dims = weights.ndim - means.ndim
+    if diff_dims:
+        for i in range(diff_dims):
+            means = np.expand_dims(means, axis=means.ndim)
+    diff_dims = weights.ndim - sigmas.ndim
+    if diff_dims:
+        for i in range(diff_dims):
+            sigmas = np.expand_dims(sigmas, axis=sigmas.ndim)
+    # Weight should have one component more (the noise component) than the
+    # means and sigmas. Check this
+    if (
+        weights.shape[0] - means.shape[0] != 1
+        or weights.shape[0] - sigmas.shape[0] != 1
+    ):
+        raise ValueError(
+            "The weights array does not have one (noise) component more than mu and sigma"
+        )
+    else:
+        # Throw away the last component, which is the noise component
+        weights = weights[:-1]
+
+    # Combine (blend) the means and sigmas
+    combined_means = np.zeros(weights.shape[1])
+    combined_sigmas = np.zeros(weights.shape[1])
+    total_weight = np.sum((weights), axis=0)
+    for i in range(weights.shape[0]):
+        combined_means += (weights[i] / total_weight) * means[i]
+        combined_sigmas += (weights[i] / total_weight) * sigmas[i]
+    # TODO: substract covariances to weigthed sigmas - still necessary?
+
+    return combined_means, combined_sigmas
+
+
+def _check_inputs(
+    precip, precip_models, velocity, velocity_models, timesteps, ar_order
+):
+    if precip.ndim != 3:
+        raise ValueError("precip must be a three-dimensional array")
+    if precip.shape[0] < ar_order + 1:
+        raise ValueError("precip.shape[0] < ar_order+1")
+    if precip_models.ndim != 2 and precip_models.ndim != 4:
+        raise ValueError(
+            "precip_models must be either a two-dimensional array containing dictionaries with decomposed model fields or a four-dimensional array containing the original (NWP) model forecasts"
+        )
+    if velocity.ndim != 3:
+        raise ValueError("velocity must be a three-dimensional array")
+    if velocity_models.ndim != 5:
+        raise ValueError("velocity_models must be a five-dimensional array")
+    if velocity.shape[0] != 2 or velocity_models.shape[2] != 2:
+        raise ValueError(
+            "velocity and velocity_models must have an x- and y-component, check the shape"
+        )
+    if precip.shape[1:3] != velocity.shape[1:3]:
+        raise ValueError(
+            "dimension mismatch between precip and velocity: shape(precip)=%s, shape(velocity)=%s"
+            % (str(precip.shape), str(velocity.shape))
+        )
+    if precip_models.shape[0] != velocity_models.shape[0]:
+        raise ValueError(
+            "precip_models and velocity_models must consist of the same number of models"
+        )
+    if isinstance(timesteps, list) and not sorted(timesteps) == timesteps:
+        raise ValueError("timesteps is not in ascending order")
+    if isinstance(timesteps, list):
+        if precip_models.shape[1] != len(timesteps) + 1:
+            raise ValueError(
+                "precip_models does not contain sufficient lead times for this forecast"
+            )
+    else:
+        if precip_models.shape[1] != timesteps + 1:
+            raise ValueError(
+                "precip_models does not contain sufficient lead times for this forecast"
+            )
+
+
+def _compute_incremental_mask(Rbin, kr, r):
+    # buffer the observation mask Rbin using the kernel kr
+    # add a grayscale rim r (for smooth rain/no-rain transition)
+
+    # buffer observation mask
+    Rbin = np.ndarray.astype(Rbin.copy(), "uint8")
+    Rd = scipy.ndimage.morphology.binary_dilation(Rbin, kr)
+
+    # add grayscale rim
+    kr1 = scipy.ndimage.generate_binary_structure(2, 1)
+    mask = Rd.astype(float)
+    for n in range(r):
+        Rd = scipy.ndimage.morphology.binary_dilation(Rd, kr1)
+        mask += Rd
+    # normalize between 0 and 1
+    return mask / mask.max()
+
+
+def _transform_to_lagrangian(
+    precip, velocity, ar_order, xy_coords, extrapolator, extrap_kwargs, num_workers
+):
+    """Advect the previous precipitation fields to the same position with the
+    most recent one (i.e. transform them into the Lagrangian coordinates).
+    """
+    extrap_kwargs = extrap_kwargs.copy()
+    extrap_kwargs["xy_coords"] = xy_coords
+    res = list()
+
+    def f(precip, i):
+        return extrapolator(
+            precip[i, :, :],
+            velocity,
+            ar_order - i,
+            "min",
+            allow_nonfinite_values=True,
+            **extrap_kwargs,
+        )[-1]
+
+    for i in range(ar_order):
+        if not DASK_IMPORTED:
+            precip[i, :, :] = f(precip, i)
+        else:
+            res.append(dask.delayed(f)(precip, i))
+
+    if DASK_IMPORTED:
+        num_workers_ = len(res) if num_workers > len(res) else num_workers
+        precip = np.stack(
+            list(dask.compute(*res, num_workers=num_workers_)) + [precip[-1, :, :]]
+        )
+
+    # replace non-finite values with the minimum value
+    precip = precip.copy()
+    for i in range(precip.shape[0]):
+        precip[i, ~np.isfinite(precip[i, :])] = np.nanmin(precip[i, :])
+    return precip
+
+
+def _init_noise(
+    precip,
+    precip_thr,
+    n_cascade_levels,
+    bp_filter,
+    decompositor,
+    fft,
+    noise_method,
+    noise_kwargs,
+    noise_stddev_adj,
+    measure_time,
+    num_workers,
+):
+    """Initialize the noise method."""
+    if noise_method is None:
+        return None, None, None
+
+    # get methods for perturbations
+    init_noise, generate_noise = noise.get_method(noise_method)
+
+    # initialize the perturbation generator for the precipitation field
+    pp = init_noise(precip, fft_method=fft, **noise_kwargs)
+
+    if noise_stddev_adj == "auto":
+        print("Computing noise adjustment coefficients... ", end="", flush=True)
+        if measure_time:
+            starttime = time.time()
+
+        R_min = np.min(precip)
+        noise_std_coeffs = noise.utils.compute_noise_stddev_adjs(
+            precip[-1, :, :],
+            precip_thr,
+            R_min,
+            bp_filter,
+            decompositor,
+            pp,
+            generate_noise,
+            20,
+            conditional=True,
+            num_workers=num_workers,
+        )
+
+        if measure_time:
+            print(f"{time.time() - starttime:.2f} seconds.")
+        else:
+            print("done.")
+    elif noise_stddev_adj == "fixed":
+        f = lambda k: 1.0 / (0.75 + 0.09 * k)
+        noise_std_coeffs = [f(k) for k in range(1, n_cascade_levels + 1)]
+    else:
+        noise_std_coeffs = np.ones(n_cascade_levels)
+
+    if noise_stddev_adj is not None:
+        print(f"noise std. dev. coeffs:   {noise_std_coeffs}")
+
+    return pp, generate_noise, noise_std_coeffs
+
+
+def _compute_cascade_decomposition_radar(
+    precip,
+    ar_order,
+    n_cascade_levels,
+    n_ens_members,
+    MASK_thr,
+    domain,
+    bp_filter,
+    decompositor,
+    fft,
+):
+    """Compute the cascade decompositions of the input precipitation fields."""
+    R_d = []
+    for i in range(ar_order + 1):
+        R_ = decompositor(
+            precip[i, :, :],
+            bp_filter,
+            mask=MASK_thr,
+            fft_method=fft,
+            output_domain=domain,
+            normalize=True,
+            compute_stats=True,
+            compact_output=True,
+        )
+        R_d.append(R_)
+
+    # Rearrange the cascaded into a four-dimensional array of shape
+    # (n_cascade_levels,ar_order+1,m,n) for the autoregressive model
+    R_c = nowcast_utils.stack_cascades(R_d, n_cascade_levels)
+
+    R_d = R_d[-1]
+    mu_extrapolation = np.array(R_d["means"])
+    sigma_extrapolation = np.array(R_d["stds"])
+    R_d = [R_d.copy() for j in range(n_ens_members)]
+    return R_c, mu_extrapolation, sigma_extrapolation
+
+
+def _compute_cascade_decomposition_nwp(
+    precip_models, bp_filter, decompositor, recompositor, fft, domain
+):
+    """If necessary, decompose (NWP) model forecasts and stack cascades."""
+    if precip_models.ndim == 4:
+        # Keep the model fields for the probability matching later on
+        precip_models_pm = precip_models.copy()
+        # Decompose the (NWP) model forecasts
+        precip_models_cascade = []
+        # Loop through the n_models
+        for i in range(precip_models.shape[0]):
+            precip_model_cascade = []
+            # Loop through the time steps
+            for j in range(precip_models.shape[1]):
+                precip_model_cascade.append(
+                    decompositor(
+                        field=precip_models[i, j, :, :],
+                        bp_filter=bp_filter,
+                        fft_method=fft,
+                        output_domain=domain,
+                        normalize=True,
+                        compute_stats=True,
+                        compact_output=True,
+                    )
+                )
+            precip_models_cascade.append(precip_model_cascade)
+        precip_models_cascade = np.array(precip_models_cascade)
+
+        precip_models, precip_model_cascade = None, None
+    else:
+        precip_models_cascade = precip_models
+        precip_models_pm = None
+
+    # Stack the (NWP) model cascades in separate normalized cascades and return
+    # the means and sigmas.
+    # The normalized model cascade should have the shape:
+    # [n_models, n_timesteps, n_cascade_levels, m, n]
+    precip_models = []
+    mu_models = []
+    sigma_models = []
+    # Stack it per model and combine that
+    for i in range(precip_models_cascade.shape[0]):
+        precip_model, mu_model, sigma_model = blending.utils.stack_cascades(
+            R_d=precip_models_cascade[i, :], donorm=False
+        )
+        precip_models.append(precip_model)
+        mu_models.append(mu_model)
+        sigma_models.append(sigma_model)
+
+    precip_models = np.stack(precip_models)
+    mu_models = np.stack(mu_models)
+    sigma_models = np.stack(sigma_models)
+
+    precip_model, mu_model, sigma_model = None, None, None
+
+    if precip_models_pm is None:
+        # Finally, recompose the (NWP) model cascades to have rainfall fields per
+        # model and time step, which will be used in the probability matching steps.
+        # Recomposed cascade will have shape: [n_models, n_timesteps, m, n]
+        precip_models_pm = []
+        for i in range(precip_models_cascade.shape[0]):
+            precip_model_pm = []
+            for time_step in range(precip_models_cascade.shape[1]):
+                precip_model_pm.append(
+                    recompositor(precip_models_cascade[i, time_step])
+                )
+            precip_models_pm.append(precip_model_pm)
+
+        precip_models_pm = np.stack(precip_models_pm)
+        precip_model_pm = None
+
+    return precip_models, mu_models, sigma_models, precip_models_pm
+
+
+def _estimate_ar_parameters_radar(R_c, ar_order, n_cascade_levels, MASK_thr):
+    """Estimate AR parameters for the radar rainfall field."""
+    # compute lag-l temporal autocorrelation coefficients for each cascade level
+    GAMMA = np.empty((n_cascade_levels, ar_order))
+    for i in range(n_cascade_levels):
+        GAMMA[i, :] = correlation.temporal_autocorrelation(R_c[i], mask=MASK_thr)
+
+    nowcast_utils.print_corrcoefs(GAMMA)
+
+    if ar_order == 2:
+        # adjust the lag-2 correlation coefficient to ensure that the AR(p)
+        # process is stationary
+        for i in range(n_cascade_levels):
+            GAMMA[i, 1] = autoregression.adjust_lag2_corrcoef2(GAMMA[i, 0], GAMMA[i, 1])
+
+    # estimate the parameters of the AR(p) model from the autocorrelation
+    # coefficients
+    PHI = np.empty((n_cascade_levels, ar_order + 1))
+    for i in range(n_cascade_levels):
+        PHI[i, :] = autoregression.estimate_ar_params_yw(GAMMA[i, :])
+
+    nowcast_utils.print_ar_params(PHI)
+    return PHI
+
+
+def _find_nwp_combination(
+    R_c,
+    precip_models,
+    R_models_pm,
+    velocity_models,
+    mu_models,
+    sigma_models,
+    n_ens_members,
+    ar_order,
+    n_cascade_levels,
+    blend_nwp_members,
+):
+    """Determine which (NWP) models will be combined with which nowcast ensemble members.
+    With the way it is implemented at this moment: n_ens_members of the output equals
+    the maximum number of (ensemble) members in the input (either the nowcasts or NWP).
+    """
+    ###
+    # First, discard all except the p-1 last cascades because they are not needed
+    # for the AR(p) model
+    R_c = [R_c[i][-ar_order:] for i in range(n_cascade_levels)]
+
+    # Check if NWP models/members should be used individually, or if all of
+    # them are blended together per nowcast ensemble member.
+    if blend_nwp_members:
+        # stack the extrapolation cascades into a list containing all ensemble members
+        R_c = [
+            [R_c[j].copy() for j in range(n_cascade_levels)]
+            for i in range(n_ens_members)
+        ]
+        R_c = np.stack(R_c)
+        n_model_indices = None
+
+    else:
+        # Start with determining the maximum and mimimum number of members/models
+        # in both input products
+        n_model_members = precip_models.shape[0]
+        n_ens_members_max = max(n_ens_members, n_model_members)
+        n_ens_members_min = min(n_ens_members, n_model_members)
+        # Also make a list of the model index numbers. These indices are needed
+        # for indexing the right climatological skill file when pysteps calculates
+        # the blended forecast in parallel.
+        if n_model_members > 1:
+            n_model_indices = np.arange(n_model_members)
+        else:
+            n_model_indices = [0]
+
+        # Now, repeat the nowcast ensemble members or the nwp models/members until
+        # it has the same amount of members as n_ens_members_max. For instance, if
+        # you have 10 ensemble nowcasts members and 3 NWP members, the output will
+        # be an ensemble of 10 members. Hence, the three NWP members are blended
+        # with the first three members of the nowcast (member one with member one,
+        # two with two, etc.), subsequently, the same NWP members are blended with
+        # the next three members (NWP member one with member 4, NWP member 2 with
+        # member 5, etc.), until 10 is reached.
+        if n_ens_members_min != n_ens_members_max:
+            if n_model_members == 1:
+                precip_models = np.repeat(
+                    precip_models[:, :, :, :, :], n_ens_members_max, axis=0
+                )
+                mu_models = np.repeat(mu_models[:, :, :], n_ens_members_max, axis=0)
+                sigma_models = np.repeat(
+                    sigma_models[:, :, :], n_ens_members_max, axis=0
+                )
+                velocity_models = np.repeat(
+                    velocity_models[:, :, :, :], n_ens_members_max, axis=0
+                )
+                # For the prob. matching
+                R_models_pm = np.repeat(
+                    R_models_pm[:, :, :, :], n_ens_members_max, axis=0
+                )
+                # Finally, for the model indices
+                n_model_indices = np.repeat(n_model_indices, n_ens_members_max, axis=0)
+
+            elif n_model_members == n_ens_members_min:
+                repeats = [
+                    (n_ens_members_max + i) // n_ens_members_min
+                    for i in range(n_ens_members_min)
+                ]
+                if n_model_members == n_ens_members_min:
+                    precip_models = np.repeat(precip_models, repeats, axis=0)
+                    mu_models = np.repeat(mu_models, repeats, axis=0)
+                    sigma_models = np.repeat(sigma_models, repeats, axis=0)
+                    velocity_models = np.repeat(velocity_models, repeats, axis=0)
+                    # For the prob. matching
+                    R_models_pm = np.repeat(R_models_pm, repeats, axis=0)
+                    # Finally, for the model indices
+                    n_model_indices = np.repeat(n_model_indices, repeats, axis=0)
+
+        R_c = [
+            [R_c[j].copy() for j in range(n_cascade_levels)]
+            for i in range(n_ens_members_max)
+        ]
+        R_c = np.stack(R_c)
+
+        n_ens_members = n_ens_members_max
+
+    return (
+        R_c,
+        precip_models,
+        R_models_pm,
+        velocity_models,
+        mu_models,
+        sigma_models,
+        n_ens_members,
+        n_model_indices,
+    )
+
+
+def _init_random_generators(
+    velocity,
+    noise_method,
+    vel_pert_method,
+    vp_par,
+    vp_perp,
+    seed,
+    n_ens_members,
+    kmperpixel,
+    timestep,
+):
+    """Initialize all the random generators."""
+    if noise_method is not None:
+        randgen_prec = []
+        randgen_motion = []
+        np.random.seed(seed)
+        for j in range(n_ens_members):
+            rs = np.random.RandomState(seed)
+            randgen_prec.append(rs)
+            seed = rs.randint(0, high=1e9)
+            rs = np.random.RandomState(seed)
+            randgen_motion.append(rs)
+            seed = rs.randint(0, high=1e9)
+
+    if vel_pert_method is not None:
+        init_vel_noise, generate_vel_noise = noise.get_method(vel_pert_method)
+
+        # initialize the perturbation generators for the motion field
+        vps = []
+        for j in range(n_ens_members):
+            kwargs = {
+                "randstate": randgen_motion[j],
+                "p_par": vp_par,
+                "p_perp": vp_perp,
+            }
+            vp_ = init_vel_noise(velocity, 1.0 / kmperpixel, timestep, **kwargs)
+            vps.append(vp_)
+    else:
+        vps, generate_vel_noise = None, None
+
+    return randgen_prec, vps, generate_vel_noise
+
+
+def _prepare_forecast_loop(
+    R_c,
+    noise_method,
+    fft_method,
+    n_cascade_levels,
+    n_ens_members,
+    mask_method,
+    mask_kwargs,
+    timestep,
+    kmperpixel,
+):
+    """Prepare for the forecast loop."""
+    # Empty arrays for the previous displacements and the forecast cascade
+    D = np.stack([np.full(n_cascade_levels, None) for j in range(n_ens_members)])
+    D_Yn = np.stack([np.full(n_cascade_levels, None) for j in range(n_ens_members)])
+    D_pb = np.stack([None for j in range(n_ens_members)])
+    R_f = [[] for j in range(n_ens_members)]
+
+    if mask_method == "incremental":
+        # get mask parameters
+        mask_rim = mask_kwargs.get("mask_rim", 10)
+        mask_f = mask_kwargs.get("mask_f", 1.0)
+        # initialize the structuring element
+        struct = scipy.ndimage.generate_binary_structure(2, 1)
+        # iterate it to expand it nxn
+        n = mask_f * timestep / kmperpixel
+        struct = scipy.ndimage.iterate_structure(struct, int((n - 1) / 2.0))
+    else:
+        mask_rim, struct = None, None
+
+    if noise_method is None:
+        R_m = [R_c[0][i].copy() for i in range(n_cascade_levels)]
+    else:
+        R_m = None
+
+    fft_objs = []
+    for i in range(n_ens_members):
+        fft_objs.append(utils.get_method(fft_method, shape=R_c.shape[-2:]))
+
+    return D, D_Yn, D_pb, R_f, R_m, mask_rim, struct, fft_objs
+
+
+def _compute_initial_nwp_skill(
+    R_c, precip_models, domain_mask, issuetime, outdir_path_skill, clim_kwargs
+):
+    """Calculate the initial skill of the (NWP) model forecasts at t=0."""
+    rho_nwp_models = [
+        blending.skill_scores.spatial_correlation(
+            obs=R_c[0, :, -1, :, :],
+            mod=precip_models[n_model, 0, :, :, :],
+            domain_mask=domain_mask,
+        )
+        for n_model in range(precip_models.shape[0])
+    ]
+    rho_nwp_models = np.stack(rho_nwp_models)
+
+    # Ensure that the model skill decreases with increasing scale level.
+    for n_model in range(precip_models.shape[0]):
+        for i in range(1, precip_models.shape[2]):
+            if rho_nwp_models[n_model, i] > rho_nwp_models[n_model, i - 1]:
+                # Set it equal to the previous scale level
+                rho_nwp_models[n_model, i] = rho_nwp_models[n_model, i - 1]
+
+    # Save this in the climatological skill file
+    blending.clim.save_skill(
+        current_skill=rho_nwp_models,
+        validtime=issuetime,
+        outdir_path=outdir_path_skill,
+        **clim_kwargs,
+    )
+    return rho_nwp_models
diff --git a/pysteps/blending/utils.py b/pysteps/blending/utils.py
new file mode 100644
index 000000000..6a63890ea
--- /dev/null
+++ b/pysteps/blending/utils.py
@@ -0,0 +1,510 @@
+# -*- coding: utf-8 -*-
+"""
+pysteps.blending.utils
+======================
+
+Module with common utilities used by the blending methods.
+
+.. autosummary::
+    :toctree: ../generated/
+
+    stack_cascades
+    blend_cascades
+    recompose_cascade
+    blend_optical_flows
+    decompose_NWP
+    compute_store_nwp_motion
+    load_NWP
+"""
+
+import datetime
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+import netCDF4
+
+from pysteps.cascade import get_method as cascade_get_method
+from pysteps.cascade.bandpass_filters import filter_gaussian
+from pysteps.utils import get_method as utils_get_method
+
+
+def stack_cascades(R_d, donorm=True):
+    """Stack the given cascades into a larger array.
+
+    Parameters
+    ----------
+    R_d : dict
+      Dictionary containing a list of cascades obtained by calling a method
+      implemented in pysteps.cascade.decomposition.
+    donorm : bool
+      If True, normalize the cascade levels before stacking.
+
+    Returns
+    -------
+    out : tuple
+      A three-element tuple containing a four-dimensional array of stacked
+      cascade levels and arrays of mean values and standard deviations for each
+      cascade level.
+    """
+    R_c = []
+    mu_c = []
+    sigma_c = []
+
+    for cascade in R_d:
+        R_ = []
+        R_i = cascade["cascade_levels"]
+        n_levels = R_i.shape[0]
+        mu_ = np.asarray(cascade["means"])
+        sigma_ = np.asarray(cascade["stds"])
+        if donorm:
+            for j in range(n_levels):
+                R__ = (R_i[j, :, :] - mu_[j]) / sigma_[j]
+                R_.append(R__)
+        else:
+            R_ = R_i
+        R_c.append(np.stack(R_))
+        mu_c.append(mu_)
+        sigma_c.append(sigma_)
+    return np.stack(R_c), np.stack(mu_c), np.stack(sigma_c)
+
+
+def blend_cascades(cascades_norm, weights):
+    """Calculate blended normalized cascades using STEPS weights following eq.
+    10 in :cite:`BPS2006`.
+
+    Parameters
+    ----------
+    cascades_norm : array-like
+      Array of shape [number_components + 1, scale_level, ...]
+      with the cascade for each component (NWP, nowcasts, noise) and scale level,
+      obtained by calling a method implemented in pysteps.blending.utils.stack_cascades
+
+    weights : array-like
+      An array of shape [number_components + 1, scale_level, ...]
+      containing the weights to be used in this routine
+      for each component plus noise, scale level, and optionally [y, x]
+      dimensions, obtained by calling a method implemented in
+      pysteps.blending.steps.calculate_weights
+
+    Returns
+    -------
+    combined_cascade : array-like
+      An array of shape [scale_level, y, x]
+      containing per scale level (cascade) the weighted combination of
+      cascades from multiple components (NWP, nowcasts and noise) to be used
+      in STEPS blending.
+    """
+    # check inputs
+    if isinstance(cascades_norm, (list, tuple)):
+        cascades_norm = np.stack(cascades_norm)
+
+    if isinstance(weights, (list, tuple)):
+        weights = np.asarray(weights)
+
+    # check weights dimensions match number of sources
+    num_sources = cascades_norm.shape[0]
+    num_sources_klevels = cascades_norm.shape[1]
+    num_weights = weights.shape[0]
+    num_weights_klevels = weights.shape[1]
+
+    if num_weights != num_sources:
+        raise ValueError(
+            "dimension mismatch between cascades and weights.\n"
+            "weights dimension must match the number of components in cascades.\n"
+            f"number of models={num_sources}, number of weights={num_weights}"
+        )
+    if num_weights_klevels != num_sources_klevels:
+        raise ValueError(
+            "dimension mismatch between cascades and weights.\n"
+            "weights cascade levels dimension must match the number of cascades in cascades_norm.\n"
+            f"number of cascade levels={num_sources_klevels}, number of weights={num_weights_klevels}"
+        )
+
+    # cascade_norm component, scales, y, x
+    # weights component, scales, ....
+    # Reshape weights to make the calculation possible with numpy
+    all_c_wn = weights.reshape(num_weights, num_weights_klevels, 1, 1) * cascades_norm
+    combined_cascade = np.sum(all_c_wn, axis=0)
+    # combined_cascade [scale, ...]
+    return combined_cascade
+
+
+def recompose_cascade(combined_cascade, combined_mean, combined_sigma):
+    """Recompose the cascades into a transformed rain rate field.
+
+
+    Parameters
+    ----------
+    combined_cascade : array-like
+      An array of shape [scale_level, y, x]
+      containing per scale level (cascade) the weighted combination of
+      cascades from multiple components (NWP, nowcasts and noise) to be used
+      in STEPS blending.
+    combined_mean : array-like
+      An array of shape [scale_level, ...]
+      similar to combined_cascade, but containing the normalization parameter
+      mean.
+    combined_sigma : array-like
+      An array of shape [scale_level, ...]
+      similar to combined_cascade, but containing the normalization parameter
+      standard deviation.
+
+    Returns
+    -------
+    out: array-like
+        A two-dimensional array containing the recomposed cascade.
+
+    """
+    # Renormalize with the blended sigma and mean values
+    renorm = (
+        combined_cascade * combined_sigma.reshape(combined_cascade.shape[0], 1, 1)
+    ) + combined_mean.reshape(combined_mean.shape[0], 1, 1)
+    # print(renorm.shape)
+    out = np.sum(renorm, axis=0)
+    # print(out.shape)
+    return out
+
+
+def blend_optical_flows(flows, weights):
+    """Combine advection fields using given weights. Following :cite:`BPS2006`
+    the second level of the cascade is used for the weights
+
+    Parameters
+    ----------
+    flows : array-like
+      A stack of multiple advenction fields having shape
+      (S, 2, m, n), where flows[N, :, :, :] contains the motion vectors
+      for source N.
+      Advection fields for each source can be obtanined by
+      calling any of the methods implemented in
+      pysteps.motion and then stack all together
+    weights : array-like
+      An array of shape [number_sources]
+      containing the weights to be used to combine
+      the advection fields of each source.
+      weights are modified to make their sum equal to one.
+    Returns
+    -------
+    out: ndarray
+        Return the blended advection field having shape
+        (2, m, n), where out[0, :, :] contains the x-components of
+        the blended motion vectors and out[1, :, :] contains the y-components.
+        The velocities are in units of pixels / timestep.
+    """
+
+    # check inputs
+    if isinstance(flows, (list, tuple)):
+        flows = np.stack(flows)
+
+    if isinstance(weights, (list, tuple)):
+        weights = np.asarray(weights)
+
+    # check weights dimensions match number of sources
+    num_sources = flows.shape[0]
+    num_weights = weights.shape[0]
+
+    if num_weights != num_sources:
+        raise ValueError(
+            "dimension mismatch between flows and weights.\n"
+            "weights dimension must match the number of flows.\n"
+            f"number of flows={num_sources}, number of weights={num_weights}"
+        )
+    # normalize weigths
+    weights = weights / np.sum(weights)
+
+    # flows dimension sources, 2, m, n
+    # weights dimension sources
+    # move source axis to last to allow broadcasting
+    # TODO: Check if broadcasting has worked well
+    all_c_wn = weights * np.moveaxis(flows, 0, -1)
+    # sum uses last axis
+    combined_flows = np.sum(all_c_wn, axis=-1)
+    # combined_flows [2, m, n]
+    return combined_flows
+
+
+def decompose_NWP(
+    R_NWP,
+    NWP_model,
+    analysis_time,
+    timestep,
+    valid_times,
+    output_path,
+    num_cascade_levels=8,
+    num_workers=1,
+    decomp_method="fft",
+    fft_method="numpy",
+    domain="spatial",
+    normalize=True,
+    compute_stats=True,
+    compact_output=True,
+):
+    """Decomposes the NWP forecast data into cascades and saves it in
+    a netCDF file
+
+    Parameters
+    ----------
+    R_NWP: array-like
+      Array of dimension (n_timesteps, x, y) containing the precipitation forecast
+      from some NWP model.
+    NWP_model: str
+      The name of the NWP model
+    analysis_time: numpy.datetime64
+      The analysis time of the NWP forecast. The analysis time is assumed to be a
+      numpy.datetime64 type as imported by the pysteps importer
+    timestep: int
+      Timestep in minutes between subsequent NWP forecast fields
+    valid_times: array_like
+      Array containing the valid times of the NWP forecast fields. The times are
+      assumed to be numpy.datetime64 types as imported by the pysteps importer.
+    output_path: str
+      The location where to save the file with the NWP cascade. Defaults to the
+      path_workdir specified in the rcparams file.
+    num_cascade_levels: int, optional
+      The number of frequency bands to use. Must be greater than 2. Defaults to 8.
+    num_workers: int, optional
+      The number of workers to use for parallel computation. Applicable if dask
+      is enabled or pyFFTW is used for computing the FFT. When num_workers>1, it
+      is advisable to disable OpenMP by setting the environment variable
+      OMP_NUM_THREADS to 1. This avoids slowdown caused by too many simultaneous
+      threads.
+
+    Other Parameters
+    ----------------
+    decomp_method: str, optional
+      A string defining the decomposition method to use. Defaults to "fft".
+    fft_method: str or tuple, optional
+      A string or a (function,kwargs) tuple defining the FFT method to use
+      (see :py:func:`pysteps.utils.interface.get_method`).
+      Defaults to "numpy". This option is not used if input_domain and
+      output_domain are both set to "spectral".
+    domain: {"spatial", "spectral"}, optional
+      If "spatial", the output cascade levels are transformed back to the
+      spatial domain by using the inverse FFT. If "spectral", the cascade is
+      kept in the spectral domain. Defaults to "spatial".
+    normalize: bool, optional
+      If True, normalize the cascade levels to zero mean and unit variance.
+      Requires that compute_stats is True. Implies that compute_stats is True.
+      Defaults to False.
+    compute_stats: bool, optional
+      If True, the output dictionary contains the keys "means" and "stds"
+      for the mean and standard deviation of each output cascade level.
+      Defaults to False.
+    compact_output: bool, optional
+      Applicable if output_domain is "spectral". If set to True, only the
+      parts of the Fourier spectrum with non-negligible filter weights are
+      stored. Defaults to False.
+
+
+    Returns
+    -------
+    Nothing
+    """
+
+    # Make a NetCDF file
+    output_date = datetime.datetime.strftime(
+        pd.to_datetime(analysis_time), "%Y%m%d%H%M%S"
+    )
+    outfn = Path(output_path) / f"cascade_{NWP_model}_{output_date}.nc"
+    ncf = netCDF4.Dataset(outfn, "w", format="NETCDF4")
+
+    # Express times relative to the zero time
+    zero_time = np.datetime64("1970-01-01T00:00:00", "ns")
+    valid_times = np.array(valid_times) - zero_time
+    analysis_time = analysis_time - zero_time
+
+    # Set attributes of decomposition method
+    ncf.domain = domain
+    ncf.normalized = int(normalize)
+    ncf.compact_output = int(compact_output)
+    ncf.analysis_time = int(analysis_time)
+    ncf.timestep = int(timestep)
+
+    # Create dimensions
+    ncf.createDimension("time", R_NWP.shape[0])
+    ncf.createDimension("cascade_levels", num_cascade_levels)
+    ncf.createDimension("x", R_NWP.shape[2])
+    ncf.createDimension("y", R_NWP.shape[1])
+
+    # Create variables (decomposed cascade, means and standard deviations)
+    R_d = ncf.createVariable(
+        "pr_decomposed",
+        np.float32,
+        ("time", "cascade_levels", "y", "x"),
+        zlib=True,
+        complevel=4,
+    )
+    means = ncf.createVariable("means", np.float64, ("time", "cascade_levels"))
+    stds = ncf.createVariable("stds", np.float64, ("time", "cascade_levels"))
+    v_times = ncf.createVariable("valid_times", np.float64, ("time",))
+    v_times.units = "nanoseconds since 1970-01-01 00:00:00"
+
+    # The valid times are saved as an array of floats, because netCDF files can't handle datetime types
+    v_times[:] = np.array([np.float64(valid_times[i]) for i in range(len(valid_times))])
+
+    # Decompose the NWP data
+    filter_g = filter_gaussian(R_NWP.shape[1:], num_cascade_levels)
+    fft = utils_get_method(fft_method, shape=R_NWP.shape[1:], n_threads=num_workers)
+    decomp_method, _ = cascade_get_method(decomp_method)
+
+    for i in range(R_NWP.shape[0]):
+        R_ = decomp_method(
+            field=R_NWP[i, :, :],
+            bp_filter=filter_g,
+            fft_method=fft,
+            input_domain=domain,
+            output_domain=domain,
+            normalize=normalize,
+            compute_stats=compute_stats,
+            compact_output=compact_output,
+        )
+
+        # Save data to netCDF file
+        # print(R_["cascade_levels"])
+        R_d[i, :, :, :] = R_["cascade_levels"]
+        means[i, :] = R_["means"]
+        stds[i, :] = R_["stds"]
+
+    # Close the file
+    ncf.close()
+
+
+def compute_store_nwp_motion(
+    precip_nwp,
+    oflow_method,
+    analysis_time,
+    nwp_model,
+    output_path,
+):
+    """Computes, per forecast lead time, the velocity field of an NWP model field.
+
+    Parameters
+    ----------
+    precip_nwp: array-like
+      Array of dimension (n_timesteps, x, y) containing the precipitation forecast
+      from some NWP model.
+    oflow_method: {'constant', 'darts', 'lucaskanade', 'proesmans', 'vet'}, optional
+      An optical flow method from pysteps.motion.get_method.
+    analysis_time: numpy.datetime64
+      The analysis time of the NWP forecast. The analysis time is assumed to be a
+      numpy.datetime64 type as imported by the pysteps importer.
+    nwp_model: str
+      The name of the NWP model.
+    output_path: str, optional
+      The location where to save the file with the NWP velocity fields. Defaults
+      to the path_workdir specified in the rcparams file.
+
+    Returns
+    -------
+    Nothing
+    """
+
+    # Set the output file
+    output_date = datetime.datetime.strftime(
+        pd.to_datetime(analysis_time), "%Y%m%d%H%M%S"
+    )
+    outfn = Path(output_path) / f"motion_{nwp_model}_{output_date}.npy"
+
+    # Get the velocity field per time step
+    v_nwp = np.zeros((precip_nwp.shape[0], 2, precip_nwp.shape[1], precip_nwp.shape[2]))
+    # Loop through the timesteps. We need two images to construct a motion
+    # field, so we can start from timestep 1.
+    for t in range(1, precip_nwp.shape[0]):
+        v_nwp[t] = oflow_method(precip_nwp[t - 1 : t + 1, :, :])
+
+    # Make timestep 0 the same as timestep 1.
+    v_nwp[0] = v_nwp[1]
+
+    assert v_nwp.ndim == 4, "v_nwp must be a four-dimensional array"
+
+    # Save it as a numpy array
+    np.save(outfn, v_nwp)
+
+
+def load_NWP(input_nc_path_decomp, input_path_velocities, start_time, n_timesteps):
+    """Loads the decomposed NWP and velocity data from the netCDF files
+
+    Parameters
+    ----------
+    input_nc_path_decomp: str
+      Path to the saved netCDF file containing the decomposed NWP data.
+    input_path_velocities: str
+        Path to the saved numpy binary file containing the estimated velocity
+        fields from the NWP data.
+    start_time: numpy.datetime64
+      The start time of the nowcasting. Assumed to be a numpy.datetime64 type
+    n_timesteps: int
+      Number of time steps to forecast
+
+    Returns
+    -------
+    R_d: list
+      A list of dictionaries with each element in the list corresponding to
+      a different time step. Each dictionary has the same structure as the
+      output of the decomposition function
+    uv: array-like
+        Array of shape (timestep,2,m,n) containing the x- and y-components
+      of the advection field for the (NWP) model field per forecast lead time.
+    """
+
+    # Open the file
+    ncf_decomp = netCDF4.Dataset(input_nc_path_decomp, "r", format="NETCDF4")
+    velocities = np.load(input_path_velocities)
+
+    # Initialise the decomposition dictionary
+    decomp_dict = dict()
+    decomp_dict["domain"] = ncf_decomp.domain
+    decomp_dict["normalized"] = bool(ncf_decomp.normalized)
+    decomp_dict["compact_output"] = bool(ncf_decomp.compact_output)
+
+    # Convert the start time and the timestep to datetime64 and timedelta64 type
+    zero_time = np.datetime64("1970-01-01T00:00:00", "ns")
+    analysis_time = np.timedelta64(int(ncf_decomp.analysis_time), "ns") + zero_time
+
+    timestep = ncf_decomp.timestep
+    timestep = np.timedelta64(timestep, "m")
+
+    valid_times = ncf_decomp.variables["valid_times"][:]
+    valid_times = np.array(
+        [np.timedelta64(int(valid_times[i]), "ns") for i in range(len(valid_times))]
+    )
+    valid_times = valid_times + zero_time
+
+    # Find the indices corresponding with the required start and end time
+    start_i = (start_time - analysis_time) // timestep
+    assert analysis_time + start_i * timestep == start_time
+    end_i = start_i + n_timesteps + 1
+
+    # Add the valid times to the output
+    decomp_dict["valid_times"] = valid_times[start_i:end_i]
+
+    # Slice the velocity fields with the start and end indices
+    uv = velocities[start_i:end_i, :, :, :]
+
+    # Initialise the list of dictionaries which will serve as the output (cf: the STEPS function)
+    R_d = list()
+
+    for i in range(start_i, end_i):
+        decomp_dict_ = decomp_dict.copy()
+
+        cascade_levels = ncf_decomp.variables["pr_decomposed"][i, :, :, :]
+
+        # In the netcdf file this is saved as a masked array, so we're checking if there is no mask
+        assert not cascade_levels.mask
+
+        means = ncf_decomp.variables["means"][i, :]
+        assert not means.mask
+
+        stds = ncf_decomp.variables["stds"][i, :]
+        assert not stds.mask
+
+        # Save the values in the dictionary as normal arrays with the filled method
+        decomp_dict_["cascade_levels"] = np.ma.filled(cascade_levels)
+        decomp_dict_["means"] = np.ma.filled(means)
+        decomp_dict_["stds"] = np.ma.filled(stds)
+
+        # Append the output list
+        R_d.append(decomp_dict_)
+
+    return R_d, uv
diff --git a/pysteps/decorators.py b/pysteps/decorators.py
index 5fdf799c8..9bfc3fad5 100644
--- a/pysteps/decorators.py
+++ b/pysteps/decorators.py
@@ -43,10 +43,8 @@ def postprocess_import(fillna=np.nan, dtype="double"):
     """
     Postprocess the imported precipitation data.
     Operations:
-
     - Allow type casting (dtype keyword)
     - Set invalid or missing data to predefined value (fillna keyword)
-
     This decorator replaces the text "{extra_kwargs}" in the function's
     docstring with the documentation of the keywords used in the postprocessing.
     The additional docstrings are added as "Other Parameters" in the importer function.
diff --git a/pysteps/extrapolation/semilagrangian.py b/pysteps/extrapolation/semilagrangian.py
index d565428ed..425fcfbf9 100644
--- a/pysteps/extrapolation/semilagrangian.py
+++ b/pysteps/extrapolation/semilagrangian.py
@@ -132,6 +132,7 @@ def extrapolate(
     n_iter = kwargs.get("n_iter", 1)
     return_displacement = kwargs.get("return_displacement", False)
     interp_order = kwargs.get("interp_order", 1)
+    map_coordinates_mode = kwargs.get("map_coordinates_mode", "constant")
 
     if precip is None and not return_displacement:
         raise ValueError("precip is None but return_displacement is False")
@@ -223,7 +224,7 @@ def interpolate_motion(displacement, velocity_inc, td):
             precip_warped = ip.map_coordinates(
                 precip,
                 coords_warped,
-                mode="constant",
+                mode=map_coordinates_mode,
                 cval=outval,
                 order=interp_order,
                 prefilter=prefilter,
@@ -233,7 +234,7 @@ def interpolate_motion(displacement, velocity_inc, td):
                 mask_warped = ip.map_coordinates(
                     mask_min,
                     coords_warped,
-                    mode="constant",
+                    mode=map_coordinates_mode,
                     cval=0,
                     order=1,
                     prefilter=False,
@@ -244,7 +245,7 @@ def interpolate_motion(displacement, velocity_inc, td):
                     mask_warped = ip.map_coordinates(
                         mask_finite,
                         coords_warped,
-                        mode="constant",
+                        mode=map_coordinates_mode,
                         cval=0,
                         order=1,
                         prefilter=False,
diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py
index 86a18b289..47ae3fe6d 100644
--- a/pysteps/io/importers.py
+++ b/pysteps/io/importers.py
@@ -90,12 +90,13 @@
 from functools import partial
 
 import numpy as np
+
+# import xarray as xr
 from matplotlib.pyplot import imread
 
 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 +235,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 +252,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 +277,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 +304,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 +334,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 +359,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)
@@ -445,11 +411,8 @@ def import_bom_rf3(filename, **kwargs):
             "but it is not installed"
         )
 
-    precip = _import_bom_rf3_data(filename)
-
-    geodata = _import_bom_rf3_geodata(filename)
+    precip, geodata = _import_bom_rf3_data(filename)
     metadata = geodata
-    # TODO(import_bom_rf3): Add missing georeferencing data.
 
     metadata["transform"] = None
     metadata["zerovalue"] = np.nanmin(precip)
@@ -460,20 +423,19 @@ def import_bom_rf3(filename, **kwargs):
 
 def _import_bom_rf3_data(filename):
     ds_rainfall = netCDF4.Dataset(filename)
+    geodata = _import_bom_rf3_geodata(ds_rainfall)
     if "precipitation" in ds_rainfall.variables.keys():
         precipitation = ds_rainfall.variables["precipitation"][:]
     else:
         precipitation = None
     ds_rainfall.close()
 
-    return precipitation
+    return precipitation, geodata
 
 
-def _import_bom_rf3_geodata(filename):
+def _import_bom_rf3_geodata(ds_rainfall):
     geodata = {}
 
-    ds_rainfall = netCDF4.Dataset(filename)
-
     if "proj" in ds_rainfall.variables.keys():
         projection = ds_rainfall.variables["proj"]
         if getattr(projection, "grid_mapping_name") == "albers_conical_equal_area":
@@ -514,7 +476,7 @@ def _import_bom_rf3_geodata(filename):
     geodata["xpixelsize"] = xpixelsize * factor_scale
     geodata["ypixelsize"] = ypixelsize * factor_scale
     geodata["cartesian_unit"] = "m"
-    geodata["yorigin"] = "upper"  # TODO(_import_bom_rf3_geodata): check this
+    geodata["yorigin"] = "upper"
 
     # get the accumulation period
     valid_time = None
@@ -549,7 +511,6 @@ def _import_bom_rf3_geodata(filename):
             geodata["unit"] = "mm"
 
     geodata["institution"] = "Commonwealth of Australia, Bureau of Meteorology"
-    ds_rainfall.close()
 
     return geodata
 
@@ -749,7 +710,7 @@ def _import_fmi_pgm_metadata(filename, gzipped=False):
 
 
 @postprocess_import()
-def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1.0, **kwargs):
+def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1000.0, **kwargs):
     """
     Import a precipitation or reflectivity field (and optionally the quality
     field) from a HDF5 file conforming to the KNMI Data Centre specification.
@@ -767,8 +728,8 @@ def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1.0, **kwargs
         is used as default, but hourly, daily and monthly accumulations
         are also available.
     pixelsize: float
-        The pixel size of a raster cell in kilometers. The default value for the
-        KNMI datasets is a 1 km grid cell size, but datasets with 2.4 km pixel
+        The pixel size of a raster cell in meters. The default value for the
+        KNMI datasets is a 1000 m grid cell size, but datasets with 2400 m pixel
         size are also available.
 
     {extra_kwargs_doc}
@@ -846,7 +807,7 @@ def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1.0, **kwargs
     # The 'where' group of mch- and Opera-data, is called 'geographic' in the
     # KNMI data.
     geographic = f["geographic"]
-    proj4str = geographic["map_projection"].attrs["projection_proj4_params"].decode()
+    proj4str = "+proj=stere +lat_0=90 +lon_0=0.0 +lat_ts=60.0 +a=6378137 +b=6356752 +x_0=0 +y_0=0"
     pr = pyproj.Proj(proj4str)
     metadata["projection"] = proj4str
 
@@ -877,7 +838,7 @@ def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1.0, **kwargs
     metadata["y2"] = y2
     metadata["xpixelsize"] = pixelsize
     metadata["ypixelsize"] = pixelsize
-    metadata["cartesian_unit"] = "km"
+    metadata["cartesian_unit"] = "m"
     metadata["yorigin"] = "upper"
     metadata["institution"] = "KNMI - Royal Netherlands Meteorological Institute"
     metadata["accutime"] = accutime
@@ -1419,9 +1380,18 @@ def import_odim_hdf5(filename, qty="RATE", **kwargs):
         x2 = ur_x
         y2 = ur_y
 
+    dataset1 = f["dataset1"]
+
     if "xscale" in where.attrs.keys() and "yscale" in where.attrs.keys():
         xpixelsize = where.attrs["xscale"]
         ypixelsize = where.attrs["yscale"]
+    elif (
+        "xscale" in dataset1["where"].attrs.keys()
+        and "yscale" in dataset1["where"].attrs.keys()
+    ):
+        where = dataset1["where"]
+        xpixelsize = where.attrs["xscale"]
+        ypixelsize = where.attrs["yscale"]
     else:
         xpixelsize = None
         ypixelsize = None
@@ -1458,6 +1428,8 @@ def import_odim_hdf5(filename, qty="RATE", **kwargs):
         "threshold": _get_threshold_value(precip),
     }
 
+    metadata.update(kwargs)
+
     f.close()
 
     return precip, quality, metadata
diff --git a/pysteps/nowcasts/interface.py b/pysteps/nowcasts/interface.py
index 1d961076a..2194d0905 100644
--- a/pysteps/nowcasts/interface.py
+++ b/pysteps/nowcasts/interface.py
@@ -32,7 +32,14 @@
 """
 
 from pysteps.extrapolation.interface import eulerian_persistence
-from pysteps.nowcasts import anvil, extrapolation, linda, sprog, steps, sseps
+from pysteps.nowcasts import (
+    anvil,
+    extrapolation,
+    linda,
+    sprog,
+    steps,
+    sseps,
+)
 from pysteps.nowcasts import lagrangian_probability
 
 _nowcast_methods = dict()
diff --git a/pysteps/pystepsrc b/pysteps/pystepsrc
index 7e2b2c9d7..0fcaff302 100644
--- a/pysteps/pystepsrc
+++ b/pysteps/pystepsrc
@@ -76,8 +76,19 @@
             "importer_kwargs": {
                 "accutime": 5,
                 "qty": "ACRR",
-                "pixelsize": 1.0
+                "pixelsize": 1000.0
 			}
+        },
+		"rmi": {
+            "root_path": "./radar/rmi/radqpe",
+            "path_fmt": "%Y%m%d",
+            "fn_pattern": "%Y%m%d%H%M00.rad.bhbjbwdnfa.comp.rate.qpe2",
+            "fn_ext": "hdf",
+            "importer": "odim_hdf5",
+            "timestep": 5,
+            "importer_kwargs": {
+				"accutime": 5.0
+            }
         },
         "saf": {
             "root_path": "./saf",
@@ -89,6 +100,39 @@
             "importer_kwargs": {
                 "gzipped": true
             }
+        },
+        "bom_nwp": {
+            "root_path": "./nwp/bom",
+            "path_fmt": "%Y/%m/%d",
+            "fn_pattern": "%Y%m%d_%H00_regrid_short",
+            "fn_ext": "nc",
+            "importer": "bom_nwp",
+            "timestep": 10,
+            "importer_kwargs": {
+                "gzipped": true
+            }
+        },
+        "rmi_nwp": {
+            "root_path": "./nwp/rmi",
+            "path_fmt": "%Y/%m/%d",
+            "fn_pattern": "ao13_%Y%m%d%H_native_5min",
+            "fn_ext": "nc",
+            "importer": "rmi_nwp",
+            "timestep": 5,
+            "importer_kwargs": {
+                "gzipped": true
+            }
+        },
+        "knmi_nwp": {
+            "root_path": "./nwp/knmi",
+            "path_fmt": "%Y/%m/%d",
+            "fn_pattern": "%Y%m%d_%H00_Pforecast_Harmonie",
+            "fn_ext": "nc",
+            "importer": "knmi_nwp",
+            "timestep": 60,
+            "importer_kwargs": {
+                "gzipped": true
+            }
         }
     }
 }
diff --git a/pysteps/utils/__init__.py b/pysteps/utils/__init__.py
index 214d9a309..9594a75ae 100644
--- a/pysteps/utils/__init__.py
+++ b/pysteps/utils/__init__.py
@@ -11,3 +11,4 @@
 from .spectral import *
 from .tapering import *
 from .transformation import *
+from .reprojection import *
diff --git a/pysteps/utils/interface.py b/pysteps/utils/interface.py
index d68ff27f9..de1b81d0b 100644
--- a/pysteps/utils/interface.py
+++ b/pysteps/utils/interface.py
@@ -18,6 +18,7 @@
 from . import fft
 from . import images
 from . import interpolate
+from . import reprojection
 from . import spectral
 from . import tapering
 from . import transformation
@@ -108,6 +109,14 @@ def get_method(name, **kwargs):
     Additional keyword arguments are passed to the initializer of the FFT
     methods, see utils.fft.
 
+    Reprojection methods:
+
+    +-------------------+-----------------------------------------------------+
+    |     Name          |              Description                            |
+    +===================+=====================================================+
+    |  reproject_grids  | Reproject grids to a destination grid.              |
+    +-------------------+-----------------------------------------------------+
+
     Spectral methods:
 
     +-------------------+-----------------------------------------------------+
@@ -188,6 +197,9 @@ def donothing(R, metadata=None, *args, **kwargs):
     methods_objects["rbfinterp2d"] = interpolate.rbfinterp2d
     methods_objects["idwinterp2d"] = interpolate.idwinterp2d
 
+    # reprojection methods
+    methods_objects["reproject_grids"] = reprojection.reproject_grids
+
     # spectral methods
     methods_objects["rapsd"] = spectral.rapsd
     methods_objects["rm_rdisc"] = spectral.remove_rain_norain_discontinuity
diff --git a/pysteps/utils/reprojection.py b/pysteps/utils/reprojection.py
new file mode 100644
index 000000000..f9bf769aa
--- /dev/null
+++ b/pysteps/utils/reprojection.py
@@ -0,0 +1,120 @@
+# -*- coding: utf-8 -*-
+"""
+pysteps.utils.reprojection
+==========================
+
+Reprojection tools to reproject grids and adjust the grid cell size of an
+input field to a destination field.
+
+.. autosummary::
+    :toctree: ../generated/
+
+    reproject_grids
+"""
+from pysteps.exceptions import MissingOptionalDependency
+
+import numpy as np
+
+try:
+    from rasterio import Affine as A
+    from rasterio.warp import reproject, Resampling
+
+    RASTERIO_IMPORTED = True
+except ImportError:
+    RASTERIO_IMPORTED = False
+
+
+def reproject_grids(src_array, dst_array, metadata_src, metadata_dst):
+    """
+    Reproject precipitation fields to the domain of another precipitation field.
+
+    Parameters
+    ----------
+    src_array: array-like
+        Three-dimensional array of shape (t, x, y) containing a time series of
+        precipitation fields. These precipitation fields will be reprojected.
+    dst_array: array-like
+        Array containing a precipitation field or a time series of precipitation
+        fields. The xarray src_array will be reprojected to the domain of
+        dst_array.
+    metadata_src: dict
+        Metadata dictionary containing the projection, x- and ypixelsize, x1 and
+        y2 attributes of the src_array as described in the documentation of
+        :py:mod:`pysteps.io.importers`.
+    metadata_dst: dict
+        Metadata dictionary containing the projection, x- and ypixelsize, x1 and
+        y2 attributes of the dst_array.
+
+    Returns
+    -------
+    r_rprj: array-like
+        Three-dimensional array of shape (t, x, y) containing the precipitation
+        fields of src_array, but reprojected to the domain of dst_array.
+    metadata: dict
+        Metadata dictionary containing the projection, x- and ypixelsize, x1 and
+        y2 attributes of the reprojected src_array.
+    """
+
+    if not RASTERIO_IMPORTED:
+        raise MissingOptionalDependency(
+            "rasterio package is required for the reprojection module, but it is "
+            "not installed"
+        )
+
+    # Extract the grid info from src_array
+    src_crs = metadata_src["projection"]
+    x1_src = metadata_src["x1"]
+    y2_src = metadata_src["y2"]
+    xpixelsize_src = metadata_src["xpixelsize"]
+    ypixelsize_src = metadata_src["ypixelsize"]
+    src_transform = A.translation(float(x1_src), float(y2_src)) * A.scale(
+        float(xpixelsize_src), float(-ypixelsize_src)
+    )
+
+    # Extract the grid info from dst_array
+    dst_crs = metadata_dst["projection"]
+    x1_dst = metadata_dst["x1"]
+    y2_dst = metadata_dst["y2"]
+    xpixelsize_dst = metadata_dst["xpixelsize"]
+    ypixelsize_dst = metadata_dst["ypixelsize"]
+    dst_transform = A.translation(float(x1_dst), float(y2_dst)) * A.scale(
+        float(xpixelsize_dst), float(-ypixelsize_dst)
+    )
+
+    # Initialise the reprojected array
+    r_rprj = np.zeros((src_array.shape[0], dst_array.shape[-2], dst_array.shape[-1]))
+
+    # For every timestep, reproject the precipitation field of src_array to
+    # the domain of dst_array
+    if metadata_src["yorigin"] != metadata_dst["yorigin"]:
+        src_array = src_array[:, ::-1, :]
+
+    for i in range(src_array.shape[0]):
+        reproject(
+            src_array[i, :, :],
+            r_rprj[i, :, :],
+            src_transform=src_transform,
+            src_crs=src_crs,
+            dst_transform=dst_transform,
+            dst_crs=dst_crs,
+            resampling=Resampling.nearest,
+            dst_nodata=np.nan,
+        )
+
+    # Update the metadata
+    metadata = metadata_src.copy()
+
+    for key in [
+        "projection",
+        "yorigin",
+        "xpixelsize",
+        "ypixelsize",
+        "x1",
+        "x2",
+        "y1",
+        "y2",
+        "cartesian_unit",
+    ]:
+        metadata[key] = metadata_dst[key]
+
+    return r_rprj, metadata
diff --git a/requirements_dev.txt b/requirements_dev.txt
index 3f18ec557..a200b7834 100644
--- a/requirements_dev.txt
+++ b/requirements_dev.txt
@@ -17,6 +17,7 @@ cartopy>=0.18
 h5py
 scikit-image
 pandas
+rasterio
 
 # Testing
 pytest
diff --git a/tox.ini b/tox.ini
index e3a46cdd8..f37444840 100644
--- a/tox.ini
+++ b/tox.ini
@@ -36,6 +36,7 @@ conda_deps =
     pyproj
     cartopy
     pygrib
+	rasterio
 conda_channels = conda-forge
 
 setenv =


Test patch:
diff --git a/.github/workflows/test_pysteps.yml b/.github/workflows/test_pysteps.yml
index da7d95da2..35f229317 100644
--- a/.github/workflows/test_pysteps.yml
+++ b/.github/workflows/test_pysteps.yml
@@ -3,9 +3,13 @@ name: Test pysteps
 on:
   # Triggers the workflow on push or pull request events to the master branch
   push:
-    branches: [ master ]
+    branches:
+      - master
+      - pysteps-v2
   pull_request:
-    branches: [ master ]
+      branches:
+        - master
+        - pysteps-v2
 
 jobs:
   unit_tests:
@@ -32,6 +36,7 @@ jobs:
       - name: Install mamba and create environment
         uses: mamba-org/provision-with-micromamba@main
         with:
+          micromamba-version: 0.20.0
           environment-file: ci/ci_test_env.yml
           environment-name: test_environment
           extra-specs: python=${{ matrix.python-version }}
diff --git a/ci/ci_test_env.yml b/ci/ci_test_env.yml
index 2866ee1fe..b999f1b83 100644
--- a/ci/ci_test_env.yml
+++ b/ci/ci_test_env.yml
@@ -26,6 +26,7 @@ dependencies:
   - PyWavelets
   - pandas
   - scikit-image
+  - rasterio
   # Test dependencies
   - pytest
   - pytest-cov
diff --git a/pysteps/tests/helpers.py b/pysteps/tests/helpers.py
index 930a04565..151157243 100644
--- a/pysteps/tests/helpers.py
+++ b/pysteps/tests/helpers.py
@@ -77,10 +77,6 @@ def get_precipitation_fields(
     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.
-
     upscale: float or None, optional
         Upscale fields in space during the pre-processing steps.
         If it is None, the precipitation field is not modified.
@@ -93,6 +89,10 @@ def get_precipitation_fields(
     log_transform: bool
         Whether to transform the output to dB.
 
+    clip: scalars (left, right, bottom, top), optional
+        The extent of the bounding box in data coordinates to be used to clip
+        the data.
+
     Other Parameters
     ----------------
 
@@ -101,9 +101,7 @@ def get_precipitation_fields(
 
     Returns
     -------
-    reference_field : array
-
-    metadata : dict
+    data_array : xr.DataArray
     """
 
     if source == "bom":
diff --git a/pysteps/tests/test_blending_clim.py b/pysteps/tests/test_blending_clim.py
new file mode 100644
index 000000000..4a85a5c1b
--- /dev/null
+++ b/pysteps/tests/test_blending_clim.py
@@ -0,0 +1,150 @@
+# -*- coding: utf-8 -*-
+
+
+from datetime import datetime, timedelta
+from os.path import join, exists
+import pickle
+import random
+
+import numpy as np
+from numpy.testing import assert_array_equal
+import pytest
+
+from pysteps.blending.clim import save_skill, calc_clim_skill
+
+
+random.seed(12356)
+n_cascade_levels = 7
+model_names = ["alaro13", "arome13"]
+default_start_skill = [0.8, 0.5]
+
+# Helper functions
+def generate_fixed_skill(n_cascade_levels, n_models=1):
+    """
+    Generate skill starting at default_start_skill which decay exponentially with scale.
+    """
+    start_skill = np.resize(default_start_skill, n_models)
+    powers = np.arange(1, n_cascade_levels + 1)
+    return pow(start_skill[:, np.newaxis], powers)
+
+
+# Test arguments
+clim_arg_names = ("startdatestr", "enddatestr", "n_models", "expected_skill_today")
+
+test_enddates = ["20210701235500", "20210702000000", "20200930235500"]
+
+clim_arg_values = [
+    (
+        "20210701230000",
+        "20210701235500",
+        1,
+        {
+            "mean_skill": generate_fixed_skill(n_cascade_levels),
+            "n": 12,
+            "last_validtime": datetime.strptime(test_enddates[0], "%Y%m%d%H%M%S"),
+        },
+    ),
+    (
+        "20210701235500",
+        "20210702000000",
+        1,
+        {
+            "mean_skill": generate_fixed_skill(n_cascade_levels),
+            "n": 1,
+            "last_validtime": datetime.strptime(test_enddates[1], "%Y%m%d%H%M%S"),
+        },
+    ),
+    (
+        "20200801000000",
+        "20200930235500",
+        1,
+        {
+            "mean_skill": generate_fixed_skill(n_cascade_levels),
+            "n": 288,
+            "last_validtime": datetime.strptime(test_enddates[2], "%Y%m%d%H%M%S"),
+        },
+    ),
+    (
+        "20210701230000",
+        "20210701235500",
+        2,
+        {
+            "mean_skill": generate_fixed_skill(n_cascade_levels, 2),
+            "n": 12,
+            "last_validtime": datetime.strptime(test_enddates[0], "%Y%m%d%H%M%S"),
+        },
+    ),
+    (
+        "20210701230000",
+        "20210702000000",
+        2,
+        {
+            "mean_skill": generate_fixed_skill(n_cascade_levels, 2),
+            "n": 1,
+            "last_validtime": datetime.strptime(test_enddates[1], "%Y%m%d%H%M%S"),
+        },
+    ),
+    (
+        "20200801000000",
+        "20200930235500",
+        2,
+        {
+            "mean_skill": generate_fixed_skill(n_cascade_levels, 2),
+            "n": 288,
+            "last_validtime": datetime.strptime(test_enddates[2], "%Y%m%d%H%M%S"),
+        },
+    ),
+]
+
+
+@pytest.mark.parametrize(clim_arg_names, clim_arg_values)
+def test_save_skill(startdatestr, enddatestr, n_models, expected_skill_today, tmpdir):
+    """Test if the skill are saved correctly and the daily average is computed"""
+
+    # get validtime
+    currentdate = datetime.strptime(startdatestr, "%Y%m%d%H%M%S")
+    enddate = datetime.strptime(enddatestr, "%Y%m%d%H%M%S")
+    timestep = timedelta(minutes=5)
+
+    outdir_path = tmpdir
+
+    while currentdate <= enddate:
+        current_skill = generate_fixed_skill(n_cascade_levels, n_models)
+        print("Saving skill: ", current_skill, currentdate, outdir_path)
+        save_skill(
+            current_skill, currentdate, outdir_path, n_models=n_models, window_length=2
+        )
+        currentdate += timestep
+
+    skill_today_file = join(outdir_path, "NWP_skill_today.pkl")
+    assert exists(skill_today_file)
+    with open(skill_today_file, "rb") as f:
+        skill_today = pickle.load(f)
+
+    # Check type
+    assert isinstance(skill_today, dict)
+    assert "mean_skill" in skill_today
+    assert "n" in skill_today
+    assert "last_validtime" in skill_today
+    assert_array_equal(skill_today["mean_skill"], expected_skill_today["mean_skill"])
+    assert skill_today["n"] == expected_skill_today["n"]
+    assert skill_today["last_validtime"] == expected_skill_today["last_validtime"]
+
+    # Finally, check if the clim skill calculation returns an array of values
+    clim_skill = calc_clim_skill(
+        outdir_path=tmpdir,
+        n_cascade_levels=n_cascade_levels,
+        n_models=n_models,
+        window_length=2,
+    )
+
+    assert clim_skill.shape[0] == n_models
+    assert clim_skill.shape[1] == n_cascade_levels
+
+
+if __name__ == "__main__":
+    save_skill(
+        generate_fixed_skill(n_cascade_levels, 1),
+        datetime.strptime("20200801000000", "%Y%m%d%H%M%S"),
+        "./tmp/",
+    )
diff --git a/pysteps/tests/test_blending_linear_blending.py b/pysteps/tests/test_blending_linear_blending.py
new file mode 100644
index 000000000..46e23dd97
--- /dev/null
+++ b/pysteps/tests/test_blending_linear_blending.py
@@ -0,0 +1,142 @@
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import pytest
+from pysteps.blending.linear_blending import forecast
+from numpy.testing import assert_array_almost_equal
+from pysteps.utils import transformation
+
+# Test function arguments
+linear_arg_values = [
+    (5, 30, 60, 20, 45, "eulerian", None, 1, "True"),
+    (5, 30, 60, 20, 45, "eulerian", None, 2, "False"),
+    (5, 30, 60, 20, 45, "eulerian", None, 0, "False"),
+    (4, 23, 33, 9, 28, "eulerian", None, 1, "False"),
+    (3, 18, 36, 13, 27, "eulerian", None, 1, "False"),
+    (7, 30, 68, 11, 49, "eulerian", None, 1, "False"),
+    (10, 100, 160, 25, 130, "eulerian", None, 1, "False"),
+    (6, 60, 180, 22, 120, "eulerian", None, 1, "False"),
+    (5, 100, 200, 40, 150, "eulerian", None, 1, "False"),
+    (5, 30, 60, 20, 45, "extrapolation", np.zeros((2, 200, 200)), 1, "False"),
+    (4, 23, 33, 9, 28, "extrapolation", np.zeros((2, 200, 200)), 1, "False"),
+    (3, 18, 36, 13, 27, "extrapolation", np.zeros((2, 200, 200)), 1, "False"),
+    (7, 30, 68, 11, 49, "extrapolation", np.zeros((2, 200, 200)), 1, "False"),
+    (10, 100, 160, 25, 130, "extrapolation", np.zeros((2, 200, 200)), 1, "False"),
+    (6, 60, 180, 22, 120, "extrapolation", np.zeros((2, 200, 200)), 1, "False"),
+    (5, 100, 200, 40, 150, "extrapolation", np.zeros((2, 200, 200)), 1, "False"),
+]
+
+
+@pytest.mark.parametrize(
+    "timestep, start_blending, end_blending, n_timesteps, controltime, nowcast_method, V, n_models, squeeze_nwp_array",
+    linear_arg_values,
+)
+def test_linear_blending(
+    timestep,
+    start_blending,
+    end_blending,
+    n_timesteps,
+    controltime,
+    nowcast_method,
+    V,
+    n_models,
+    squeeze_nwp_array,
+):
+    """Tests if the linear blending function is correct. For the nowcast data a precipitation field
+    which is constant over time is taken. One half of the field has no rain and the other half
+    has a set value. For the NWP data a similar field is taken, the only difference
+    being that now the other half of the field is zero. The blended field should have a
+    constant value over the entire field at the timestep right in the middle between the start
+    of the blending and the end of the blending. This assertion is checked to see if the
+    linear blending function works well."""
+
+    # The argument controltime gives the timestep at which the field is assumed to be
+    # entirely constant
+
+    # Assert that the control time step is in the range of the forecasted time steps
+    assert controltime <= (
+        n_timesteps * timestep
+    ), "Control time needs to be within reach of forecasts, controltime = {} and n_timesteps = {}".format(
+        controltime, n_timesteps
+    )
+
+    # Assert that the start time of the blending comes before the end time of the blending
+    assert (
+        start_blending < end_blending
+    ), "Start time of blending needs to be smaller than end time of blending"
+
+    # Assert that the control time is a multiple of the time step
+    assert (
+        not controltime % timestep
+    ), "Control time needs to be a multiple of the time step"
+
+    # Initialise dummy NWP data
+    if n_models == 0:
+        r_nwp = None
+    else:
+        r_nwp = np.zeros((n_models, n_timesteps, 200, 200))
+
+        for i in range(100):
+            r_nwp[:, :, i, :] = 11.0
+
+        if squeeze_nwp_array:
+            r_nwp = np.squeeze(r_nwp)
+
+    # Define nowcast input data
+    r_input = np.zeros((200, 200))
+
+    for i in range(100, 200):
+        r_input[i, :] = 11.0
+
+    # Transform from mm/h to dB
+    r_input, _ = transformation.dB_transform(
+        r_input, None, threshold=0.1, zerovalue=-15.0
+    )
+
+    # Calculate the blended field
+    r_blended = forecast(
+        r_input,
+        dict({"unit": "mm/h", "transform": "dB"}),
+        V,
+        n_timesteps,
+        timestep,
+        nowcast_method,
+        r_nwp,
+        dict({"unit": "mm/h", "transform": None}),
+        start_blending=start_blending,
+        end_blending=end_blending,
+    )
+
+    # Assert that the blended field has the expected dimension
+    if n_models > 1:
+        assert r_blended.shape == (
+            n_models,
+            n_timesteps,
+            200,
+            200,
+        ), "The shape of the blended array does not have the expected value. The shape is {}".format(
+            r_blended.shape
+        )
+    else:
+        assert r_blended.shape == (
+            n_timesteps,
+            200,
+            200,
+        ), "The shape of the blended array does not have the expected value. The shape is {}".format(
+            r_blended.shape
+        )
+
+    # Assert that the blended field at the control time step is equal to
+    # a constant field with the expected value.
+    if n_models > 1:
+        assert_array_almost_equal(
+            r_blended[0, controltime // timestep - 1],
+            np.ones((200, 200)) * 5.5,
+            err_msg="The blended array does not have the expected value",
+        )
+    elif n_models > 0:
+        assert_array_almost_equal(
+            r_blended[controltime // timestep - 1],
+            np.ones((200, 200)) * 5.5,
+            err_msg="The blended array does not have the expected value",
+        )
diff --git a/pysteps/tests/test_blending_skill_scores.py b/pysteps/tests/test_blending_skill_scores.py
new file mode 100644
index 000000000..b8b648c1b
--- /dev/null
+++ b/pysteps/tests/test_blending_skill_scores.py
@@ -0,0 +1,294 @@
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import pytest
+from numpy.testing import assert_array_almost_equal
+
+from pysteps.blending.skill_scores import (
+    spatial_correlation,
+    lt_dependent_cor_nwp,
+    lt_dependent_cor_extrapolation,
+    clim_regr_values,
+)
+
+
+# Set the climatological correlations values
+clim_cor_values_8lev = np.array(
+    [0.848, 0.537, 0.237, 0.065, 0.020, 0.0044, 0.0052, 0.0040]
+)
+clim_cor_values_6lev = np.array([0.848, 0.537, 0.237, 0.065, 0.020, 0.0044])
+clim_cor_values_9lev = np.array(
+    [0.848, 0.537, 0.237, 0.065, 0.020, 0.0044, 0.0052, 0.0040, 1e-4]
+)
+
+# Set the regression values
+regr_pars_8lev = np.array(
+    [
+        [130.0, 165.0, 120.0, 55.0, 50.0, 15.0, 15.0, 10.0],
+        [155.0, 220.0, 200.0, 75.0, 10e4, 10e4, 10e4, 10e4],
+    ]
+)
+regr_pars_6lev = np.array(
+    [
+        [130.0, 165.0, 120.0, 55.0, 50.0, 15.0],
+        [155.0, 220.0, 200.0, 75.0, 10e4, 10e4],
+    ]
+)
+regr_pars_9lev = np.array(
+    [
+        [130.0, 165.0, 120.0, 55.0, 50.0, 15.0, 15.0, 10.0, 10.0],
+        [155.0, 220.0, 200.0, 75.0, 10e4, 10e4, 10e4, 10e4, 10e4],
+    ]
+)
+
+# Set the dummy observation and model values
+dummy_2d_array = np.array([[1.0, 2.0], [3.0, 4.0]])
+obs_8lev = np.repeat(dummy_2d_array[None, :, :], 8, axis=0)
+obs_6lev = np.repeat(dummy_2d_array[None, :, :], 6, axis=0)
+obs_9lev = np.repeat(dummy_2d_array[None, :, :], 9, axis=0)
+mod_8lev = np.repeat(dummy_2d_array[None, :, :], 8, axis=0)
+mod_6lev = np.repeat(dummy_2d_array[None, :, :], 6, axis=0)
+mod_9lev = np.repeat(dummy_2d_array[None, :, :], 9, axis=0)
+
+# Gives some dummy values to PHI
+dummy_phi = np.array([0.472650, 0.523825, 0.103454])
+PHI_8lev = np.repeat(dummy_phi[None, :], 8, axis=0)
+PHI_6lev = np.repeat(dummy_phi[None, :], 6, axis=0)
+PHI_9lev = np.repeat(dummy_phi[None, :], 9, axis=0)
+
+# Test function arguments
+skill_scores_arg_names = (
+    "obs",
+    "mod",
+    "lt",
+    "PHI",
+    "cor_prev",
+    "clim_cor_values",
+    "regr_pars",
+    "n_cascade_levels",
+    "expected_cor_t0",
+    "expected_cor_nwp_lt",
+    "expected_cor_nowcast_lt",
+)
+
+# Test function values
+skill_scores_arg_values = [
+    (
+        obs_8lev,
+        mod_8lev,
+        60,
+        PHI_8lev,
+        None,
+        clim_cor_values_8lev,
+        regr_pars_8lev,
+        8,
+        np.repeat(1.0, 8),
+        np.array(
+            [
+                0.97455941,
+                0.9356775,
+                0.81972779,
+                0.55202975,
+                0.31534738,
+                0.02264599,
+                0.02343133,
+                0.00647032,
+            ]
+        ),
+        np.array(
+            [
+                0.996475,
+                0.996475,
+                0.996475,
+                0.996475,
+                0.996475,
+                0.996475,
+                0.996475,
+                0.996475,
+            ]
+        ),
+    ),
+    (
+        obs_6lev,
+        mod_6lev,
+        60,
+        PHI_6lev,
+        None,
+        clim_cor_values_6lev,
+        regr_pars_6lev,
+        6,
+        np.repeat(1.0, 6),
+        np.array(
+            [0.97455941, 0.9356775, 0.81972779, 0.55202975, 0.31534738, 0.02264599]
+        ),
+        np.array([0.996475, 0.996475, 0.996475, 0.996475, 0.996475, 0.996475]),
+    ),
+    (
+        obs_9lev,
+        mod_9lev,
+        60,
+        PHI_9lev,
+        None,
+        clim_cor_values_9lev,
+        regr_pars_9lev,
+        9,
+        np.repeat(1.0, 9),
+        np.array(
+            [
+                0.97455941,
+                0.9356775,
+                0.81972779,
+                0.55202975,
+                0.31534738,
+                0.02264599,
+                0.02343133,
+                0.00647032,
+                0.00347776,
+            ]
+        ),
+        np.array(
+            [
+                0.996475,
+                0.996475,
+                0.996475,
+                0.996475,
+                0.996475,
+                0.996475,
+                0.996475,
+                0.996475,
+                0.996475,
+            ]
+        ),
+    ),
+    (
+        obs_8lev,
+        mod_8lev,
+        0,
+        PHI_8lev,
+        None,
+        clim_cor_values_8lev,
+        regr_pars_8lev,
+        8,
+        np.repeat(1.0, 8),
+        np.repeat(1.0, 8),
+        np.array(
+            [
+                0.996475,
+                0.996475,
+                0.996475,
+                0.996475,
+                0.996475,
+                0.996475,
+                0.996475,
+                0.996475,
+            ]
+        ),
+    ),
+]
+
+
+# The test
+@pytest.mark.parametrize(skill_scores_arg_names, skill_scores_arg_values)
+
+# The test function to be used
+def test_blending_skill_scores(
+    obs,
+    mod,
+    lt,
+    PHI,
+    cor_prev,
+    clim_cor_values,
+    regr_pars,
+    n_cascade_levels,
+    expected_cor_t0,
+    expected_cor_nwp_lt,
+    expected_cor_nowcast_lt,
+):
+    """Tests if the skill_score functions behave correctly. A dummy gridded
+    model and observation field should be given for n_cascade_levels, which
+    leads to a given spatial correlation per cascade level. Then, the function
+    tests if the correlation regresses towards the climatological values given
+    lead time lt for the NWP fields or given the PHI-values for the
+    extrapolation field.
+
+    """
+    domain_mask = np.full(obs[0, :, :].shape, False, dtype=bool)
+
+    # Calculate the spatial correlation of the given model field
+    correlations_t0 = np.array(spatial_correlation(obs, mod, domain_mask))
+
+    # Check if the field has the same number of cascade levels as the model
+    # field and as the given n_cascade_levels
+    assert (
+        correlations_t0.shape[0] == mod.shape[0]
+    ), "Number of cascade levels should be the same as in the model field"
+    assert (
+        correlations_t0.shape[0] == n_cascade_levels
+    ), "Number of cascade levels should be the same as n_cascade_levels"
+
+    # Check if the returned values are as expected
+    assert_array_almost_equal(
+        correlations_t0,
+        expected_cor_t0,
+        decimal=3,
+        err_msg="Returned spatial correlation is not the same as the expected value",
+    )
+
+    # Test if the NWP correlation regresses towards the correct value given
+    # a lead time in minutes
+    # First, check if the climatological values are returned correctly
+    correlations_clim, regr_clim = clim_regr_values(
+        n_cascade_levels=n_cascade_levels, outdir_path="./tmp/"
+    )
+    assert (
+        correlations_clim.shape[0] == n_cascade_levels
+    ), "Number of cascade levels should be the same as n_cascade_levels"
+    assert_array_almost_equal(
+        correlations_clim,
+        clim_cor_values,
+        decimal=3,
+        err_msg="Not the correct climatological correlations were returned",
+    )
+    assert_array_almost_equal(
+        regr_clim,
+        regr_pars,
+        decimal=3,
+        err_msg="Not the correct regression parameters were returned",
+    )
+
+    # Then, check the regression of the correlation values
+    correlations_nwp_lt = lt_dependent_cor_nwp(
+        lt=lt, correlations=correlations_t0, outdir_path="./tmp/"
+    )
+    assert (
+        correlations_nwp_lt.shape[0] == mod.shape[0]
+    ), "Number of cascade levels should be the same as in the model field"
+    assert (
+        correlations_nwp_lt.shape[0] == n_cascade_levels
+    ), "Number of cascade levels should be the same as n_cascade_levels"
+    assert_array_almost_equal(
+        correlations_nwp_lt,
+        expected_cor_nwp_lt,
+        decimal=3,
+        err_msg="Correlations of NWP not equal to the expected correlations",
+    )
+
+    # Finally, make sure nowcast correlation regresses towards the correct
+    # value given some PHI-values.
+    correlations_nowcast_lt, __ = lt_dependent_cor_extrapolation(
+        PHI, correlations_t0, cor_prev
+    )
+
+    print(correlations_nowcast_lt)
+    assert (
+        correlations_nowcast_lt.shape[0] == mod.shape[0]
+    ), "Number of cascade levels should be the same as in the model field"
+    assert (
+        correlations_nowcast_lt.shape[0] == n_cascade_levels
+    ), "Number of cascade levels should be the same as n_cascade_levels"
+    assert_array_almost_equal(
+        correlations_nowcast_lt,
+        expected_cor_nowcast_lt,
+        decimal=3,
+        err_msg="Correlations of nowcast not equal to the expected correlations",
+    )
diff --git a/pysteps/tests/test_blending_steps.py b/pysteps/tests/test_blending_steps.py
new file mode 100644
index 000000000..201fb04b3
--- /dev/null
+++ b/pysteps/tests/test_blending_steps.py
@@ -0,0 +1,271 @@
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import datetime
+import pytest
+import pysteps
+from pysteps import cascade, blending
+
+
+steps_arg_names = (
+    "n_models",
+    "n_timesteps",
+    "n_ens_members",
+    "n_cascade_levels",
+    "mask_method",
+    "probmatching_method",
+    "blend_nwp_members",
+    "weights_method",
+    "decomposed_nwp",
+    "expected_n_ens_members",
+)
+
+steps_arg_values = [
+    (1, 3, 4, 8, None, None, False, "spn", True, 4),
+    (1, 3, 4, 8, "obs", None, False, "spn", True, 4),
+    (1, 3, 4, 8, "incremental", None, False, "spn", True, 4),
+    (1, 3, 4, 8, None, "mean", False, "spn", True, 4),
+    (1, 3, 4, 8, None, "cdf", False, "spn", True, 4),
+    (1, 3, 4, 8, "incremental", "cdf", False, "spn", True, 4),
+    (1, 3, 4, 6, "incremental", "cdf", False, "bps", True, 4),
+    (1, 3, 4, 6, "incremental", "cdf", False, "bps", False, 4),
+    (1, 3, 4, 9, "incremental", "cdf", False, "spn", True, 4),
+    (2, 3, 10, 8, "incremental", "cdf", False, "spn", True, 10),
+    (5, 3, 4, 8, "incremental", "cdf", False, "spn", True, 5),
+    (1, 10, 1, 8, "incremental", "cdf", False, "spn", True, 1),
+    (5, 3, 2, 8, "incremental", "cdf", True, "spn", True, 2),
+]
+
+
+@pytest.mark.parametrize(steps_arg_names, steps_arg_values)
+def test_steps_blending(
+    n_models,
+    n_timesteps,
+    n_ens_members,
+    n_cascade_levels,
+    mask_method,
+    probmatching_method,
+    blend_nwp_members,
+    weights_method,
+    decomposed_nwp,
+    expected_n_ens_members,
+):
+    ###
+    # The input data
+    ###
+    # Initialise dummy NWP data
+    nwp_precip = np.zeros((n_models, n_timesteps + 1, 200, 200))
+
+    for n_model in range(n_models):
+        for i in range(nwp_precip.shape[1]):
+            nwp_precip[n_model, i, 30:185, 30 + 1 * (i + 1) * n_model] = 0.1
+            nwp_precip[n_model, i, 30:185, 31 + 1 * (i + 1) * n_model] = 0.1
+            nwp_precip[n_model, i, 30:185, 32 + 1 * (i + 1) * n_model] = 1.0
+            nwp_precip[n_model, i, 30:185, 33 + 1 * (i + 1) * n_model] = 5.0
+            nwp_precip[n_model, i, 30:185, 34 + 1 * (i + 1) * n_model] = 5.0
+            nwp_precip[n_model, i, 30:185, 35 + 1 * (i + 1) * n_model] = 4.5
+            nwp_precip[n_model, i, 30:185, 36 + 1 * (i + 1) * n_model] = 4.5
+            nwp_precip[n_model, i, 30:185, 37 + 1 * (i + 1) * n_model] = 4.0
+            nwp_precip[n_model, i, 30:185, 38 + 1 * (i + 1) * n_model] = 2.0
+            nwp_precip[n_model, i, 30:185, 39 + 1 * (i + 1) * n_model] = 1.0
+            nwp_precip[n_model, i, 30:185, 40 + 1 * (i + 1) * n_model] = 0.5
+            nwp_precip[n_model, i, 30:185, 41 + 1 * (i + 1) * n_model] = 0.1
+
+    # Define dummy nowcast input data
+    radar_precip = np.zeros((3, 200, 200))
+
+    for i in range(2):
+        radar_precip[i, 5:150, 30 + 1 * i] = 0.1
+        radar_precip[i, 5:150, 31 + 1 * i] = 0.5
+        radar_precip[i, 5:150, 32 + 1 * i] = 0.5
+        radar_precip[i, 5:150, 33 + 1 * i] = 5.0
+        radar_precip[i, 5:150, 34 + 1 * i] = 5.0
+        radar_precip[i, 5:150, 35 + 1 * i] = 4.5
+        radar_precip[i, 5:150, 36 + 1 * i] = 4.5
+        radar_precip[i, 5:150, 37 + 1 * i] = 4.0
+        radar_precip[i, 5:150, 38 + 1 * i] = 1.0
+        radar_precip[i, 5:150, 39 + 1 * i] = 0.5
+        radar_precip[i, 5:150, 40 + 1 * i] = 0.5
+        radar_precip[i, 5:150, 41 + 1 * i] = 0.1
+    radar_precip[2, 30:155, 30 + 1 * 2] = 0.1
+    radar_precip[2, 30:155, 31 + 1 * 2] = 0.1
+    radar_precip[2, 30:155, 32 + 1 * 2] = 1.0
+    radar_precip[2, 30:155, 33 + 1 * 2] = 5.0
+    radar_precip[2, 30:155, 34 + 1 * 2] = 5.0
+    radar_precip[2, 30:155, 35 + 1 * 2] = 4.5
+    radar_precip[2, 30:155, 36 + 1 * 2] = 4.5
+    radar_precip[2, 30:155, 37 + 1 * 2] = 4.0
+    radar_precip[2, 30:155, 38 + 1 * 2] = 2.0
+    radar_precip[2, 30:155, 39 + 1 * 2] = 1.0
+    radar_precip[2, 30:155, 40 + 1 * 3] = 0.5
+    radar_precip[2, 30:155, 41 + 1 * 3] = 0.1
+
+    metadata = dict()
+    metadata["unit"] = "mm"
+    metadata["transformation"] = "dB"
+    metadata["accutime"] = 5.0
+    metadata["transform"] = "dB"
+    metadata["zerovalue"] = 0.0
+    metadata["threshold"] = 0.01
+    metadata["zr_a"] = 200.0
+    metadata["zr_b"] = 1.6
+
+    # Also set the outdir_path and clim_kwargs
+    outdir_path_skill = "./tmp/"
+    clim_kwargs = dict({"n_models": n_models, "window_length": 30})
+
+    ###
+    # First threshold the data and convert it to dBR
+    ###
+    # threshold the data
+    radar_precip[radar_precip < metadata["threshold"]] = 0.0
+    nwp_precip[nwp_precip < metadata["threshold"]] = 0.0
+
+    # convert the data
+    converter = pysteps.utils.get_method("mm/h")
+    radar_precip, _ = converter(radar_precip, metadata)
+    nwp_precip, metadata = converter(nwp_precip, metadata)
+
+    # transform the data
+    transformer = pysteps.utils.get_method(metadata["transformation"])
+    radar_precip, _ = transformer(radar_precip, metadata)
+    nwp_precip, metadata = transformer(nwp_precip, metadata)
+
+    # set NaN equal to zero
+    radar_precip[~np.isfinite(radar_precip)] = metadata["zerovalue"]
+    nwp_precip[~np.isfinite(nwp_precip)] = metadata["zerovalue"]
+
+    assert (
+        np.any(~np.isfinite(radar_precip)) == False
+    ), "There are still infinite values in the input radar data"
+    assert (
+        np.any(~np.isfinite(nwp_precip)) == False
+    ), "There are still infinite values in the NWP data"
+
+    ###
+    # Decompose the R_NWP data
+    ###
+
+    # Initial decomposition settings
+    decomp_method, _ = cascade.get_method("fft")
+    bandpass_filter_method = "gaussian"
+    precip_shape = radar_precip.shape[1:]
+    filter_method = cascade.get_method(bandpass_filter_method)
+    bp_filter = filter_method(precip_shape, n_cascade_levels)
+
+    # If we only use one model:
+    if nwp_precip.ndim == 3:
+        nwp_precip = nwp_precip[None, :]
+
+    if decomposed_nwp:
+        nwp_precip_decomp = []
+        # Loop through the n_models
+        for i in range(nwp_precip.shape[0]):
+            R_d_models_ = []
+            # Loop through the time steps
+            for j in range(nwp_precip.shape[1]):
+                R_ = decomp_method(
+                    field=nwp_precip[i, j, :, :],
+                    bp_filter=bp_filter,
+                    normalize=True,
+                    compute_stats=True,
+                    compact_output=True,
+                )
+                R_d_models_.append(R_)
+            nwp_precip_decomp.append(R_d_models_)
+
+        nwp_precip_decomp = np.array(nwp_precip_decomp)
+
+        assert nwp_precip_decomp.ndim == 2, "Wrong number of dimensions in R_d_models"
+
+    else:
+        nwp_precip_decomp = nwp_precip.copy()
+
+        assert nwp_precip_decomp.ndim == 4, "Wrong number of dimensions in R_d_models"
+
+    ###
+    # Determine the velocity fields
+    ###
+    oflow_method = pysteps.motion.get_method("lucaskanade")
+    radar_velocity = oflow_method(radar_precip)
+    nwp_velocity = []
+    # Loop through the models
+    for n_model in range(nwp_precip.shape[0]):
+        # Loop through the timesteps. We need two images to construct a motion
+        # field, so we can start from timestep 1. Timestep 0 will be the same
+        # as timestep 0.
+        _V_NWP_ = []
+        for t in range(1, nwp_precip.shape[1]):
+            V_NWP_ = oflow_method(nwp_precip[n_model, t - 1 : t + 1, :])
+            _V_NWP_.append(V_NWP_)
+            V_NWP_ = None
+        _V_NWP_ = np.insert(_V_NWP_, 0, _V_NWP_[0], axis=0)
+        nwp_velocity.append(_V_NWP_)
+
+    nwp_velocity = np.stack(nwp_velocity)
+
+    assert nwp_velocity.ndim == 5, "nwp_velocity must be a five-dimensional array"
+
+    ###
+    # The nowcasting
+    ###
+    precip_forecast = blending.steps.forecast(
+        precip=radar_precip,
+        precip_models=nwp_precip_decomp,
+        velocity=radar_velocity,
+        velocity_models=nwp_velocity,
+        timesteps=n_timesteps,
+        timestep=5.0,
+        issuetime=datetime.datetime.strptime("202112012355", "%Y%m%d%H%M"),
+        n_ens_members=n_ens_members,
+        n_cascade_levels=n_cascade_levels,
+        blend_nwp_members=blend_nwp_members,
+        precip_thr=metadata["threshold"],
+        kmperpixel=1.0,
+        extrap_method="semilagrangian",
+        decomp_method="fft",
+        bandpass_filter_method="gaussian",
+        noise_method="nonparametric",
+        noise_stddev_adj="auto",
+        ar_order=2,
+        vel_pert_method=None,
+        weights_method=weights_method,
+        conditional=False,
+        probmatching_method=probmatching_method,
+        mask_method=mask_method,
+        callback=None,
+        return_output=True,
+        seed=None,
+        num_workers=1,
+        fft_method="numpy",
+        domain="spatial",
+        outdir_path_skill=outdir_path_skill,
+        extrap_kwargs=None,
+        filter_kwargs=None,
+        noise_kwargs=None,
+        vel_pert_kwargs=None,
+        clim_kwargs=clim_kwargs,
+        mask_kwargs=None,
+        measure_time=False,
+    )
+
+    assert precip_forecast.ndim == 4, "Wrong amount of dimensions in forecast output"
+    assert (
+        precip_forecast.shape[0] == expected_n_ens_members
+    ), "Wrong amount of output ensemble members in forecast output"
+    assert (
+        precip_forecast.shape[1] == n_timesteps
+    ), "Wrong amount of output time steps in forecast output"
+
+    # Transform the data back into mm/h
+    precip_forecast, _ = converter(precip_forecast, metadata)
+
+    assert (
+        precip_forecast.ndim == 4
+    ), "Wrong amount of dimensions in converted forecast output"
+    assert (
+        precip_forecast.shape[0] == expected_n_ens_members
+    ), "Wrong amount of output ensemble members in converted forecast output"
+    assert (
+        precip_forecast.shape[1] == n_timesteps
+    ), "Wrong amount of output time steps in converted forecast output"
diff --git a/pysteps/tests/test_blending_utils.py b/pysteps/tests/test_blending_utils.py
new file mode 100644
index 000000000..40f138d9c
--- /dev/null
+++ b/pysteps/tests/test_blending_utils.py
@@ -0,0 +1,380 @@
+# -*- coding: utf-8 -*-
+
+import os
+
+import numpy as np
+from numpy.testing import assert_array_almost_equal
+import pytest
+
+import pysteps
+from pysteps.blending.utils import (
+    stack_cascades,
+    blend_cascades,
+    recompose_cascade,
+    blend_optical_flows,
+    decompose_NWP,
+    compute_store_nwp_motion,
+    load_NWP,
+)
+
+pytest.importorskip("netCDF4")
+
+precip_nwp = np.zeros((24, 564, 564))
+
+for t in range(precip_nwp.shape[0]):
+    precip_nwp[t, 30 + t : 185 + t, 30 + 2 * t] = 0.1
+    precip_nwp[t, 30 + t : 185 + t, 31 + 2 * t] = 0.1
+    precip_nwp[t, 30 + t : 185 + t, 32 + 2 * t] = 1.0
+    precip_nwp[t, 30 + t : 185 + t, 33 + 2 * t] = 5.0
+    precip_nwp[t, 30 + t : 185 + t, 34 + 2 * t] = 5.0
+    precip_nwp[t, 30 + t : 185 + t, 35 + 2 * t] = 4.5
+    precip_nwp[t, 30 + t : 185 + t, 36 + 2 * t] = 4.5
+    precip_nwp[t, 30 + t : 185 + t, 37 + 2 * t] = 4.0
+    precip_nwp[t, 30 + t : 185 + t, 38 + 2 * t] = 2.0
+    precip_nwp[t, 30 + t : 185 + t, 39 + 2 * t] = 1.0
+    precip_nwp[t, 30 + t : 185 + t, 40 + 2 * t] = 0.5
+    precip_nwp[t, 30 + t : 185 + t, 41 + 2 * t] = 0.1
+
+nwp_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.000000001"
+)
+
+nwp_metadata = dict(
+    projection=nwp_proj,
+    institution="Royal Meteorological Institute of Belgium",
+    transform=None,
+    zerovalue=0.0,
+    threshold=0,
+    unit="mm",
+    accutime=5,
+    xpixelsize=1300.0,
+    ypixelsize=1300.0,
+    yorigin="upper",
+    cartesian_unit="m",
+    x1=0.0,
+    x2=731900.0,
+    y1=-731900.0,
+    y2=0.0,
+)
+
+# Get the analysis time and valid time
+times_nwp = np.array(
+    [
+        "2021-07-04T16:05:00.000000000",
+        "2021-07-04T16:10:00.000000000",
+        "2021-07-04T16:15:00.000000000",
+        "2021-07-04T16:20:00.000000000",
+        "2021-07-04T16:25:00.000000000",
+        "2021-07-04T16:30:00.000000000",
+        "2021-07-04T16:35:00.000000000",
+        "2021-07-04T16:40:00.000000000",
+        "2021-07-04T16:45:00.000000000",
+        "2021-07-04T16:50:00.000000000",
+        "2021-07-04T16:55:00.000000000",
+        "2021-07-04T17:00:00.000000000",
+        "2021-07-04T17:05:00.000000000",
+        "2021-07-04T17:10:00.000000000",
+        "2021-07-04T17:15:00.000000000",
+        "2021-07-04T17:20:00.000000000",
+        "2021-07-04T17:25:00.000000000",
+        "2021-07-04T17:30:00.000000000",
+        "2021-07-04T17:35:00.000000000",
+        "2021-07-04T17:40:00.000000000",
+        "2021-07-04T17:45:00.000000000",
+        "2021-07-04T17:50:00.000000000",
+        "2021-07-04T17:55:00.000000000",
+        "2021-07-04T18:00:00.000000000",
+    ],
+    dtype="datetime64[ns]",
+)
+
+
+# Prepare input NWP files
+# Convert to rain rates [mm/h]
+converter = pysteps.utils.get_method("mm/h")
+precip_nwp, nwp_metadata = converter(precip_nwp, nwp_metadata)
+
+# Threshold the data
+precip_nwp[precip_nwp < 0.1] = 0.0
+nwp_metadata["threshold"] = 0.1
+
+# Transform the data
+transformer = pysteps.utils.get_method("dB")
+precip_nwp, nwp_metadata = transformer(
+    precip_nwp, nwp_metadata, threshold=nwp_metadata["threshold"]
+)
+
+# Set two issue times for testing
+issue_time_first = times_nwp[0]
+issue_time_second = times_nwp[3]
+
+# Set the blending weights (we'll blend with a 50-50 weight)
+weights = np.full((2, 8), fill_value=0.5)
+
+# Set the testing arguments
+# Test function arguments
+utils_arg_names = (
+    "precip_nwp",
+    "nwp_model",
+    "issue_times",
+    "timestep",
+    "n_timesteps",
+    "valid_times",
+    "shape",
+    "weights",
+)
+
+# Test function values
+utils_arg_values = [
+    (
+        precip_nwp,
+        "test",
+        [issue_time_first, issue_time_second],
+        5.0,
+        3,
+        times_nwp,
+        precip_nwp.shape[1:],
+        weights,
+    )
+]
+
+
+###
+# The test
+###
+@pytest.mark.parametrize(utils_arg_names, utils_arg_values)
+
+# The test function to be used
+def test_blending_utils(
+    precip_nwp,
+    nwp_model,
+    issue_times,
+    timestep,
+    n_timesteps,
+    valid_times,
+    shape,
+    weights,
+):
+    """Tests if all blending utils functions behave correctly."""
+
+    # First, make the output path if it does not exist yet
+    tmpdir = "./tmp/"
+    os.makedirs(tmpdir, exist_ok=True)
+
+    # Get the optical flow method
+    oflow_method = pysteps.motion.get_method("lucaskanade")
+
+    ###
+    # Compute and store the motion
+    ###
+    compute_store_nwp_motion(
+        precip_nwp=precip_nwp,
+        oflow_method=oflow_method,
+        analysis_time=valid_times[0],
+        nwp_model=nwp_model,
+        output_path=tmpdir,
+    )
+
+    # Check if file exists
+    date_string = np.datetime_as_string(valid_times[0], "s")
+    motion_file = os.path.join(
+        tmpdir,
+        "motion_"
+        + nwp_model
+        + "_"
+        + date_string[:4]
+        + date_string[5:7]
+        + date_string[8:10]
+        + date_string[11:13]
+        + date_string[14:16]
+        + date_string[17:19]
+        + ".npy",
+    )
+    assert os.path.exists(motion_file)
+
+    ###
+    # Decompose and store NWP forecast
+    ###
+    decompose_NWP(
+        R_NWP=precip_nwp,
+        NWP_model=nwp_model,
+        analysis_time=valid_times[0],
+        timestep=timestep,
+        valid_times=valid_times,
+        num_cascade_levels=8,
+        num_workers=1,
+        output_path=tmpdir,
+        decomp_method="fft",
+        fft_method="numpy",
+        domain="spatial",
+        normalize=True,
+        compute_stats=True,
+        compact_output=False,
+    )
+
+    # Check if file exists
+    decomp_file = os.path.join(
+        tmpdir,
+        "cascade_"
+        + nwp_model
+        + "_"
+        + date_string[:4]
+        + date_string[5:7]
+        + date_string[8:10]
+        + date_string[11:13]
+        + date_string[14:16]
+        + date_string[17:19]
+        + ".nc",
+    )
+    assert os.path.exists(decomp_file)
+
+    ###
+    # Now check if files load correctly for two different issue times
+    ###
+    precip_decomposed_nwp_first, v_nwp_first = load_NWP(
+        input_nc_path_decomp=os.path.join(decomp_file),
+        input_path_velocities=os.path.join(motion_file),
+        start_time=issue_times[0],
+        n_timesteps=n_timesteps,
+    )
+
+    precip_decomposed_nwp_second, v_nwp_second = load_NWP(
+        input_nc_path_decomp=os.path.join(decomp_file),
+        input_path_velocities=os.path.join(motion_file),
+        start_time=issue_times[1],
+        n_timesteps=n_timesteps,
+    )
+
+    # Check if the output type and shapes are correct
+    assert isinstance(precip_decomposed_nwp_first, list)
+    assert isinstance(precip_decomposed_nwp_second, list)
+    assert isinstance(precip_decomposed_nwp_first[0], dict)
+    assert isinstance(precip_decomposed_nwp_second[0], dict)
+
+    assert "domain" in precip_decomposed_nwp_first[0]
+    assert "normalized" in precip_decomposed_nwp_first[0]
+    assert "compact_output" in precip_decomposed_nwp_first[0]
+    assert "valid_times" in precip_decomposed_nwp_first[0]
+    assert "cascade_levels" in precip_decomposed_nwp_first[0]
+    assert "means" in precip_decomposed_nwp_first[0]
+    assert "stds" in precip_decomposed_nwp_first[0]
+
+    assert precip_decomposed_nwp_first[0]["cascade_levels"].shape == (
+        8,
+        shape[0],
+        shape[1],
+    )
+    assert precip_decomposed_nwp_first[0]["domain"] == "spatial"
+    assert precip_decomposed_nwp_first[0]["normalized"] == True
+    assert precip_decomposed_nwp_first[0]["compact_output"] == False
+    assert len(precip_decomposed_nwp_first) == n_timesteps + 1
+    assert len(precip_decomposed_nwp_second) == n_timesteps + 1
+    assert precip_decomposed_nwp_first[0]["means"].shape[0] == 8
+    assert precip_decomposed_nwp_first[0]["stds"].shape[0] == 8
+
+    assert np.array(v_nwp_first).shape == (n_timesteps + 1, 2, shape[0], shape[1])
+    assert np.array(v_nwp_second).shape == (n_timesteps + 1, 2, shape[0], shape[1])
+
+    # Check if the right times are loaded
+    assert (
+        precip_decomposed_nwp_first[0]["valid_times"][0] == valid_times[0]
+    ), "Not the right valid times were loaded for the first forecast"
+    assert (
+        precip_decomposed_nwp_second[0]["valid_times"][0] == valid_times[3]
+    ), "Not the right valid times were loaded for the second forecast"
+
+    # Check, for a sample, if the stored motion fields are as expected
+    assert_array_almost_equal(
+        v_nwp_first[1],
+        oflow_method(precip_nwp[0:2, :, :]),
+        decimal=3,
+        err_msg="Stored motion field of first forecast not equal to expected motion field",
+    )
+    assert_array_almost_equal(
+        v_nwp_second[1],
+        oflow_method(precip_nwp[3:5, :, :]),
+        decimal=3,
+        err_msg="Stored motion field of second forecast not equal to expected motion field",
+    )
+
+    ###
+    # Stack the cascades
+    ###
+    precip_decomposed_first_stack, mu_first_stack, sigma_first_stack = stack_cascades(
+        R_d=precip_decomposed_nwp_first, donorm=False
+    )
+
+    print(precip_decomposed_nwp_first)
+    print(precip_decomposed_first_stack)
+    print(mu_first_stack)
+
+    (
+        precip_decomposed_second_stack,
+        mu_second_stack,
+        sigma_second_stack,
+    ) = stack_cascades(R_d=precip_decomposed_nwp_second, donorm=False)
+
+    # Check if the array shapes are still correct
+    assert precip_decomposed_first_stack.shape == (
+        n_timesteps + 1,
+        8,
+        shape[0],
+        shape[1],
+    )
+    assert mu_first_stack.shape == (n_timesteps + 1, 8)
+    assert sigma_first_stack.shape == (n_timesteps + 1, 8)
+
+    ###
+    # Blend the cascades
+    ###
+    precip_decomposed_blended = blend_cascades(
+        cascades_norm=np.stack(
+            (precip_decomposed_first_stack[0], precip_decomposed_second_stack[0])
+        ),
+        weights=weights,
+    )
+
+    assert precip_decomposed_blended.shape == precip_decomposed_first_stack[0].shape
+
+    ###
+    # Blend the optical flow fields
+    ###
+    v_nwp_blended = blend_optical_flows(
+        flows=np.stack((v_nwp_first[1], v_nwp_second[1])), weights=weights[:, 1]
+    )
+
+    assert v_nwp_blended.shape == v_nwp_first[1].shape
+    assert_array_almost_equal(
+        v_nwp_blended,
+        (oflow_method(precip_nwp[0:2, :, :]) + oflow_method(precip_nwp[3:5, :, :])) / 2,
+        decimal=3,
+        err_msg="Blended motion field does not equal average of the two motion fields",
+    )
+
+    ###
+    # Recompose the fields (the non-blended fields are used for this here)
+    ###
+    precip_recomposed_first = recompose_cascade(
+        combined_cascade=precip_decomposed_first_stack[0],
+        combined_mean=mu_first_stack[0],
+        combined_sigma=sigma_first_stack[0],
+    )
+    precip_recomposed_second = recompose_cascade(
+        combined_cascade=precip_decomposed_second_stack[0],
+        combined_mean=mu_second_stack[0],
+        combined_sigma=sigma_second_stack[0],
+    )
+
+    assert_array_almost_equal(
+        precip_recomposed_first,
+        precip_nwp[0, :, :],
+        decimal=3,
+        err_msg="Recomposed field of first forecast does not equal original field",
+    )
+    assert_array_almost_equal(
+        precip_recomposed_second,
+        precip_nwp[3, :, :],
+        decimal=3,
+        err_msg="Recomposed field of second forecast does not equal original field",
+    )
diff --git a/pysteps/tests/test_datasets.py b/pysteps/tests/test_datasets.py
index 93eefaaaa..72d6f3660 100644
--- a/pysteps/tests/test_datasets.py
+++ b/pysteps/tests/test_datasets.py
@@ -13,6 +13,7 @@
 )
 from pysteps.exceptions import DirectoryNotEmpty
 
+
 _datasets_opt_deps = dict(
     fmi=["pyproj"],
     mch=["PIL"],
diff --git a/pysteps/tests/test_interfaces.py b/pysteps/tests/test_interfaces.py
index 49b6bb69b..e2d499543 100644
--- a/pysteps/tests/test_interfaces.py
+++ b/pysteps/tests/test_interfaces.py
@@ -19,6 +19,25 @@ def _generic_interface_test(method_getter, valid_names_func_pair, invalid_names)
             method_getter(invalid_name)
 
 
+def test_nowcasts_interface():
+    """Test the nowcasts module interface."""
+
+    from pysteps.blending import (
+        linear_blending,
+        steps,
+    )
+
+    method_getter = pysteps.nowcasts.interface.get_method
+
+    valid_names_func_pair = [
+        ("linear_blending", linear_blending.forecast),
+        ("steps", steps.forecast),
+    ]
+
+    invalid_names = ["linear", "step", "blending", "pysteps"]
+    _generic_interface_test(method_getter, valid_names_func_pair, invalid_names)
+
+
 def test_cascade_interface():
     """Test the cascade module interface."""
 
@@ -251,6 +270,7 @@ def test_utils_interface():
     from pysteps.utils import dimension
     from pysteps.utils import images
     from pysteps.utils import interpolate
+    from pysteps.utils import reprojection
     from pysteps.utils import spectral
     from pysteps.utils import transformation
 
@@ -272,6 +292,7 @@ def test_utils_interface():
         ("upscale", dimension.aggregate_fields_space),
         ("morph_opening", images.morph_opening),
         ("rbfinterp2d", interpolate.rbfinterp2d),
+        ("reproject_grids", reprojection.reproject_grids),
         ("rapsd", spectral.rapsd),
         ("rm_rdisc", spectral.remove_rain_norain_discontinuity),
         ("boxcox", transformation.boxcox_transform),
diff --git a/pysteps/tests/test_io_archive.py b/pysteps/tests/test_io_archive.py
new file mode 100644
index 000000000..74a0567ec
--- /dev/null
+++ b/pysteps/tests/test_io_archive.py
@@ -0,0 +1,35 @@
+from datetime import datetime
+
+import pytest
+
+import pysteps
+
+
+def test_find_by_date_mch():
+
+    pytest.importorskip("PIL")
+
+    date = datetime.strptime("201505151630", "%Y%m%d%H%M")
+    data_source = pysteps.rcparams.data_sources["mch"]
+    root_path = data_source["root_path"]
+    path_fmt = data_source["path_fmt"]
+    fn_pattern = data_source["fn_pattern"]
+    fn_ext = data_source["fn_ext"]
+    timestep = data_source["timestep"]
+
+    fns = pysteps.io.archive.find_by_date(
+        date,
+        root_path,
+        path_fmt,
+        fn_pattern,
+        fn_ext,
+        timestep=timestep,
+        num_prev_files=1,
+        num_next_files=1,
+    )
+
+    assert len(fns) == 2
+    assert len(fns[0]) == 3
+    assert len(fns[1]) == 3
+    assert isinstance(fns[0][0], str)
+    assert isinstance(fns[1][0], datetime)
diff --git a/pysteps/tests/test_io_bom_rf3.py b/pysteps/tests/test_io_bom_rf3.py
index 6c2f17062..8d6a3cec5 100644
--- a/pysteps/tests/test_io_bom_rf3.py
+++ b/pysteps/tests/test_io_bom_rf3.py
@@ -7,23 +7,13 @@
 import pysteps
 from pysteps.tests.helpers import smart_assert
 
-pytest.importorskip("netCDF4")
-
-
-def test_io_import_bom_rf3_shape():
-    """Test the importer Bom RF3."""
-    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_100000.prcp-cscn.nc")
-    precip, _, _ = pysteps.io.import_bom_rf3(filename)
-    assert precip.shape == (512, 512)
-
+netCDF4 = pytest.importorskip("netCDF4")
 
+# Test import_bom_rf3 function
 expected_proj1 = (
     "+proj=aea  +lon_0=144.752 +lat_0=-37.852 " "+lat_1=-18.000 +lat_2=-36.000"
 )
 
-# test_metadata_bom: list of (variable,expected,tolerance) tuples
 test_metadata_bom = [
     ("transform", None, None),
     ("zerovalue", 0.0, 0.1),
@@ -48,10 +38,22 @@ def test_io_import_bom_rf3_metadata(variable, expected, tolerance):
     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_100000.prcp-cscn.nc")
-    _, _, metadata = pysteps.io.import_bom_rf3(filename)
+    precip, _, metadata = pysteps.io.import_bom_rf3(filename)
     smart_assert(metadata[variable], expected, tolerance)
+    assert precip.shape == (512, 512)
 
 
+# Test _import_bom_rf3_data function
+def test_io_import_bom_rf3_shape():
+    """Test the importer Bom RF3."""
+    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_100000.prcp-cscn.nc")
+    precip, _ = pysteps.io.importers._import_bom_rf3_data(filename)
+    assert precip.shape == (512, 512)
+
+
+# Test _import_bom_rf3_geodata function
 expected_proj2 = (
     "+proj=aea  +lon_0=144.752 +lat_0=-37.852 " "+lat_1=-18.000 +lat_2=-36.000"
 )
@@ -78,5 +80,8 @@ def test_io_import_bom_rf3_geodata(variable, expected, tolerance):
     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_100000.prcp-cscn.nc")
-    geodata = pysteps.io.importers._import_bom_rf3_geodata(filename)
+    ds_rainfall = netCDF4.Dataset(filename)
+    geodata = pysteps.io.importers._import_bom_rf3_geodata(ds_rainfall)
     smart_assert(geodata[variable], expected, tolerance)
+
+    ds_rainfall.close()
diff --git a/pysteps/tests/test_io_fmi_pgm.py b/pysteps/tests/test_io_fmi_pgm.py
index 58dfd5137..1658db957 100644
--- a/pysteps/tests/test_io_fmi_pgm.py
+++ b/pysteps/tests/test_io_fmi_pgm.py
@@ -10,62 +10,65 @@
 pytest.importorskip("pyproj")
 
 
-def test_io_import_fmi_pgm_shape():
-    """Test the importer FMI PGM."""
-    root_path = pysteps.rcparams.data_sources["fmi"]["root_path"]
-    filename = os.path.join(
-        root_path,
-        "20160928",
-        "201609281600_fmi.radar.composite.lowest_FIN_SUOMI1.pgm.gz",
-    )
-    R, _, _ = pysteps.io.import_fmi_pgm(filename, gzipped=True)
-    assert R.shape == (1226, 760)
-
-
-# test_metadata: list of (variable,expected) tuples
-test_metadata = [
-    ("composite_area", ["FIN"]),
-    ("projection_name", ["SUOMI1"]),
-    ("radar", ["LUO", "1", "26.9008", "67.1386"]),
-    ("obstime", ["201609281600"]),
-    ("producttype", ["CAPPI"]),
-    ("productname", ["LOWEST"]),
-    ("param", ["CorrectedReflectivity"]),
-    ("metersperpixel_x", ["999.674053"]),
-    ("metersperpixel_y", ["999.62859"]),
-    ("projection", ["radar", "{"]),
-    ("type", ["stereographic"]),
-    ("centrallongitude", ["25"]),
-    ("centrallatitude", ["90"]),
-    ("truelatitude", ["60"]),
-    ("bottomleft", ["18.600000", "57.930000"]),
-    ("topright", ["34.903000", "69.005000"]),
-    ("missingval", 255),
-]
+root_path = pysteps.rcparams.data_sources["fmi"]["root_path"]
+filename = os.path.join(
+    root_path,
+    "20160928",
+    "201609281600_fmi.radar.composite.lowest_FIN_SUOMI1.pgm.gz",
+)
+precip, _, metadata = pysteps.io.import_fmi_pgm(filename, gzipped=True)
 
 
-@pytest.mark.parametrize("variable,expected", test_metadata)
-def test_io_import_fmi_pmg_metadata(variable, expected):
-    """Test the importer FMI PMG."""
-    root_path = pysteps.rcparams.data_sources["fmi"]["root_path"]
-    filename = os.path.join(
-        root_path,
-        "20160928",
-        "201609281600_fmi.radar.composite.lowest_FIN_SUOMI1.pgm.gz",
-    )
-    metadata = pysteps.io.importers._import_fmi_pgm_metadata(filename, gzipped=True)
-    assert metadata[variable] == expected
+def test_io_import_fmi_pgm_shape():
+    """Test the importer FMI PGM."""
+    assert precip.shape == (1226, 760)
 
 
-expected_proj1 = (
+expected_proj = (
     "+proj=stere  +lon_0=25E +lat_0=90N "
     "+lat_ts=60 +a=6371288 +x_0=380886.310 "
     "+y_0=3395677.920 +no_defs"
 )
 
+test_attrs = [
+    ("projection", expected_proj, None),
+    ("institution", "Finnish Meteorological Institute", None),
+    # ("composite_area", ["FIN"]),
+    # ("projection_name", ["SUOMI1"]),
+    # ("radar", ["LUO", "1", "26.9008", "67.1386"]),
+    # ("obstime", ["201609281600"]),
+    # ("producttype", ["CAPPI"]),
+    # ("productname", ["LOWEST"]),
+    # ("param", ["CorrectedReflectivity"]),
+    # ("metersperpixel_x", ["999.674053"]),
+    # ("metersperpixel_y", ["999.62859"]),
+    # ("projection", ["radar", "{"]),
+    # ("type", ["stereographic"]),
+    # ("centrallongitude", ["25"]),
+    # ("centrallatitude", ["90"]),
+    # ("truelatitude", ["60"]),
+    # ("bottomleft", ["18.600000", "57.930000"]),
+    # ("topright", ["34.903000", "69.005000"]),
+    # ("missingval", 255),
+    ("accutime", 5.0, 0.1),
+    ("unit", "dBZ", None),
+    ("transform", "dB", None),
+    ("zerovalue", -32.0, 0.1),
+    ("threshold", -31.5, 0.1),
+    ("zr_a", 223.0, 0.1),
+    ("zr_b", 1.53, 0.1),
+]
+
+
+@pytest.mark.parametrize("variable, expected, tolerance", test_attrs)
+def test_io_import_mch_gif_dataset_attrs(variable, expected, tolerance):
+    """Test the importer FMI PMG."""
+    smart_assert(metadata[variable], expected, tolerance)
+
+
 # test_geodata: list of (variable,expected,tolerance) tuples
 test_geodata = [
-    ("projection", expected_proj1, None),
+    ("projection", expected_proj, None),
     ("x1", 0.0049823258887045085, 1e-20),
     ("x2", 759752.2852757066, 1e-10),
     ("y1", 0.009731985162943602, 1e-18),
diff --git a/pysteps/tests/test_io_knmi_hdf5.py b/pysteps/tests/test_io_knmi_hdf5.py
index 44e777a11..3e30cb575 100644
--- a/pysteps/tests/test_io_knmi_hdf5.py
+++ b/pysteps/tests/test_io_knmi_hdf5.py
@@ -10,35 +10,32 @@
 pytest.importorskip("h5py")
 
 
-def test_io_import_knmi_hdf5_shape():
-    """Test the importer KNMI HDF5."""
-
-    root_path = pysteps.rcparams.data_sources["knmi"]["root_path"]
+root_path = pysteps.rcparams.data_sources["knmi"]["root_path"]
+filename = os.path.join(root_path, "2010/08", "RAD_NL25_RAP_5min_201008260000.h5")
+precip, _, metadata = pysteps.io.import_knmi_hdf5(filename)
 
-    filename = os.path.join(root_path, "2010/08", "RAD_NL25_RAP_5min_201008260000.h5")
-
-    precip, __, __ = pysteps.io.import_knmi_hdf5(filename)
 
+def test_io_import_knmi_hdf5_shape():
+    """Test the importer KNMI HDF5."""
     assert precip.shape == (765, 700)
 
 
 # test_metadata: list of (variable,expected, tolerance) tuples
 
-expected_proj1 = (
-    "+proj=stere +lat_0=90 +lon_0=0.0 "
-    "+lat_ts=60.0 +a=6378.137 +b=6356.752 +x_0=0 "
-    "+y_0=0"
+expected_proj = (
+    "+proj=stere +lat_0=90 +lon_0=0.0 +lat_ts=60.0 +a=6378137 +b=6356752 +x_0=0 +y_0=0"
 )
 
-test_metadata = [
-    ("projection", expected_proj1, None),
+# list of (variable,expected,tolerance) tuples
+test_attrs = [
+    ("projection", expected_proj, None),
     ("x1", 0.0, 1e-10),
-    ("y1", -4415.038179210632, 1e-10),
-    ("x2", 699.9842646331593, 1e-10),
-    ("y2", -3649.950360247753, 1e-10),
-    ("xpixelsize", 1.0, 1e-10),
-    ("xpixelsize", 1.0, 1e-10),
-    ("cartesian_unit", "km", None),
+    ("y1", -4415038.179210632, 1e-10),
+    ("x2", 699984.2646331593, 1e-10),
+    ("y2", -3649950.360247753, 1e-10),
+    ("xpixelsize", 1000.0, 1e-10),
+    ("xpixelsize", 1000.0, 1e-10),
+    ("cartesian_unit", "m", None),
     ("accutime", 5.0, 1e-10),
     ("yorigin", "upper", None),
     ("unit", "mm", None),
@@ -51,13 +48,7 @@ def test_io_import_knmi_hdf5_shape():
 ]
 
 
-@pytest.mark.parametrize("variable,expected,tolerance", test_metadata)
+@pytest.mark.parametrize("variable,expected,tolerance", test_attrs)
 def test_io_import_knmi_hdf5_metadata(variable, expected, tolerance):
     """Test the importer KNMI HDF5."""
-    root_path = pysteps.rcparams.data_sources["knmi"]["root_path"]
-
-    filename = os.path.join(root_path, "2010/08", "RAD_NL25_RAP_5min_201008260000.h5")
-
-    __, __, metadata = pysteps.io.import_knmi_hdf5(filename)
-
     smart_assert(metadata[variable], expected, tolerance)
diff --git a/pysteps/tests/test_io_mch_gif.py b/pysteps/tests/test_io_mch_gif.py
index f569dd4a2..903033988 100644
--- a/pysteps/tests/test_io_mch_gif.py
+++ b/pysteps/tests/test_io_mch_gif.py
@@ -9,16 +9,17 @@
 
 pytest.importorskip("PIL")
 
+root_path = pysteps.rcparams.data_sources["mch"]["root_path"]
+filename = os.path.join(root_path, "20170131", "AQC170310945F_00005.801.gif")
+precip, _, metadata = pysteps.io.import_mch_gif(filename, "AQC", "mm", 5.0)
+
 
 def test_io_import_mch_gif_shape():
     """Test the importer MCH GIF."""
-    root_path = pysteps.rcparams.data_sources["mch"]["root_path"]
-    filename = os.path.join(root_path, "20170131", "AQC170310945F_00005.801.gif")
-    precip, _, metadata = pysteps.io.import_mch_gif(filename, "AQC", "mm", 5.0)
     assert precip.shape == (640, 710)
 
 
-expected_proj1 = (
+expected_proj = (
     "+proj=somerc  +lon_0=7.43958333333333 "
     "+lat_0=46.9524055555556 +k_0=1 "
     "+x_0=600000 +y_0=200000 +ellps=bessel "
@@ -26,9 +27,17 @@ def test_io_import_mch_gif_shape():
     "+units=m +no_defs"
 )
 
-# test_metadata: list of (variable,expected,tolerance) tuples
-test_metadata = [
-    ("projection", expected_proj1, None),
+# list of (variable,expected,tolerance) tuples
+test_attrs = [
+    ("projection", expected_proj, None),
+    ("institution", "MeteoSwiss", None),
+    ("accutime", 5.0, 0.1),
+    ("unit", "mm", None),
+    ("transform", None, None),
+    ("zerovalue", 0.0, 0.1),
+    ("threshold", 0.0009628129986471908, 1e-19),
+    ("zr_a", 316.0, 0.1),
+    ("zr_b", 1.5, 0.1),
     ("x1", 255000.0, 0.1),
     ("y1", -160000.0, 0.1),
     ("x2", 965000.0, 0.1),
@@ -37,37 +46,18 @@ def test_io_import_mch_gif_shape():
     ("ypixelsize", 1000.0, 0.1),
     ("cartesian_unit", "m", None),
     ("yorigin", "upper", None),
-    ("accutime", 5.0, 0.1),
-    ("unit", "mm", None),
-    ("transform", None, None),
-    ("zerovalue", 0.0, None),
-    ("threshold", 0.0009628129986471908, 1e-19),
-    ("institution", "MeteoSwiss", None),
-    ("product", "AQC", None),
-    ("zr_a", 316.0, None),
-    ("zr_b", 1.5, None),
 ]
 
 
-@pytest.mark.parametrize("variable, expected, tolerance", test_metadata)
-def test_io_import_mch_gif_metadata(variable, expected, tolerance):
+@pytest.mark.parametrize("variable, expected, tolerance", test_attrs)
+def test_io_import_mch_gif_dataset_attrs(variable, expected, tolerance):
     """Test the importer MCH GIF."""
-    root_path = pysteps.rcparams.data_sources["mch"]["root_path"]
-    filename = os.path.join(root_path, "20170131", "AQC170310945F_00005.801.gif")
-    _, _, metadata = pysteps.io.import_mch_gif(filename, "AQC", "mm", 5.0)
     smart_assert(metadata[variable], expected, tolerance)
 
 
-expected_proj2 = (
-    "+proj=somerc  +lon_0=7.43958333333333 +lat_0=46.9524055555556 "
-    "+k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel "
-    "+towgs84=674.374,15.056,405.346,0,0,0,0 "
-    "+units=m +no_defs"
-)
-
 # test_geodata: list of (variable,expected,tolerance) tuples
 test_geodata = [
-    ("projection", expected_proj2, None),
+    ("projection", expected_proj, None),
     ("x1", 255000.0, 0.1),
     ("y1", -160000.0, 0.1),
     ("x2", 965000.0, 0.1),
diff --git a/pysteps/tests/test_io_mrms_grib.py b/pysteps/tests/test_io_mrms_grib.py
index 65dd0caeb..cc7e20a84 100644
--- a/pysteps/tests/test_io_mrms_grib.py
+++ b/pysteps/tests/test_io_mrms_grib.py
@@ -15,18 +15,16 @@ def test_io_import_mrms_grib():
     """Test the importer for NSSL data."""
 
     root_path = pysteps.rcparams.data_sources["mrms"]["root_path"]
-
     filename = os.path.join(
         root_path, "2019/06/10/", "PrecipRate_00.00_20190610-000000.grib2"
     )
+    precip, _, metadata = pysteps.io.import_mrms_grib(filename, fillna=0, window_size=1)
 
-    precip_full, _, metadata = pysteps.io.import_mrms_grib(
-        filename, fillna=0, window_size=1
-    )
-    assert precip_full.shape == (3500, 7000)
-    assert precip_full.dtype == "single"
+    assert precip.shape == (3500, 7000)
+    assert precip.dtype == "single"
 
     expected_metadata = {
+        "institution": "NOAA National Severe Storms Laboratory",
         "xpixelsize": 0.01,
         "ypixelsize": 0.01,
         "unit": "mm/h",
@@ -50,64 +48,36 @@ def test_io_import_mrms_grib():
 
     x = np.arange(metadata["x1"], metadata["x2"], metadata["xpixelsize"])
     y = np.arange(metadata["y1"], metadata["y2"], metadata["ypixelsize"])
-    assert y.size == precip_full.shape[0]
-    assert x.size == precip_full.shape[1]
+
+    assert y.size == precip.shape[0]
+    assert x.size == precip.shape[1]
 
     # The full latitude range is (20.005, 54.995)
     # The full longitude range is (230.005, 299.995)
 
     # Test that if the bounding box is larger than the domain, all the points are returned.
-    precip_full2 = pysteps.io.import_mrms_grib(
+    precip2, _, _ = pysteps.io.import_mrms_grib(
         filename, fillna=0, extent=(220, 300, 20, 55), window_size=1
-    )[0]
-    assert precip_full2.shape == (3500, 7000)
+    )
+    assert precip2.shape == (3500, 7000)
 
-    assert_array_almost_equal(precip_full, precip_full2)
+    assert_array_almost_equal(precip, precip2)
 
-    del precip_full2
+    del precip2
 
     # Test that a portion of the domain is returned correctly
-    precip_clipped = pysteps.io.import_mrms_grib(
+    precip3, _, _ = pysteps.io.import_mrms_grib(
         filename, fillna=0, extent=(250, 260, 30, 35), window_size=1
-    )[0]
-
-    assert precip_clipped.shape == (500, 1000)
-    assert_array_almost_equal(precip_clipped, precip_full[2000:2500, 2000:3000])
-    del precip_clipped
-
-    precip_single = pysteps.io.import_mrms_grib(filename, dtype="double", fillna=0)[0]
-    assert precip_single.dtype == "double"
-    del precip_single
-
-    precip_single = pysteps.io.import_mrms_grib(filename, dtype="single", fillna=0)[0]
-    assert precip_single.dtype == "single"
-    del precip_single
-
-    precip_donwscaled = pysteps.io.import_mrms_grib(
-        filename, dtype="single", fillna=0, window_size=2
-    )[0]
-    assert precip_donwscaled.shape == (3500 / 2, 7000 / 2)
-
-    precip_donwscaled, _, metadata = 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,
     )
 
-    # 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)
+    assert precip3.shape == (500, 1000)
+    assert_array_almost_equal(precip3, precip[2000:2500, 2000:3000])
+    del precip3
+
+    precip4, _, _ = pysteps.io.import_mrms_grib(filename, dtype="double", fillna=0)
+    assert precip4.dtype == "double"
+    del precip4
 
-    for key in expected_metadata.keys():
-        assert metadata[key] == expected_metadata[key]
-    assert precip_donwscaled.shape == (3500 // 3, 7000 // 3)
+    precip5, _, _ = pysteps.io.import_mrms_grib(filename, dtype="single", fillna=0)
+    assert precip5.dtype == "single"
+    del precip5
diff --git a/pysteps/tests/test_io_opera_hdf5.py b/pysteps/tests/test_io_opera_hdf5.py
index 6411bcd33..6b2214d49 100644
--- a/pysteps/tests/test_io_opera_hdf5.py
+++ b/pysteps/tests/test_io_opera_hdf5.py
@@ -10,26 +10,28 @@
 pytest.importorskip("h5py")
 
 
+root_path = pysteps.rcparams.data_sources["opera"]["root_path"]
+filename = os.path.join(root_path, "20180824", "T_PAAH21_C_EUOC_20180824180000.hdf")
+precip, _, metadata = pysteps.io.import_opera_hdf5(filename)
+
+
 def test_io_import_opera_hdf5_shape():
     """Test the importer OPERA HDF5."""
-
-    root_path = pysteps.rcparams.data_sources["opera"]["root_path"]
-    filename = os.path.join(root_path, "20180824", "T_PAAH21_C_EUOC_20180824180000.hdf")
-    R, _, _ = pysteps.io.import_opera_hdf5(filename)
-    assert R.shape == (2200, 1900)
+    assert precip.shape == (2200, 1900)
 
 
 # test_metadata: list of (variable,expected, tolerance) tuples
 
-expected_proj1 = (
+expected_proj = (
     "+proj=laea +lat_0=55.0 +lon_0=10.0 "
     "+x_0=1950000.0 "
     "+y_0=-2100000.0 "
     "+units=m +ellps=WGS84"
 )
 
-test_metadata = [
-    ("projection", expected_proj1, None),
+# list of (variable,expected,tolerance) tuples
+test_attrs = [
+    ("projection", expected_proj, None),
     ("ll_lon", -10.434576838640398, 1e-10),
     ("ll_lat", 31.746215319325056, 1e-10),
     ("ur_lon", 57.81196475014995, 1e-10),
@@ -51,12 +53,7 @@ def test_io_import_opera_hdf5_shape():
 ]
 
 
-@pytest.mark.parametrize("variable,expected,tolerance", test_metadata)
-def test_io_import_opera_hdf5_metadata(variable, expected, tolerance):
+@pytest.mark.parametrize("variable, expected, tolerance", test_attrs)
+def test_io_import_mch_gif_dataset_attrs(variable, expected, tolerance):
     """Test the importer OPERA HDF5."""
-    root_path = pysteps.rcparams.data_sources["opera"]["root_path"]
-    filename = os.path.join(root_path, "20180824", "T_PAAH21_C_EUOC_20180824180000.hdf")
-    _, _, metadata = pysteps.io.import_opera_hdf5(filename)
-    print(metadata)
-    # assert metadata[variable] == expected
     smart_assert(metadata[variable], expected, tolerance)
diff --git a/pysteps/tests/test_io_readers.py b/pysteps/tests/test_io_readers.py
new file mode 100644
index 000000000..7decde461
--- /dev/null
+++ b/pysteps/tests/test_io_readers.py
@@ -0,0 +1,39 @@
+from datetime import datetime
+
+import numpy as np
+import pytest
+
+import pysteps
+
+
+def test_read_timeseries_mch():
+
+    pytest.importorskip("PIL")
+
+    date = datetime.strptime("201505151630", "%Y%m%d%H%M")
+    data_source = pysteps.rcparams.data_sources["mch"]
+    root_path = data_source["root_path"]
+    path_fmt = data_source["path_fmt"]
+    fn_pattern = data_source["fn_pattern"]
+    fn_ext = data_source["fn_ext"]
+    importer_name = data_source["importer"]
+    importer_kwargs = data_source["importer_kwargs"]
+    timestep = data_source["timestep"]
+
+    fns = pysteps.io.archive.find_by_date(
+        date,
+        root_path,
+        path_fmt,
+        fn_pattern,
+        fn_ext,
+        timestep=timestep,
+        num_prev_files=1,
+        num_next_files=1,
+    )
+
+    importer = pysteps.io.get_method(importer_name, "importer")
+    precip, _, metadata = pysteps.io.read_timeseries(fns, importer, **importer_kwargs)
+
+    assert isinstance(precip, np.ndarray)
+    assert isinstance(metadata, dict)
+    assert precip.shape[0] == 3
diff --git a/pysteps/tests/test_io_saf_crri.py b/pysteps/tests/test_io_saf_crri.py
index 66faa29d2..f9b9c249d 100644
--- a/pysteps/tests/test_io_saf_crri.py
+++ b/pysteps/tests/test_io_saf_crri.py
@@ -36,29 +36,30 @@ def test_io_import_saf_crri_geodata(variable, expected, tolerance):
         root_path, rel_path, "S_NWC_CRR_MSG4_Europe-VISIR_20180601T070000Z.nc"
     )
     geodata = pysteps.io.importers._import_saf_crri_geodata(filename)
-    print(variable, geodata[variable], expected, tolerance)
     smart_assert(geodata[variable], expected, tolerance)
 
 
-test_metadata_crri = [
+root_path = pysteps.rcparams.data_sources["saf"]["root_path"]
+rel_path = "20180601/CRR"
+filename = os.path.join(
+    root_path, rel_path, "S_NWC_CRR_MSG4_Europe-VISIR_20180601T070000Z.nc"
+)
+_, _, metadata = pysteps.io.import_saf_crri(filename)
+
+# list of (variable,expected,tolerance) tuples
+test_attrs = [
+    ("projection", expected_proj, None),
+    ("institution", "Agencia Estatal de Meteorología (AEMET)", None),
     ("transform", None, None),
     ("zerovalue", 0.0, 0.1),
     ("unit", "mm/h", None),
     ("accutime", None, None),
-    ("institution", "Agencia Estatal de Meteorología (AEMET)", None),
 ]
 
 
-@pytest.mark.parametrize("variable, expected, tolerance", test_metadata_crri)
-def test_io_import_saf_crri_metadata(variable, expected, tolerance):
+@pytest.mark.parametrize("variable, expected, tolerance", test_attrs)
+def test_io_import_saf_crri_attrs(variable, expected, tolerance):
     """Test the importer SAF CRRI."""
-    root_path = pysteps.rcparams.data_sources["saf"]["root_path"]
-    rel_path = "20180601/CRR"
-    filename = os.path.join(
-        root_path, rel_path, "S_NWC_CRR_MSG4_Europe-VISIR_20180601T070000Z.nc"
-    )
-    _, _, metadata = pysteps.io.import_saf_crri(filename)
-    print(variable, metadata[variable], expected, tolerance)
     smart_assert(metadata[variable], expected, tolerance)
 
 
@@ -85,20 +86,5 @@ def test_io_import_saf_crri_extent(extent, expected_extent, expected_shape, tole
     )
     precip, _, metadata = pysteps.io.import_saf_crri(filename, extent=extent)
     extent_out = (metadata["x1"], metadata["x2"], metadata["y1"], metadata["y2"])
-
-    print(extent, extent_out, expected_extent)
     smart_assert(extent_out, expected_extent, tolerance)
-    print(precip.shape, expected_shape)
     smart_assert(precip.shape, expected_shape, tolerance)
-
-
-if __name__ == "__main__":
-
-    for i, args in enumerate(test_geodata_crri):
-        test_io_import_saf_crri_geodata(args[0], args[1], args[2])
-
-    for i, args in enumerate(test_metadata_crri):
-        test_io_import_saf_crri_metadata(args[0], args[1], args[2])
-
-    for i, args in enumerate(test_extent_crri):
-        test_io_import_saf_crri_extent(args[0], args[1], args[2], args[3])
diff --git a/pysteps/tests/test_nowcasts_steps.py b/pysteps/tests/test_nowcasts_steps.py
index ec01998db..09833e2a8 100644
--- a/pysteps/tests/test_nowcasts_steps.py
+++ b/pysteps/tests/test_nowcasts_steps.py
@@ -9,6 +9,7 @@
 from pysteps import io, motion, nowcasts, verification
 from pysteps.tests.helpers import get_precipitation_fields
 
+
 steps_arg_names = (
     "n_ens_members",
     "n_cascade_levels",
diff --git a/pysteps/tests/test_plugins_support.py b/pysteps/tests/test_plugins_support.py
index 91f9f8de5..bbf3bcf5b 100644
--- a/pysteps/tests/test_plugins_support.py
+++ b/pysteps/tests/test_plugins_support.py
@@ -10,6 +10,7 @@
 import sys
 import tempfile
 
+
 __ = pytest.importorskip("cookiecutter")
 from cookiecutter.main import cookiecutter
 
diff --git a/pysteps/tests/test_utils_reprojection.py b/pysteps/tests/test_utils_reprojection.py
new file mode 100644
index 000000000..b80ff4d65
--- /dev/null
+++ b/pysteps/tests/test_utils_reprojection.py
@@ -0,0 +1,113 @@
+# -*- coding: utf-8 -*-
+
+import os
+import numpy as np
+import pytest
+import pysteps
+from pysteps.utils import reprojection as rpj
+
+pytest.importorskip("rasterio")
+
+root_path_radar = pysteps.rcparams.data_sources["rmi"]["root_path"]
+
+rel_path_radar = "20210704"  # Different date, but that does not matter for the tester
+
+filename_radar = os.path.join(
+    root_path_radar, rel_path_radar, "20210704180500.rad.best.comp.rate.qpe.hdf"
+)
+
+# Open the radar data
+radar_array, _, metadata_dst = pysteps.io.importers.import_odim_hdf5(filename_radar)
+
+# Initialise dummy NWP data
+nwp_array = np.zeros((24, 564, 564))
+
+for t in range(nwp_array.shape[0]):
+    nwp_array[t, 30 + t : 185 + t, 30 + 2 * t] = 0.1
+    nwp_array[t, 30 + t : 185 + t, 31 + 2 * t] = 0.1
+    nwp_array[t, 30 + t : 185 + t, 32 + 2 * t] = 1.0
+    nwp_array[t, 30 + t : 185 + t, 33 + 2 * t] = 5.0
+    nwp_array[t, 30 + t : 185 + t, 34 + 2 * t] = 5.0
+    nwp_array[t, 30 + t : 185 + t, 35 + 2 * t] = 4.5
+    nwp_array[t, 30 + t : 185 + t, 36 + 2 * t] = 4.5
+    nwp_array[t, 30 + t : 185 + t, 37 + 2 * t] = 4.0
+    nwp_array[t, 30 + t : 185 + t, 38 + 2 * t] = 2.0
+    nwp_array[t, 30 + t : 185 + t, 39 + 2 * t] = 1.0
+    nwp_array[t, 30 + t : 185 + t, 40 + 2 * t] = 0.5
+    nwp_array[t, 30 + t : 185 + t, 41 + 2 * t] = 0.1
+
+nwp_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.000000001"
+)
+
+metadata_src = dict(
+    projection=nwp_proj,
+    institution="Royal Meteorological Institute of Belgium",
+    transform=None,
+    zerovalue=0.0,
+    threshold=0,
+    unit="mm",
+    accutime=5,
+    xpixelsize=1300.0,
+    ypixelsize=1300.0,
+    yorigin="upper",
+    cartesian_unit="m",
+    x1=0.0,
+    x2=731900.0,
+    y1=-731900.0,
+    y2=0.0,
+)
+
+steps_arg_names = (
+    "radar_array",
+    "nwp_array",
+    "metadata_src",
+    "metadata_dst",
+)
+
+steps_arg_values = [
+    (radar_array, nwp_array, metadata_src, metadata_dst),
+]
+
+
+@pytest.mark.parametrize(steps_arg_names, steps_arg_values)
+def test_utils_reproject_grids(
+    radar_array,
+    nwp_array,
+    metadata_src,
+    metadata_dst,
+):
+
+    # Reproject
+    nwp_array_reproj, metadata_reproj = rpj.reproject_grids(
+        nwp_array, radar_array, metadata_src, metadata_dst
+    )
+
+    # The tests
+    assert (
+        nwp_array_reproj.shape[0] == nwp_array.shape[0]
+    ), "Time dimension has not the same length as source"
+    assert (
+        nwp_array_reproj.shape[1] == radar_array.shape[0]
+    ), "y dimension has not the same length as radar composite"
+    assert (
+        nwp_array_reproj.shape[2] == radar_array.shape[1]
+    ), "x dimension has not the same length as radar composite"
+
+    assert (
+        metadata_reproj["x1"] == metadata_dst["x1"]
+    ), "x-value lower left corner is not equal to radar composite"
+    assert (
+        metadata_reproj["x2"] == metadata_dst["x2"]
+    ), "x-value upper right corner is not equal to radar composite"
+    assert (
+        metadata_reproj["y1"] == metadata_dst["y1"]
+    ), "y-value lower left corner is not equal to radar composite"
+    assert (
+        metadata_reproj["y2"] == metadata_dst["y2"]
+    ), "y-value upper right corner is not equal to radar composite"
+
+    assert (
+        metadata_reproj["projection"] == metadata_dst["projection"]
+    ), "projection is different than destionation projection"
diff --git a/pysteps/tests/test_verification_detcatscores.py b/pysteps/tests/test_verification_detcatscores.py
index 2126df6bb..ba14b23e6 100644
--- a/pysteps/tests/test_verification_detcatscores.py
+++ b/pysteps/tests/test_verification_detcatscores.py
@@ -6,6 +6,7 @@
 
 from pysteps.verification import det_cat_fct
 
+
 # CREATE A LARGE DATASET TO MATCH
 # EXAMPLES IN
 # http://www.cawcr.gov.au/projects/verification/