diff --git a/pysteps/io/nowcast_importers.py b/pysteps/io/nowcast_importers.py
index 8f5e8ac85..2253c46e0 100755
--- a/pysteps/io/nowcast_importers.py
+++ b/pysteps/io/nowcast_importers.py
@@ -72,6 +72,7 @@
 
 from pysteps.decorators import postprocess_import
 from pysteps.exceptions import MissingOptionalDependency, DataModelError
+import xarray as xr
 
 try:
     import netCDF4
@@ -81,7 +82,6 @@
     NETCDF4_IMPORTED = False
 
 
-@postprocess_import(dtype="single")
 def import_netcdf_pysteps(filename, onerror="warn", **kwargs):
     """
     Read a nowcast or an ensemble of nowcasts from a NetCDF file conforming
@@ -119,122 +119,12 @@ def import_netcdf_pysteps(filename, onerror="warn", **kwargs):
     onerror = onerror.lower()
     if onerror not in ["warn", "raise"]:
         raise ValueError("'onerror' keyword must be 'warn' or 'raise'.")
-
     try:
-        ds = netCDF4.Dataset(filename, "r")
-
-        var_names = list(ds.variables.keys())
-
-        if "precip_intensity" in var_names:
-            precip = ds.variables["precip_intensity"]
-            unit = "mm/h"
-            accutime = None
-            transform = None
-        elif "precip_accum" in var_names:
-            precip = ds.variables["precip_accum"]
-            unit = "mm"
-            accutime = None
-            transform = None
-        elif "hourly_precip_accum" in var_names:
-            precip = ds.variables["hourly_precip_accum"]
-            unit = "mm"
-            accutime = 60.0
-            transform = None
-        elif "reflectivity" in var_names:
-            precip = ds.variables["reflectivity"]
-            unit = "dBZ"
-            accutime = None
-            transform = "dB"
-        else:
-            raise DataModelError(
-                "Non CF compilant file: "
-                "the netCDF file does not contain any "
-                "supported variable name.\n"
-                "Supported names: 'precip_intensity', 'hourly_precip_accum', "
-                "or 'reflectivity'\n"
-                "file: " + filename
-            )
-
-        precip = precip[...].squeeze().astype(float)
-
-        if isinstance(precip, np.ma.MaskedArray):
-            invalid_mask = np.ma.getmaskarray(precip)
-            precip = precip.data
-            precip[invalid_mask] = np.nan
-
-        metadata = {}
-
-        time_var = ds.variables["time"]
-        leadtimes = time_var[:] / 60.0  # minutes leadtime
-        metadata["leadtimes"] = leadtimes
-        timestamps = netCDF4.num2date(time_var[:], time_var.units)
-        metadata["timestamps"] = timestamps
-
-        if "polar_stereographic" in var_names:
-            vn = "polar_stereographic"
-
-            attr_dict = {}
-            for attr_name in ds.variables[vn].ncattrs():
-                attr_dict[attr_name] = ds[vn].getncattr(attr_name)
-
-            proj_str = _convert_grid_mapping_to_proj4(attr_dict)
-            metadata["projection"] = proj_str
-
-        # geodata
-        metadata["xpixelsize"] = abs(ds.variables["x"][1] - ds.variables["x"][0])
-        metadata["ypixelsize"] = abs(ds.variables["y"][1] - ds.variables["y"][0])
-
-        xmin = np.min(ds.variables["x"]) - 0.5 * metadata["xpixelsize"]
-        xmax = np.max(ds.variables["x"]) + 0.5 * metadata["xpixelsize"]
-        ymin = np.min(ds.variables["y"]) - 0.5 * metadata["ypixelsize"]
-        ymax = np.max(ds.variables["y"]) + 0.5 * metadata["ypixelsize"]
-
-        # TODO: this is only a quick solution
-        metadata["x1"] = xmin
-        metadata["y1"] = ymin
-        metadata["x2"] = xmax
-        metadata["y2"] = ymax
-
-        metadata["yorigin"] = "upper"  # TODO: check this
-
-        # TODO: Read the metadata to the dictionary.
-        if (accutime is None) and (leadtimes.size > 1):
-            accutime = leadtimes[1] - leadtimes[0]
-        metadata["accutime"] = accutime
-        metadata["unit"] = unit
-        metadata["transform"] = transform
-        metadata["zerovalue"] = np.nanmin(precip)
-        if metadata["zerovalue"] == np.nanmax(precip):
-            metadata["threshold"] = metadata["zerovalue"]
-        else:
-            metadata["threshold"] = np.nanmin(precip[precip > metadata["zerovalue"]])
-
-        ds.close()
-
-        return precip, metadata
+        dataset = xr.open_dataset(filename)
+        return dataset
     except Exception as er:
         if onerror == "warn":
             print("There was an error processing the file", er)
             return None, None
         else:
             raise er
-
-
-def _convert_grid_mapping_to_proj4(grid_mapping):
-    gm_keys = list(grid_mapping.keys())
-
-    # TODO: implement more projection types here
-    if grid_mapping["grid_mapping_name"] == "polar_stereographic":
-        proj_str = "+proj=stere"
-        proj_str += " +lon_0=%s" % grid_mapping["straight_vertical_longitude_from_pole"]
-        proj_str += " +lat_0=%s" % grid_mapping["latitude_of_projection_origin"]
-        if "standard_parallel" in gm_keys:
-            proj_str += " +lat_ts=%s" % grid_mapping["standard_parallel"]
-        if "scale_factor_at_projection_origin" in gm_keys:
-            proj_str += " +k_0=%s" % grid_mapping["scale_factor_at_projection_origin"]
-        proj_str += " +x_0=%s" % grid_mapping["false_easting"]
-        proj_str += " +y_0=%s" % grid_mapping["false_northing"]
-
-        return proj_str
-    else:
-        return None
diff --git a/pysteps/pystepsrc b/pysteps/pystepsrc
index 3df9ec288..47f4cabcf 100644
--- a/pysteps/pystepsrc
+++ b/pysteps/pystepsrc
@@ -49,7 +49,7 @@
         "fmi_geotiff": {
             "root_path": "./radar/fmi/geotiff",
             "path_fmt": "%Y%m%d",
-            "fn_pattern": "%Y%m%d%H%M_FINUTM.tif",
+            "fn_pattern": "%Y%m%d%H%M_FINUTM",
             "fn_ext": "tif",
             "importer": "geotiff",
             "timestep": 5,
diff --git a/pysteps/tests/helpers.py b/pysteps/tests/helpers.py
index 52fe940dc..be4805003 100644
--- a/pysteps/tests/helpers.py
+++ b/pysteps/tests/helpers.py
@@ -19,12 +19,14 @@
 _reference_dates = dict()
 _reference_dates["bom"] = datetime(2018, 6, 16, 10, 0)
 _reference_dates["fmi"] = datetime(2016, 9, 28, 16, 0)
+_reference_dates["fmi_geotiff"] = datetime(2016, 9, 28, 16, 0)
 _reference_dates["knmi"] = datetime(2010, 8, 26, 0, 0)
 _reference_dates["mch"] = datetime(2015, 5, 15, 16, 30)
 _reference_dates["dwd"] = datetime(2025, 6, 4, 17, 0)
 _reference_dates["opera"] = datetime(2018, 8, 24, 18, 0)
 _reference_dates["saf"] = datetime(2018, 6, 1, 7, 0)
 _reference_dates["mrms"] = datetime(2019, 6, 10, 0, 0)
+_reference_dates["rmi"] = datetime(2021, 7, 4, 18, 5)
 
 
 def assert_dataset_equivalent(dataset1: xr.Dataset, dataset2: xr.Dataset) -> None:
@@ -65,11 +67,14 @@ def get_precipitation_fields(
     Get a precipitation field from the archive to be used as reference.
 
     Source: bom
-    Reference time: 2018/06/16 10000 UTC
+    Reference time: 2018/06/16 1000 UTC
 
     Source: fmi
     Reference time: 2016/09/28 1600 UTC
 
+    Source: fmi_geotiff
+    Reference time: 2016/09/28 1600 UTC
+
     Source: knmi
     Reference time: 2010/08/26 0000 UTC
 
@@ -88,6 +93,9 @@ def get_precipitation_fields(
     Source: mrms
     Reference time: 2019/06/10 0000 UTC
 
+    Source: rmi
+    Reference time: 2021/07/04 1805 UTC
+
     Parameters
     ----------
 
@@ -110,7 +118,7 @@ def get_precipitation_fields(
         If it is a float, represents the length of the space window that is
         used to upscale the fields.
 
-    source: {"bom", "fmi" , "knmi", "mch", "opera", "saf", "mrms"}, optional
+    source: {"bom", "fmi" , "fmi_geotiff", "knmi", "mch", "dwd", "opera", "saf", "mrms", "rmi"}, optional
         Name of the data source to be used.
 
     log_transform: bool
@@ -138,6 +146,10 @@ def get_precipitation_fields(
     if source == "fmi":
         pytest.importorskip("pyproj")
 
+    if source == "fmi_geotiff":
+        pytest.importorskip("pyproj")
+        pytest.importorskip("osgeo")
+
     if source == "knmi":
         pytest.importorskip("h5py")
 
@@ -156,6 +168,10 @@ def get_precipitation_fields(
     if source == "mrms":
         pytest.importorskip("pygrib")
 
+    if source == "rmi":
+        pytest.importorskip("rasterio")
+        pytest.importorskip("pyproj")
+
     try:
         date = _reference_dates[source]
     except KeyError:
diff --git a/pysteps/tests/test_ensscores.py b/pysteps/tests/test_ensscores.py
index e1e7b89e8..ea00d1f34 100644
--- a/pysteps/tests/test_ensscores.py
+++ b/pysteps/tests/test_ensscores.py
@@ -7,9 +7,14 @@
 from pysteps.tests.helpers import get_precipitation_fields
 from pysteps.verification import ensscores
 
-precip = get_precipitation_fields(num_next_files=10, return_raw=True)
+precip_dataset = get_precipitation_fields(num_next_files=10, return_raw=True)
 np.random.seed(42)
 
+# XR: extract values from array because score functions are not xarray compatible
+precip_var = precip_dataset.attrs["precip_var"]
+precip = precip_dataset[precip_var].values
+
+
 # rankhist
 test_data = [
     (precip[:10], precip[-1], None, True, 11),
diff --git a/pysteps/tests/test_io_dwd_hdf5.py b/pysteps/tests/test_io_dwd_hdf5.py
index e86a22f01..59de1692d 100644
--- a/pysteps/tests/test_io_dwd_hdf5.py
+++ b/pysteps/tests/test_io_dwd_hdf5.py
@@ -2,27 +2,26 @@
 
 import pytest
 
-import pysteps
 from pysteps.tests.helpers import smart_assert, get_precipitation_fields
 
-pytest.importorskip("h5py")
-
 # Test for RADOLAN RY product
-
-precip_ry, metadata_ry = get_precipitation_fields(
+precip_dataset = get_precipitation_fields(
     num_prev_files=0,
     num_next_files=0,
-    return_raw=False,
+    return_raw=True,
     metadata=True,
     source="dwd",
     log_transform=False,
     importer_kwargs=dict(qty="RATE"),
 )
 
+precip_var = precip_dataset.attrs["precip_var"]
+precip_dataarray = precip_dataset[precip_var]
+
 
 def test_io_import_dwd_hdf5_ry_shape():
     """Test the importer DWD HDF5."""
-    assert precip_ry.shape == (1200, 1100)
+    assert precip_dataarray.shape == (1, 1200, 1100)
 
 
 # Test_metadata
@@ -36,29 +35,34 @@ def test_io_import_dwd_hdf5_ry_shape():
 
 # List of (variable,expected,tolerance) tuples
 test_ry_attrs = [
-    ("projection", expected_proj, None),
-    ("ll_lon", 3.566994635, 1e-10),
-    ("ll_lat", 45.69642538, 1e-10),
-    ("ur_lon", 18.73161645, 1e-10),
-    ("ur_lat", 55.84543856, 1e-10),
-    ("x1", -500.0, 1e-6),
-    ("y1", -1199500.0, 1e-6),
-    ("x2", 1099500.0, 1e-6),
-    ("y2", 500.0, 1e-6),
-    ("xpixelsize", 1000.0, 1e-10),
-    ("xpixelsize", 1000.0, 1e-10),
-    ("cartesian_unit", "m", None),
-    ("yorigin", "upper", None),
-    ("institution", "ORG:78,CTY:616,CMT:Deutscher Wetterdienst radolan@dwd.de", None),
-    ("accutime", 5.0, 1e-10),
-    ("unit", "mm/h", None),
-    ("transform", None, None),
-    ("zerovalue", 0.0, 1e-6),
-    ("threshold", 0.12, 1e-6),
+    (precip_dataset.attrs["projection"], expected_proj, None),
+    (float(precip_dataset.lon.isel(x=0, y=0).values), 3.57220017, 1e-8),
+    (float(precip_dataset.lat.isel(x=0, y=0).values), 45.70099971, 1e-8),
+    (float(precip_dataset.lon.isel(x=-1, y=-1).values), 18.72270377, 1e-8),
+    (float(precip_dataset.lat.isel(x=-1, y=-1).values), 55.84175857, 1e-8),
+    (precip_dataset.x.isel(x=0).values, 0.0, 1e-3),
+    (precip_dataset.y.isel(y=0).values, -1199000.0, 1e-6),
+    (precip_dataset.x.isel(x=-1).values, 1099000.0, 1e-6),
+    (precip_dataset.y.isel(y=-1).values, 0.0, 1e-3),
+    (precip_dataset.x.attrs["stepsize"], 1000.0, 1e-10),
+    (precip_dataset.y.attrs["stepsize"], 1000.0, 1e-10),
+    (precip_dataset.x.attrs["units"], "m", None),
+    (precip_dataset.y.attrs["units"], "m", None),
+    (
+        precip_dataset.attrs["institution"],
+        "ORG:78,CTY:616,CMT:Deutscher Wetterdienst radolan@dwd.de",
+        None,
+    ),
+    (precip_dataarray.attrs["accutime"], 5.0, 1e-10),
+    (precip_dataset.time.attrs["stepsize"], 5.0, 1e-10),
+    (precip_dataarray.attrs["units"], "mm/h", None),
+    (precip_dataarray.attrs["transform"], None, None),
+    (precip_dataarray.attrs["zerovalue"], 0.0, 1e-6),
+    (precip_dataarray.attrs["threshold"], 0.12, 1e-6),
 ]
 
 
 @pytest.mark.parametrize("variable, expected, tolerance", test_ry_attrs)
 def test_io_import_dwd_hdf5_ry_metadata(variable, expected, tolerance):
     """Test the importer OPERA HDF5."""
-    smart_assert(metadata_ry[variable], expected, tolerance)
+    smart_assert(variable, expected, tolerance)
diff --git a/pysteps/tests/test_io_fmi_geotiff.py b/pysteps/tests/test_io_fmi_geotiff.py
index fbcc3153c..2a07f03b0 100644
--- a/pysteps/tests/test_io_fmi_geotiff.py
+++ b/pysteps/tests/test_io_fmi_geotiff.py
@@ -1,25 +1,25 @@
-import os
+# -*- coding: utf-8 -*-
 
 import pytest
 
-import pysteps
-from pysteps.tests.helpers import smart_assert
+from pysteps.tests.helpers import smart_assert, get_precipitation_fields
 
-pytest.importorskip("pyproj")
-pytest.importorskip("osgeo")
-
-root_path = pysteps.rcparams.data_sources["fmi_geotiff"]["root_path"]
-filename = os.path.join(
-    root_path,
-    "20160928",
-    "201609281600_FINUTM.tif",
+precip_dataset = get_precipitation_fields(
+    num_prev_files=0,
+    num_next_files=0,
+    return_raw=True,
+    metadata=True,
+    source="fmi_geotiff",
+    log_transform=False,
 )
-precip, _, metadata = pysteps.io.import_fmi_geotiff(filename)
+
+precip_var = precip_dataset.attrs["precip_var"]
+precip_dataarray = precip_dataset[precip_var]
 
 
 def test_io_import_fmi_geotiff_shape():
     """Test the shape of the read file."""
-    assert precip.shape == (7316, 4963)
+    assert precip_dataarray.shape == (1, 7316, 4963)
 
 
 expected_proj = (
@@ -28,19 +28,20 @@ def test_io_import_fmi_geotiff_shape():
 
 # test_geodata: list of (variable,expected,tolerance) tuples
 test_geodata = [
-    ("projection", expected_proj, None),
-    ("x1", -196593.0043142295908183, 1e-10),
-    ("x2", 1044176.9413554778, 1e-10),
-    ("y1", 6255329.6988206729292870, 1e-10),
-    ("y2", 8084432.005259146, 1e-10),
-    ("xpixelsize", 250.0040188736061566, 1e-6),
-    ("ypixelsize", 250.0139839309011904, 1e-6),
-    ("cartesian_unit", "m", None),
-    ("yorigin", "upper", None),
+    (precip_dataset.attrs["projection"], expected_proj, None),
+    (precip_dataset.attrs["institution"], "Finnish Meteorological Institute", None),
+    (precip_dataset.x.isel(x=0).values, -196468.00230479, 1e-10),
+    (precip_dataset.y.isel(y=0).values, 6255454.70581264, 1e-10),
+    (precip_dataset.x.isel(x=-1).values, 1044051.93934604, 1e-10),
+    (precip_dataset.y.isel(y=-1).values, 8084306.99826718, 1e-10),
+    (precip_dataset.x.attrs["stepsize"], 250.0040188736061566, 1e-10),
+    (precip_dataset.y.attrs["stepsize"], 250.0139839309011904, 1e-10),
+    (precip_dataset.x.attrs["units"], "m", None),
+    (precip_dataset.y.attrs["units"], "m", None),
 ]
 
 
 @pytest.mark.parametrize("variable, expected, tolerance", test_geodata)
 def test_io_import_fmi_pgm_geodata(variable, expected, tolerance):
     """Test the GeoTIFF and metadata reading."""
-    smart_assert(metadata[variable], expected, tolerance)
+    smart_assert(variable, expected, tolerance)
diff --git a/pysteps/tests/test_io_fmi_pgm.py b/pysteps/tests/test_io_fmi_pgm.py
index f91fbc8a9..a704e0e50 100644
--- a/pysteps/tests/test_io_fmi_pgm.py
+++ b/pysteps/tests/test_io_fmi_pgm.py
@@ -1,25 +1,24 @@
-import os
-
+# -*- coding: utf-8 -*-
 import pytest
 
-import pysteps
-from pysteps.tests.helpers import smart_assert
-
-pytest.importorskip("pyproj")
-
+from pysteps.tests.helpers import smart_assert, get_precipitation_fields
 
-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_dataset = get_precipitation_fields(
+    num_prev_files=0,
+    num_next_files=0,
+    return_raw=True,
+    metadata=True,
+    source="fmi",
+    log_transform=False,
 )
-precip, _, metadata = pysteps.io.import_fmi_pgm(filename, gzipped=True)
+
+precip_var = precip_dataset.attrs["precip_var"]
+precip_dataarray = precip_dataset[precip_var]
 
 
 def test_io_import_fmi_pgm_shape():
     """Test the importer FMI PGM."""
-    assert precip.shape == (1226, 760)
+    assert precip_dataarray.shape == (1, 1226, 760)
 
 
 expected_proj = (
@@ -28,66 +27,39 @@ def test_io_import_fmi_pgm_shape():
     "+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),
+    (precip_dataset.attrs["projection"], expected_proj, None),
+    (precip_dataset.attrs["institution"], "Finnish Meteorological Institute", None),
+    (precip_dataarray.attrs["accutime"], 5.0, 1e-10),
+    (precip_dataarray.attrs["units"], "dBZ", None),
+    (precip_dataarray.attrs["transform"], "dB", None),
+    (precip_dataarray.attrs["zerovalue"], -32.0, 1e-6),
+    (precip_dataarray.attrs["threshold"], -31.5, 1e-6),
+    (precip_dataarray.attrs["zr_a"], 223.0, 1e-6),
+    (precip_dataarray.attrs["zr_b"], 1.53, 1e-6),
 ]
 
 
 @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)
+    smart_assert(variable, expected, tolerance)
 
 
 # test_geodata: list of (variable,expected,tolerance) tuples
 test_geodata = [
-    ("projection", expected_proj, None),
-    ("x1", 0.0049823258887045085, 1e-20),
-    ("x2", 759752.2852757066, 1e-10),
-    ("y1", 0.009731985162943602, 1e-18),
-    ("y2", 1225544.6588913496, 1e-10),
-    ("xpixelsize", 999.674053, 1e-6),
-    ("ypixelsize", 999.62859, 1e-5),
-    ("cartesian_unit", "m", None),
-    ("yorigin", "upper", None),
+    (precip_dataset.x.isel(x=0).values, 499.84200883, 1e-10),
+    (precip_dataset.y.isel(y=0).values, 499.8240261, 1e-10),
+    (precip_dataset.x.isel(x=-1).values, 759252.4482492, 1e-10),
+    (precip_dataset.y.isel(y=-1).values, 1225044.84459724, 1e-10),
+    (precip_dataset.x.attrs["stepsize"], 999.674053, 1e-8),
+    (precip_dataset.y.attrs["stepsize"], 999.62859, 1e-8),
+    (precip_dataset.x.attrs["units"], "m", None),
+    (precip_dataset.y.attrs["units"], "m", None),
 ]
 
 
 @pytest.mark.parametrize("variable, expected, tolerance", test_geodata)
 def test_io_import_fmi_pgm_geodata(variable, expected, tolerance):
-    """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",
-    )
-    metadata = pysteps.io.importers._import_fmi_pgm_metadata(filename, gzipped=True)
-    geodata = pysteps.io.importers._import_fmi_pgm_geodata(metadata)
-
-    smart_assert(geodata[variable], expected, tolerance)
+    smart_assert(variable, expected, tolerance)
diff --git a/pysteps/tests/test_io_knmi_hdf5.py b/pysteps/tests/test_io_knmi_hdf5.py
index 3e30cb575..25b96acdc 100644
--- a/pysteps/tests/test_io_knmi_hdf5.py
+++ b/pysteps/tests/test_io_knmi_hdf5.py
@@ -1,23 +1,26 @@
 # -*- coding: utf-8 -*-
 
-import os
-
 import pytest
 
-import pysteps
-from pysteps.tests.helpers import smart_assert
-
-pytest.importorskip("h5py")
+from pysteps.tests.helpers import smart_assert, get_precipitation_fields
 
+precip_dataset = get_precipitation_fields(
+    num_prev_files=0,
+    num_next_files=0,
+    return_raw=True,
+    metadata=True,
+    source="knmi",
+    log_transform=False,
+    importer_kwargs=dict(qty="ACRR"),
+)
 
-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)
+precip_var = precip_dataset.attrs["precip_var"]
+precip_dataarray = precip_dataset[precip_var]
 
 
 def test_io_import_knmi_hdf5_shape():
     """Test the importer KNMI HDF5."""
-    assert precip.shape == (765, 700)
+    assert precip_dataarray.shape == (1, 765, 700)
 
 
 # test_metadata: list of (variable,expected, tolerance) tuples
@@ -28,27 +31,32 @@ def test_io_import_knmi_hdf5_shape():
 
 # list of (variable,expected,tolerance) tuples
 test_attrs = [
-    ("projection", expected_proj, None),
-    ("x1", 0.0, 1e-10),
-    ("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),
-    ("institution", "KNMI - Royal Netherlands Meteorological Institute", None),
-    ("transform", None, None),
-    ("zerovalue", 0.0, 1e-10),
-    ("threshold", 0.01, 1e-10),
-    ("zr_a", 200.0, None),
-    ("zr_b", 1.6, None),
+    (precip_dataset.attrs["projection"], expected_proj, None),
+    (precip_dataset.x.isel(x=0).values, 499.98876045, 1e-10),
+    (precip_dataset.y.isel(y=0).values, -4414538.12181262, 1e-10),
+    (precip_dataset.x.isel(x=-1).values, 699484.27587271, 1e-10),
+    (precip_dataset.y.isel(y=-1).values, -3650450.41764577, 1e-10),
+    (precip_dataset.x.attrs["stepsize"], 1000.0, 1e-10),
+    (precip_dataset.y.attrs["stepsize"], 1000.0, 1e-10),
+    (precip_dataarray.attrs["accutime"], 5.0, 1e-10),
+    (precip_dataset.time.attrs["stepsize"], 5.0, 1e-10),
+    (precip_dataarray.attrs["units"], "mm/h", None),
+    (precip_dataset.x.attrs["units"], "m", None),
+    (precip_dataset.y.attrs["units"], "m", None),
+    (
+        precip_dataset.attrs["institution"],
+        "KNMI - Royal Netherlands Meteorological Institute",
+        None,
+    ),
+    (precip_dataarray.attrs["transform"], None, None),
+    (precip_dataarray.attrs["zerovalue"], 0.0, 1e-6),
+    (precip_dataarray.attrs["threshold"], 0.01, 1e-6),
+    (precip_dataarray.attrs["zr_a"], 200.0, None),
+    (precip_dataarray.attrs["zr_b"], 1.6, None),
 ]
 
 
 @pytest.mark.parametrize("variable,expected,tolerance", test_attrs)
 def test_io_import_knmi_hdf5_metadata(variable, expected, tolerance):
     """Test the importer KNMI HDF5."""
-    smart_assert(metadata[variable], expected, tolerance)
+    smart_assert(variable, expected, tolerance)
diff --git a/pysteps/tests/test_io_mch_gif.py b/pysteps/tests/test_io_mch_gif.py
index 903033988..13d71f9a5 100644
--- a/pysteps/tests/test_io_mch_gif.py
+++ b/pysteps/tests/test_io_mch_gif.py
@@ -1,22 +1,26 @@
 # -*- coding: utf-8 -*-
 
-import os
-
 import pytest
 
-import pysteps
-from pysteps.tests.helpers import smart_assert
+from pysteps.tests.helpers import smart_assert, get_precipitation_fields
 
-pytest.importorskip("PIL")
+precip_dataset = get_precipitation_fields(
+    num_prev_files=0,
+    num_next_files=0,
+    return_raw=True,
+    metadata=True,
+    source="mch",
+    log_transform=False,
+    importer_kwargs=dict(qty="AQC"),
+)
 
-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)
+precip_var = precip_dataset.attrs["precip_var"]
+precip_dataarray = precip_dataset[precip_var]
 
 
 def test_io_import_mch_gif_shape():
     """Test the importer MCH GIF."""
-    assert precip.shape == (640, 710)
+    assert precip_dataarray.shape == (1, 640, 710)
 
 
 expected_proj = (
@@ -29,48 +33,28 @@ def test_io_import_mch_gif_shape():
 
 # 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),
-    ("y2", 480000.0, 0.1),
-    ("xpixelsize", 1000.0, 0.1),
-    ("ypixelsize", 1000.0, 0.1),
-    ("cartesian_unit", "m", None),
-    ("yorigin", "upper", None),
+    (precip_dataset.attrs["projection"], expected_proj, None),
+    (precip_dataset.attrs["institution"], "MeteoSwiss", None),
+    (precip_dataarray.attrs["accutime"], 5.0, 1e-10),
+    (precip_dataset.time.attrs["stepsize"], 5.0, 1e-10),
+    (precip_dataarray.attrs["units"], "mm", None),
+    (precip_dataarray.attrs["transform"], None, None),
+    (precip_dataarray.attrs["zerovalue"], 0.0, 1e-6),
+    (precip_dataarray.attrs["threshold"], 0.0008258007600496956, 1e-19),
+    (precip_dataarray.attrs["zr_a"], 316.0, None),
+    (precip_dataarray.attrs["zr_b"], 1.5, None),
+    (precip_dataset.x.isel(x=0).values, 255500.0, 1e-10),
+    (precip_dataset.y.isel(y=0).values, -159500.0, 1e-10),
+    (precip_dataset.x.isel(x=-1).values, 964500.0, 1e-10),
+    (precip_dataset.y.isel(y=-1).values, 479500.0, 1e-10),
+    (precip_dataset.x.attrs["stepsize"], 1000.0, 0.1),
+    (precip_dataset.y.attrs["stepsize"], 1000.0, 0.1),
+    (precip_dataset.x.attrs["units"], "m", None),
+    (precip_dataset.y.attrs["units"], "m", None),
 ]
 
 
 @pytest.mark.parametrize("variable, expected, tolerance", test_attrs)
 def test_io_import_mch_gif_dataset_attrs(variable, expected, tolerance):
     """Test the importer MCH GIF."""
-    smart_assert(metadata[variable], expected, tolerance)
-
-
-# test_geodata: list of (variable,expected,tolerance) tuples
-test_geodata = [
-    ("projection", expected_proj, None),
-    ("x1", 255000.0, 0.1),
-    ("y1", -160000.0, 0.1),
-    ("x2", 965000.0, 0.1),
-    ("y2", 480000.0, 0.1),
-    ("xpixelsize", 1000.0, 0.1),
-    ("ypixelsize", 1000.0, 0.1),
-    ("cartesian_unit", "m", None),
-    ("yorigin", "upper", None),
-]
-
-
-@pytest.mark.parametrize("variable, expected, tolerance", test_geodata)
-def test_io_import_mch_geodata(variable, expected, tolerance):
-    """Test the importer MCH geodata."""
-    geodata = pysteps.io.importers._import_mch_geodata()
-    smart_assert(geodata[variable], expected, tolerance)
+    smart_assert(variable, expected, tolerance)
diff --git a/pysteps/tests/test_io_nowcast_importers.py b/pysteps/tests/test_io_nowcast_importers.py
index d99414f14..e35a5c26b 100644
--- a/pysteps/tests/test_io_nowcast_importers.py
+++ b/pysteps/tests/test_io_nowcast_importers.py
@@ -1,10 +1,12 @@
 import numpy as np
 import pytest
 
+from tempfile import NamedTemporaryFile
 from pysteps import io
 from pysteps.tests.helpers import get_precipitation_fields
+import xarray as xr
 
-precip, metadata = get_precipitation_fields(
+precip_dataset = get_precipitation_fields(
     num_prev_files=1,
     num_next_files=0,
     return_raw=False,
@@ -12,36 +14,57 @@
     upscale=2000,
 )
 
+# XR: The following lines are currently needed as we don't have a dedicated xarray writer, mayve we should fix the xarray reader to disallow non types in attrs and to store datetime units as encoding instead of attributes
+precip_dataset.time.encoding = {"units": precip_dataset.time.attrs["units"]}
+del precip_dataset.time.attrs["units"]
+
+precip_var = precip_dataset.attrs["precip_var"]
+precip_dataarray = precip_dataset[precip_var]
+
+zeros = np.zeros(precip_dataarray.shape, dtype=np.float32)
+
+zero_dataset = xr.Dataset(
+    data_vars={
+        "precip_intensity": (
+            ("time", "y", "x"),
+            zeros,
+            {
+                "long_name": "Precipitation intensity",
+                "units": "mm hr-1",  # keep attrs simple types, no None
+                "_FillValue": np.float32(-9999),  # valid NetCDF fill value
+                # omit standard_name unless you have a CF-valid value
+            },
+        )
+    },
+    coords={
+        "time": ("time", precip_dataarray["time"].values),
+        "y": ("y", precip_dataarray["y"].values),
+        "x": ("x", precip_dataarray["x"].values),
+    },
+    attrs={"precip_var": "precip_intensity"},  # simple, serializable globals
+)
+
 
 @pytest.mark.parametrize(
-    "precip, metadata",
-    [(precip, metadata), (np.zeros_like(precip), metadata)],
+    "precip_dataset",
+    [(precip_dataset), (zero_dataset)],
 )
-def test_import_netcdf(precip, metadata, tmp_path):
-
-    pytest.importorskip("pyproj")
-
-    field_shape = (precip.shape[1], precip.shape[2])
-    startdate = metadata["timestamps"][-1]
-    timestep = metadata["accutime"]
-    exporter = io.exporters.initialize_forecast_exporter_netcdf(
-        outpath=tmp_path.as_posix(),
-        outfnprefix="test",
-        startdate=startdate,
-        timestep=timestep,
-        n_timesteps=precip.shape[0],
-        shape=field_shape,
-        metadata=metadata,
-    )
-    io.exporters.export_forecast_dataset(precip, exporter)
-    io.exporters.close_forecast_files(exporter)
-
-    tmp_file = tmp_path / "test.nc"
-    precip_netcdf, metadata_netcdf = io.import_netcdf_pysteps(tmp_file, dtype="float64")
-
-    assert isinstance(precip_netcdf, np.ndarray)
-    assert isinstance(metadata_netcdf, dict)
-    assert precip_netcdf.ndim == precip.ndim, "Wrong number of dimensions"
-    assert precip_netcdf.shape[0] == precip.shape[0], "Wrong number of lead times"
-    assert precip_netcdf.shape[1:] == field_shape, "Wrong field shape"
-    assert np.allclose(precip_netcdf, precip)
+def test_import_netcdf(precip_dataset):
+    # XR: this test might not make that much sense in the future
+    with NamedTemporaryFile() as tempfile:
+        precip_var = precip_dataset.attrs["precip_var"]
+        precip_dataarray = precip_dataset[precip_var]
+        field_shape = (precip_dataarray.shape[1], precip_dataarray.shape[2])
+
+        precip_dataset.to_netcdf(tempfile.name)
+        precip_netcdf = io.import_netcdf_pysteps(tempfile.name, dtype="float64")
+
+        assert isinstance(precip_netcdf, xr.Dataset)
+        assert (
+            precip_netcdf[precip_var].ndim == precip_dataarray.ndim
+        ), "Wrong number of dimensions"
+        assert (
+            precip_netcdf[precip_var].shape[0] == precip_dataarray.shape[0]
+        ), "Wrong number of lead times"
+        assert precip_netcdf[precip_var].shape[1:] == field_shape, "Wrong field shape"
+        assert np.allclose(precip_netcdf[precip_var].values, precip_dataarray.values)
diff --git a/pysteps/tests/test_io_opera_hdf5.py b/pysteps/tests/test_io_opera_hdf5.py
index 52f86f7b8..dc62e637b 100644
--- a/pysteps/tests/test_io_opera_hdf5.py
+++ b/pysteps/tests/test_io_opera_hdf5.py
@@ -14,49 +14,55 @@
 # CIRRUS max. reflectivity composites
 # NIMBUS rain rate composites
 
+# XR: since the pysteps.datasets module does not support all the OPERA data sources below, we dont use get_precipitation_fields
 root_path = pysteps.rcparams.data_sources["opera"]["root_path"]
 
 filename = os.path.join(root_path, "20180824", "T_PAAH21_C_EUOC_20180824180000.hdf")
-precip_odyssey, _, metadata_odyssey = pysteps.io.import_opera_hdf5(filename, qty="RATE")
+precip_odyssey = pysteps.io.import_opera_hdf5(filename, qty="RATE")
+precip_var = precip_odyssey.attrs["precip_var"]
+precip_odyssey_dataarray = precip_odyssey[precip_var]
+
 
 filename = os.path.join(
     root_path, "20241126", "CIRRUS", "T_PABV21_C_EUOC_20241126010000.hdf"
 )
-precip_cirrus, _, metadata_cirrus = pysteps.io.import_opera_hdf5(filename, qty="DBZH")
+precip_cirrus = pysteps.io.import_opera_hdf5(filename, qty="DBZH")
+precip_var = precip_cirrus.attrs["precip_var"]
+precip_cirrus_dataarray = precip_cirrus[precip_var]
 
 filename = os.path.join(
     root_path, "20241126", "NIMBUS", "T_PAAH22_C_EUOC_20241126010000.hdf"
 )
-precip_nimbus_rain_rate, _, metadata_nimbus_rain_rate = pysteps.io.import_opera_hdf5(
-    filename, qty="RATE"
-)
+precip_nimbus_rain_rate = pysteps.io.import_opera_hdf5(filename, qty="RATE")
+precip_var = precip_nimbus_rain_rate.attrs["precip_var"]
+precip_nimbus_rain_rate_dataarray = precip_nimbus_rain_rate[precip_var]
 
 filename = os.path.join(
     root_path, "20241126", "NIMBUS", "T_PASH22_C_EUOC_20241126010000.hdf"
 )
-precip_nimbus_rain_accum, _, metadata_nimbus_rain_accum = pysteps.io.import_opera_hdf5(
-    filename, qty="ACRR"
-)
+precip_nimbus_rain_accum = pysteps.io.import_opera_hdf5(filename, qty="ACRR")
+precip_var = precip_nimbus_rain_accum.attrs["precip_var"]
+precip_nimbus_rain_accum_dataarray = precip_nimbus_rain_accum[precip_var]
 
 
 def test_io_import_opera_hdf5_odyssey_shape():
     """Test the importer OPERA HDF5."""
-    assert precip_odyssey.shape == (2200, 1900)
+    assert precip_odyssey_dataarray.shape == (2200, 1900)
 
 
 def test_io_import_opera_hdf5_cirrus_shape():
     """Test the importer OPERA HDF5."""
-    assert precip_cirrus.shape == (4400, 3800)
+    assert precip_cirrus_dataarray.shape == (4400, 3800)
 
 
 def test_io_import_opera_hdf5_nimbus_rain_rate_shape():
     """Test the importer OPERA HDF5."""
-    assert precip_nimbus_rain_rate.shape == (2200, 1900)
+    assert precip_nimbus_rain_rate_dataarray.shape == (2200, 1900)
 
 
 def test_io_import_opera_hdf5_nimbus_rain_accum_shape():
     """Test the importer OPERA HDF5."""
-    assert precip_nimbus_rain_accum.shape == (2200, 1900)
+    assert precip_nimbus_rain_accum_dataarray.shape == (2200, 1900)
 
 
 # test_metadata: list of (variable,expected, tolerance) tuples
@@ -69,85 +75,87 @@ def test_io_import_opera_hdf5_nimbus_rain_accum_shape():
 
 # list of (variable,expected,tolerance) tuples
 test_odyssey_attrs = [
-    ("projection", expected_proj, None),
-    ("ll_lon", -10.434576838640398, 1e-10),
-    ("ll_lat", 31.746215319325056, 1e-10),
-    ("ur_lon", 57.81196475014995, 1e-10),
-    ("ur_lat", 67.62103710275053, 1e-10),
-    ("x1", -0.0004161088727414608, 1e-6),
-    ("y1", -4400000.001057557, 1e-10),
-    ("x2", 3800000.0004256153, 1e-10),
-    ("y2", -0.0004262728616595268, 1e-6),
-    ("xpixelsize", 2000.0, 1e-10),
-    ("xpixelsize", 2000.0, 1e-10),
-    ("cartesian_unit", "m", None),
-    ("accutime", 15.0, 1e-10),
-    ("yorigin", "upper", None),
-    ("unit", "mm/h", None),
-    ("institution", "Odyssey datacentre", None),
-    ("transform", None, None),
-    ("zerovalue", 0.0, 1e-10),
-    ("threshold", 0.01, 1e-10),
+    (precip_odyssey.attrs["projection"], expected_proj, None),
+    (float(precip_odyssey.lon.isel(x=0, y=0).values), -10.4268122372, 1e-10),
+    (float(precip_odyssey.lat.isel(x=0, y=0).values), 31.7575305091, 1e-10),
+    (float(precip_odyssey.lon.isel(x=-1, y=-1).values), 57.7778944303, 1e-10),
+    (float(precip_odyssey.lat.isel(x=-1, y=-1).values), 67.6204665961, 1e-10),
+    (precip_odyssey.x.isel(x=0).values, 999.999583891, 1e-6),
+    (precip_odyssey.y.isel(y=0).values, -4399000.00106, 1e-10),
+    (precip_odyssey.x.isel(x=-1).values, 3799000.00043, 1e-10),
+    (precip_odyssey.y.isel(y=-1).values, -1000.00042627, 1e-6),
+    (precip_odyssey.x.attrs["stepsize"], 2000.0, 1e-10),
+    (precip_odyssey.y.attrs["stepsize"], 2000.0, 1e-10),
+    (precip_odyssey.x.attrs["units"], "m", None),
+    (precip_odyssey.y.attrs["units"], "m", None),
+    (precip_odyssey_dataarray.attrs["accutime"], 15.0, 1e-10),
+    # ("yorigin", "upper", None),
+    (precip_odyssey_dataarray.attrs["units"], "mm/h", None),
+    (precip_odyssey.attrs["institution"], "Odyssey datacentre", None),
+    (precip_odyssey_dataarray.attrs["transform"], None, None),
+    (precip_odyssey_dataarray.attrs["zerovalue"], 0.0, 1e-6),
+    (precip_odyssey_dataarray.attrs["threshold"], 0.01, 1e-6),
 ]
 
 
 @pytest.mark.parametrize("variable, expected, tolerance", test_odyssey_attrs)
 def test_io_import_opera_hdf5_odyssey_dataset_attrs(variable, expected, tolerance):
     """Test the importer OPERA HDF5."""
-    smart_assert(metadata_odyssey[variable], expected, tolerance)
+    smart_assert(variable, expected, tolerance)
 
 
 # list of (variable,expected,tolerance) tuples
 test_cirrus_attrs = [
-    ("projection", expected_proj, None),
-    ("ll_lon", -10.4345768386404, 1e-10),
-    ("ll_lat", 31.7462153182675, 1e-10),
-    ("ur_lon", 57.8119647501499, 1e-10),
-    ("ur_lat", 67.6210371071631, 1e-10),
-    ("x1", -0.00027143326587975025, 1e-6),
-    ("y1", -4400000.00116988, 1e-10),
-    ("x2", 3800000.0000817003, 1e-10),
-    ("y2", -8.761277422308922e-05, 1e-6),
-    ("xpixelsize", 1000.0, 1e-10),
-    ("ypixelsize", 1000.0, 1e-10),
-    ("cartesian_unit", "m", None),
-    ("accutime", 15.0, 1e-10),
-    ("yorigin", "upper", None),
-    ("unit", "dBZ", None),
-    ("institution", "Odyssey datacentre", None),
-    ("transform", "dB", None),
-    ("zerovalue", -32.0, 1e-10),
-    ("threshold", -31.5, 1e-10),
+    (precip_cirrus.attrs["projection"], expected_proj, None),
+    (float(precip_cirrus.lon.isel(x=0, y=0).values), -10.4306947565, 1e-10),
+    (float(precip_cirrus.lat.isel(x=0, y=0).values), 31.7518730135, 1e-10),
+    (float(precip_cirrus.lon.isel(x=-1, y=-1).values), 57.7949292793, 1e-10),
+    (float(precip_cirrus.lat.isel(x=-1, y=-1).values), 67.6207527344, 1e-10),
+    (precip_cirrus.x.isel(x=0).values, 499.99972864, 1e-6),
+    (precip_cirrus.y.isel(y=0).values, -4399500.00116976, 1e-10),
+    (precip_cirrus.x.isel(x=-1).values, 3799500.00025612, 1e-10),
+    (precip_cirrus.y.isel(y=-1).values, -500.00008774, 1e-6),
+    (precip_cirrus.x.attrs["stepsize"], 1000.0, 1e-10),
+    (precip_cirrus.y.attrs["stepsize"], 1000.0, 1e-10),
+    (precip_cirrus.x.attrs["units"], "m", None),
+    (precip_cirrus.y.attrs["units"], "m", None),
+    (precip_cirrus_dataarray.attrs["accutime"], 15.0, 1e-10),
+    # ("yorigin", "upper", None),
+    (precip_cirrus_dataarray.attrs["units"], "dBZ", None),
+    (precip_cirrus.attrs["institution"], "Odyssey datacentre", None),
+    (precip_cirrus_dataarray.attrs["transform"], "dB", None),
+    (precip_cirrus_dataarray.attrs["zerovalue"], -32.0, 1e-10),
+    (precip_cirrus_dataarray.attrs["threshold"], -31.5, 1e-10),
 ]
 
 
 @pytest.mark.parametrize("variable, expected, tolerance", test_cirrus_attrs)
 def test_io_import_opera_hdf5_cirrus_dataset_attrs(variable, expected, tolerance):
     """Test OPERA HDF5 importer: max. reflectivity composites from CIRRUS."""
-    smart_assert(metadata_cirrus[variable], expected, tolerance)
+    smart_assert(variable, expected, tolerance)
 
 
 # list of (variable,expected,tolerance) tuples
 test_nimbus_rain_rate_attrs = [
-    ("projection", expected_proj, None),
-    ("ll_lon", -10.434599999137568, 1e-10),
-    ("ll_lat", 31.74619995126678, 1e-10),
-    ("ur_lon", 57.8119032106317, 1e-10),
-    ("ur_lat", 67.62104536996274, 1e-10),
-    ("x1", -2.5302714337594807, 1e-6),
-    ("y1", -4400001.031169886, 1e-10),
-    ("x2", 3799997.4700817037, 1e-10),
-    ("y2", -1.0300876162946224, 1e-6),
-    ("xpixelsize", 2000.0, 1e-10),
-    ("ypixelsize", 2000.0, 1e-10),
-    ("cartesian_unit", "m", None),
-    ("accutime", 15.0, 1e-10),
-    ("yorigin", "upper", None),
-    ("unit", "mm/h", None),
-    ("institution", "Odyssey datacentre", None),
-    ("transform", None, None),
-    ("zerovalue", 0.0, 1e-10),
-    ("threshold", 0.01, 1e-10),
+    (precip_nimbus_rain_rate.attrs["projection"], expected_proj, None),
+    (float(precip_nimbus_rain_rate.lon.isel(x=0, y=0).values), -10.4268354001, 1e-10),
+    (float(precip_nimbus_rain_rate.lat.isel(x=0, y=0).values), 31.7575151437, 1e-10),
+    (float(precip_nimbus_rain_rate.lon.isel(x=-1, y=-1).values), 57.7778328845, 1e-10),
+    (float(precip_nimbus_rain_rate.lat.isel(x=-1, y=-1).values), 67.6204748496, 1e-10),
+    (precip_nimbus_rain_rate.x.isel(x=0).values, 997.46972871, 1e-6),
+    (precip_nimbus_rain_rate.y.isel(y=0).values, -4399001.03116964, 1e-10),
+    (precip_nimbus_rain_rate.x.isel(x=-1).values, 3798997.47025605, 1e-10),
+    (precip_nimbus_rain_rate.y.isel(y=-1).values, -1001.03008786, 1e-6),
+    (precip_nimbus_rain_rate.x.attrs["stepsize"], 2000.0, 1e-10),
+    (precip_nimbus_rain_rate.y.attrs["stepsize"], 2000.0, 1e-10),
+    (precip_nimbus_rain_rate.x.attrs["units"], "m", None),
+    (precip_nimbus_rain_rate.y.attrs["units"], "m", None),
+    (precip_nimbus_rain_rate_dataarray.attrs["accutime"], 15.0, 1e-10),
+    (precip_nimbus_rain_rate_dataarray.attrs["units"], "mm/h", None),
+    (precip_nimbus_rain_rate.attrs["institution"], "Odyssey datacentre", None),
+    (precip_nimbus_rain_rate_dataarray.attrs["transform"], None, None),
+    (precip_nimbus_rain_rate_dataarray.attrs["zerovalue"], 0.0, 1e-10),
+    (precip_nimbus_rain_rate_dataarray.attrs["threshold"], 0.01, 1e-10),
 ]
 
 
@@ -156,30 +164,30 @@ def test_io_import_opera_hdf5_nimbus_rain_rate_dataset_attrs(
     variable, expected, tolerance
 ):
     """Test OPERA HDF5 importer: rain rate composites from NIMBUS."""
-    smart_assert(metadata_nimbus_rain_rate[variable], expected, tolerance)
+    smart_assert(variable, expected, tolerance)
 
 
 # list of (variable,expected,tolerance) tuples
 test_nimbus_rain_accum_attrs = [
-    ("projection", expected_proj, None),
-    ("ll_lon", -10.434599999137568, 1e-10),
-    ("ll_lat", 31.74619995126678, 1e-10),
-    ("ur_lon", 57.8119032106317, 1e-10),
-    ("ur_lat", 67.62104536996274, 1e-10),
-    ("x1", -2.5302714337594807, 1e-6),
-    ("y1", -4400001.031169886, 1e-10),
-    ("x2", 3799997.4700817037, 1e-10),
-    ("y2", -1.0300876162946224, 1e-6),
-    ("xpixelsize", 2000.0, 1e-10),
-    ("ypixelsize", 2000.0, 1e-10),
-    ("cartesian_unit", "m", None),
-    ("accutime", 15.0, 1e-10),
-    ("yorigin", "upper", None),
-    ("unit", "mm", None),
-    ("institution", "Odyssey datacentre", None),
-    ("transform", None, None),
-    ("zerovalue", 0.0, 1e-10),
-    ("threshold", 0.01, 1e-10),
+    (precip_nimbus_rain_accum.attrs["projection"], expected_proj, None),
+    (float(precip_nimbus_rain_accum.lon.isel(x=0, y=0).values), -10.4268354001, 1e-10),
+    (float(precip_nimbus_rain_accum.lat.isel(x=0, y=0).values), 31.7575151437, 1e-10),
+    (float(precip_nimbus_rain_accum.lon.isel(x=-1, y=-1).values), 57.7778328845, 1e-10),
+    (float(precip_nimbus_rain_accum.lat.isel(x=-1, y=-1).values), 67.6204748496, 1e-10),
+    (precip_nimbus_rain_accum.x.isel(x=0).values, 997.46972871, 1e-6),
+    (precip_nimbus_rain_accum.y.isel(y=0).values, -4399001.03116964, 1e-10),
+    (precip_nimbus_rain_accum.x.isel(x=-1).values, 3798997.47025605, 1e-10),
+    (precip_nimbus_rain_accum.y.isel(y=-1).values, -1001.03008786, 1e-6),
+    (precip_nimbus_rain_accum.x.attrs["stepsize"], 2000.0, 1e-10),
+    (precip_nimbus_rain_accum.y.attrs["stepsize"], 2000.0, 1e-10),
+    (precip_nimbus_rain_accum.x.attrs["units"], "m", None),
+    (precip_nimbus_rain_accum.y.attrs["units"], "m", None),
+    (precip_nimbus_rain_accum_dataarray.attrs["accutime"], 15.0, 1e-10),
+    (precip_nimbus_rain_accum_dataarray.attrs["units"], "mm", None),
+    (precip_nimbus_rain_accum.attrs["institution"], "Odyssey datacentre", None),
+    (precip_nimbus_rain_accum_dataarray.attrs["transform"], None, None),
+    (precip_nimbus_rain_accum_dataarray.attrs["zerovalue"], 0.0, 1e-10),
+    (precip_nimbus_rain_accum_dataarray.attrs["threshold"], 0.01, 1e-10),
 ]
 
 
@@ -188,4 +196,4 @@ def test_io_import_opera_hdf5_nimbus_rain_accum_dataset_attrs(
     variable, expected, tolerance
 ):
     """Test OPERA HDF5 importer: rain accumulation composites from NIMBUS."""
-    smart_assert(metadata_nimbus_rain_accum[variable], expected, tolerance)
+    smart_assert(variable, expected, tolerance)
diff --git a/pysteps/tests/test_io_saf_crri.py b/pysteps/tests/test_io_saf_crri.py
index f9b9c249d..074e39f70 100644
--- a/pysteps/tests/test_io_saf_crri.py
+++ b/pysteps/tests/test_io_saf_crri.py
@@ -1,90 +1,97 @@
 # -*- coding: utf-8 -*-
 
-import os
-
 import pytest
 
-import pysteps
-from pysteps.tests.helpers import smart_assert
+from pysteps.tests.helpers import smart_assert, get_precipitation_fields
 
-pytest.importorskip("netCDF4")
+precip_dataset = get_precipitation_fields(
+    num_prev_files=0,
+    num_next_files=0,
+    return_raw=True,
+    metadata=True,
+    source="saf",
+    log_transform=False,
+)
 
+precip_var = precip_dataset.attrs["precip_var"]
+precip_dataarray = precip_dataset[precip_var]
 
-expected_proj = (
-    "+proj=geos +a=6378137.000000 +b=6356752.300000 "
-    "+lon_0=0.000000 +h=35785863.000000"
-)
-test_geodata_crri = [
-    ("projection", expected_proj, None),
-    ("x1", -3301500.0, 0.1),
-    ("x2", 3298500.0, 0.1),
-    ("y1", 2512500.0, 0.1),
-    ("y2", 5569500.0, 0.1),
-    ("xpixelsize", 3000.0, 0.1),
-    ("ypixelsize", 3000.0, 0.1),
-    ("cartesian_unit", "m", None),
-    ("yorigin", "upper", None),
+test_extent_crri = [
+    (None, (-3300000.0, 3297000.0, 2514000.0, 5568000.0), (1, 1019, 2200), None),
+    (
+        (-1980000.0, 1977000.0, 2514000.0, 4818000.0),
+        (-1977000.0, 1974000.0, 2517000.0, 4815000.0),
+        (1, 767, 1318),
+        None,
+    ),
 ]
 
 
-@pytest.mark.parametrize("variable, expected, tolerance", test_geodata_crri)
-def test_io_import_saf_crri_geodata(variable, expected, tolerance):
+@pytest.mark.parametrize(
+    "extent, expected_extent, expected_shape, tolerance", test_extent_crri
+)
+def test_io_import_saf_crri_extent(extent, expected_extent, expected_shape, 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"
+
+    precip_dataset_reduced_domain = get_precipitation_fields(
+        num_prev_files=0,
+        num_next_files=0,
+        return_raw=True,
+        metadata=True,
+        source="saf",
+        log_transform=False,
+        extent=extent,
     )
-    geodata = pysteps.io.importers._import_saf_crri_geodata(filename)
-    smart_assert(geodata[variable], expected, tolerance)
+    precip_var = precip_dataset_reduced_domain.attrs["precip_var"]
+    precip_dataarray_reduced_domain = precip_dataset_reduced_domain[precip_var]
+    x_min = float(precip_dataset_reduced_domain.x.isel(x=0).values)
+    x_max = float(precip_dataset_reduced_domain.x.isel(x=-1).values)
+    y_min = float(precip_dataset_reduced_domain.y.isel(y=0).values)
+    y_max = float(precip_dataset_reduced_domain.y.isel(y=-1).values)
+    extent_out = (x_min, x_max, y_min, y_max)
+    smart_assert(extent_out, expected_extent, tolerance)
+    smart_assert(precip_dataarray_reduced_domain.shape, expected_shape, tolerance)
 
 
-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"
+expected_proj = (
+    "+proj=geos +a=6378137.000000 +b=6356752.300000 "
+    "+lon_0=0.000000 +h=35785863.000000"
 )
-_, _, 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),
+test_geodata_crri = [
+    (precip_dataset.attrs["projection"], expected_proj, None),
+    (precip_dataset.x.isel(x=0).values, -3300000.0, 1e-6),
+    (precip_dataset.y.isel(y=0).values, 2514000.0, 1e-6),
+    (precip_dataset.x.isel(x=-1).values, 3297000.0, 1e-6),
+    (precip_dataset.y.isel(y=-1).values, 5568000.0, 1e-9),
+    (precip_dataset.x.attrs["stepsize"], 3000.0, 1e-10),
+    (precip_dataset.y.attrs["stepsize"], 3000.0, 1e-10),
+    (precip_dataset.x.attrs["units"], "m", None),
+    (precip_dataset.y.attrs["units"], "m", None),
 ]
 
 
-@pytest.mark.parametrize("variable, expected, tolerance", test_attrs)
-def test_io_import_saf_crri_attrs(variable, expected, tolerance):
-    """Test the importer SAF CRRI."""
-    smart_assert(metadata[variable], expected, tolerance)
+@pytest.mark.parametrize("variable, expected, tolerance", test_geodata_crri)
+def test_io_import_saf_crri_geodata(variable, expected, tolerance):
+    smart_assert(variable, expected, tolerance)
 
 
-test_extent_crri = [
-    (None, (-3301500.0, 3298500.0, 2512500.0, 5569500.0), (1019, 2200), None),
+# list of (variable,expected,tolerance) tuples
+test_attrs = [
+    (precip_dataset.attrs["projection"], expected_proj, None),
     (
-        (-1980000.0, 1977000.0, 2514000.0, 4818000.0),
-        (-1978500.0, 1975500.0, 2515500.0, 4816500.0),
-        (767, 1318),
+        precip_dataset.attrs["institution"],
+        "Agencia Estatal de Meteorología (AEMET)",
         None,
     ),
+    (precip_dataarray.attrs["accutime"], None, None),
+    (precip_dataarray.attrs["units"], "mm/h", None),
+    (precip_dataarray.attrs["transform"], None, None),
+    (precip_dataarray.attrs["zerovalue"], 0.0, 1e-6),
 ]
 
 
-@pytest.mark.parametrize(
-    "extent, expected_extent, expected_shape, tolerance", test_extent_crri
-)
-def test_io_import_saf_crri_extent(extent, expected_extent, expected_shape, 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"
-    )
-    precip, _, metadata = pysteps.io.import_saf_crri(filename, extent=extent)
-    extent_out = (metadata["x1"], metadata["x2"], metadata["y1"], metadata["y2"])
-    smart_assert(extent_out, expected_extent, tolerance)
-    smart_assert(precip.shape, expected_shape, tolerance)
+    smart_assert(variable, expected, tolerance)
diff --git a/pysteps/tests/test_noise_fftgenerators.py b/pysteps/tests/test_noise_fftgenerators.py
index cecaf8ca4..cc27258e1 100644
--- a/pysteps/tests/test_noise_fftgenerators.py
+++ b/pysteps/tests/test_noise_fftgenerators.py
@@ -4,18 +4,23 @@
 from pysteps.tests.helpers import get_precipitation_fields
 
 
-PRECIP = get_precipitation_fields(
+precip_dataset = get_precipitation_fields(
     num_prev_files=0,
     num_next_files=0,
     return_raw=False,
     metadata=False,
     upscale=2000,
 )
-PRECIP = PRECIP.filled()
 
+precip_var = precip_dataset.attrs["precip_var"]
+precip_dataarray = precip_dataset[precip_var]
 
+
+# XR: all tests assume a 2D field, so we select the first timestep, these tests need to be changed when fftgenerators support xarray DataArrays
 def test_noise_param_2d_fft_filter():
-    fft_filter = fftgenerators.initialize_param_2d_fft_filter(PRECIP)
+    fft_filter = fftgenerators.initialize_param_2d_fft_filter(
+        precip_dataarray.isel(time=0).values
+    )
 
     assert isinstance(fft_filter, dict)
     assert all([key in fft_filter for key in ["field", "input_shape", "model", "pars"]])
@@ -23,11 +28,13 @@ def test_noise_param_2d_fft_filter():
     out = fftgenerators.generate_noise_2d_fft_filter(fft_filter)
 
     assert isinstance(out, np.ndarray)
-    assert out.shape == PRECIP.shape
+    assert out.shape == precip_dataarray.isel(time=0).shape
 
 
 def test_noise_nonparam_2d_fft_filter():
-    fft_filter = fftgenerators.initialize_nonparam_2d_fft_filter(PRECIP)
+    fft_filter = fftgenerators.initialize_nonparam_2d_fft_filter(
+        precip_dataarray.isel(time=0).values
+    )
 
     assert isinstance(fft_filter, dict)
     assert all([key in fft_filter for key in ["field", "input_shape"]])
@@ -35,11 +42,13 @@ def test_noise_nonparam_2d_fft_filter():
     out = fftgenerators.generate_noise_2d_fft_filter(fft_filter)
 
     assert isinstance(out, np.ndarray)
-    assert out.shape == PRECIP.shape
+    assert out.shape == precip_dataarray.isel(time=0).shape
 
 
 def test_noise_nonparam_2d_ssft_filter():
-    fft_filter = fftgenerators.initialize_nonparam_2d_ssft_filter(PRECIP)
+    fft_filter = fftgenerators.initialize_nonparam_2d_ssft_filter(
+        precip_dataarray.isel(time=0).values
+    )
 
     assert isinstance(fft_filter, dict)
     assert all([key in fft_filter for key in ["field", "input_shape"]])
@@ -47,11 +56,13 @@ def test_noise_nonparam_2d_ssft_filter():
     out = fftgenerators.generate_noise_2d_ssft_filter(fft_filter)
 
     assert isinstance(out, np.ndarray)
-    assert out.shape == PRECIP.shape
+    assert out.shape == precip_dataarray.isel(time=0).shape
 
 
 def test_noise_nonparam_2d_nested_filter():
-    fft_filter = fftgenerators.initialize_nonparam_2d_nested_filter(PRECIP)
+    fft_filter = fftgenerators.initialize_nonparam_2d_nested_filter(
+        precip_dataarray.isel(time=0).values
+    )
 
     assert isinstance(fft_filter, dict)
     assert all([key in fft_filter for key in ["field", "input_shape"]])
@@ -59,4 +70,4 @@ def test_noise_nonparam_2d_nested_filter():
     out = fftgenerators.generate_noise_2d_ssft_filter(fft_filter)
 
     assert isinstance(out, np.ndarray)
-    assert out.shape == PRECIP.shape
+    assert out.shape == precip_dataarray.isel(time=0).shape
diff --git a/pysteps/tests/test_utils_reprojection.py b/pysteps/tests/test_utils_reprojection.py
index 84b0f177b..2fcb4f7aa 100644
--- a/pysteps/tests/test_utils_reprojection.py
+++ b/pysteps/tests/test_utils_reprojection.py
@@ -3,23 +3,99 @@
 import os
 import numpy as np
 import pytest
-import pysteps
+import xarray as xr
+
 from pysteps.utils import reprojection as rpj
+from pysteps.tests.helpers import get_precipitation_fields
+
+
+def build_precip_dataset(
+    data: np.ndarray,  # shape (time, y, x)
+    *,
+    projection: str = "EPSG:3035",  # PROJ4/EPSG string
+    cartesian_unit: str = "m",  # 'm' or 'km'
+    institution: str = "rmi",
+    precip_var_name: str = "precip_intensity",  # or 'precip_accum' / 'reflectivity'
+    # grid + time spec (regular spacing)
+    nx: int | None = None,
+    ny: int | None = None,
+    dx: float = 1000.0,  # x stepsize, in cartesian_unit
+    dy: float = 1000.0,  # y stepsize, in cartesian_unit
+    x0: float | None = None,  # center of pixel (0,0) in cartesian_unit
+    y0: float | None = None,  # center of pixel (0,0) in cartesian_unit
+    nt: int | None = None,
+    t0: int = 0,  # seconds since forecast start
+    dt: int = 3600,  # timestep in seconds
+    # precip variable attrs
+    units: str = "mm/h",
+    transform_attr: str | None = None,  # e.g. 'dB', 'Box-Cox', or None
+    accutime_min: float = 60.0,
+    threshold: float = 0.1,
+    zerovalue: float = 0.0,
+    zr_a: float = 200.0,
+    zr_b: float = 1.6,
+) -> xr.Dataset:
+    assert data.ndim == 3, "data must be (time, y, x)"
+    nt_ = data.shape[0] if nt is None else nt
+    ny_ = data.shape[1] if ny is None else ny
+    nx_ = data.shape[2] if nx is None else nx
+
+    # Build regular coords (centers). If x0/y0 are not given, start at half a pixel.
+    if x0 is None:
+        x0 = 0.5 * dx
+    if y0 is None:
+        y0 = 0.5 * dy
+
+    x = x0 + dx * np.arange(nx_)
+    y = y0 + dy * np.arange(ny_)
+    time = t0 + dt * np.arange(nt_)
+
+    da = xr.DataArray(
+        data,
+        dims=("time", "y", "x"),
+        coords={"time": time, "y": y, "x": x},
+        name=precip_var_name,
+        attrs={
+            "units": units,
+            "transform": transform_attr,
+            "accutime": float(accutime_min),
+            "threshold": float(threshold),
+            "zerovalue": float(zerovalue),
+            "zr_a": float(zr_a),
+            "zr_b": float(zr_b),
+        },
+    )
 
-pytest.importorskip("rasterio")
-pytest.importorskip("pyproj")
+    # stepsize attrs on coords (required by your spec)
+    da.coords["time"].attrs["stepsize"] = int(dt)  # seconds
+    da.coords["time"].attrs["standard_name"] = "time"
+    da.coords["x"].attrs["stepsize"] = float(dx)  # in cartesian_unit
+    da.coords["x"].attrs["units"] = cartesian_unit
+    da.coords["y"].attrs["stepsize"] = float(dy)  # in cartesian_unit
+    da.coords["y"].attrs["units"] = cartesian_unit
+
+    ds = xr.Dataset({precip_var_name: da})
+    ds.attrs.update(
+        {
+            "projection": projection,  # PROJ string or EPSG code
+            "institution": institution,
+            "precip_var": precip_var_name,
+            "cartesian_unit": cartesian_unit,
+        }
+    )
 
-root_path_radar = pysteps.rcparams.data_sources["rmi"]["root_path"]
+    return ds
 
-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"
+precip_dataset = get_precipitation_fields(
+    num_prev_files=0,
+    num_next_files=0,
+    source="rmi",
+    return_raw=True,
+    metadata=True,
+    log_transform=False,
 )
 
-# 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))
 
@@ -42,72 +118,64 @@
     "+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",
+nwp_dataset = build_precip_dataset(
+    nwp_array,
+    projection=nwp_proj,  # ETRS89 / LAEA Europe (meters)
     cartesian_unit="m",
-    x1=0.0,
-    x2=731900.0,
-    y1=-731900.0,
-    y2=0.0,
+    precip_var_name="precip_intensity",
+    dx=1300.0,
+    dy=1300.0,  # 1 km grid
+    dt=300,  # hourly
+    accutime_min=5.0,  # accumulation window (min)
+    threshold=0.1,  # mm/h rain/no-rain threshold
+    zerovalue=0.0,
 )
 
 steps_arg_names = (
-    "radar_array",
-    "nwp_array",
-    "metadata_src",
-    "metadata_dst",
+    "precip_dataset",
+    "nwp_dataset",
 )
 
 steps_arg_values = [
-    (radar_array, nwp_array, metadata_src, metadata_dst),
+    (precip_dataset, nwp_dataset),
 ]
 
 
+# XR: since reproject_grids is not xarray compatible yet, we cannot use xarray DataArrays in the tests
 @pytest.mark.parametrize(steps_arg_names, steps_arg_values)
-def test_utils_reproject_grids(
-    radar_array,
-    nwp_array,
-    metadata_src,
-    metadata_dst,
-):
+def test_utils_reproject_grids(precip_dataset, nwp_dataset):
     # Reproject
-    nwp_array_reproj, metadata_reproj = rpj.reproject_grids(
-        nwp_array, radar_array, metadata_src, metadata_dst
-    )
+    nwp_dataset_reproj = rpj.reproject_grids(nwp_dataset, precip_dataset)
 
+    nwp_dataset_reproj_dataarray = nwp_dataset_reproj[
+        nwp_dataset_reproj.attrs["precip_var"]
+    ]
+    nwp_dataarray = nwp_dataset[nwp_dataset.attrs["precip_var"]]
+    precip_dataarray = precip_dataset[precip_dataset.attrs["precip_var"]]
     # The tests
     assert (
-        nwp_array_reproj.shape[0] == nwp_array.shape[0]
+        nwp_dataset_reproj_dataarray.shape[0] == nwp_dataarray.shape[0]
     ), "Time dimension has not the same length as source"
     assert (
-        nwp_array_reproj.shape[1] == radar_array.shape[0]
+        nwp_dataset_reproj_dataarray.shape[1] == precip_dataarray.shape[1]
     ), "y dimension has not the same length as radar composite"
     assert (
-        nwp_array_reproj.shape[2] == radar_array.shape[1]
+        nwp_dataset_reproj_dataarray.shape[2] == precip_dataarray.shape[2]
     ), "x dimension has not the same length as radar composite"
 
-    assert (
-        metadata_reproj["x1"] == metadata_dst["x1"]
+    assert float(nwp_dataset_reproj_dataarray.x.isel(x=0).values) == float(
+        precip_dataarray.x.isel(x=0).values
     ), "x-value lower left corner is not equal to radar composite"
-    assert (
-        metadata_reproj["x2"] == metadata_dst["x2"]
+    assert float(nwp_dataset_reproj_dataarray.x.isel(x=-1).values) == float(
+        precip_dataarray.x.isel(x=-1).values
     ), "x-value upper right corner is not equal to radar composite"
-    assert (
-        metadata_reproj["y1"] == metadata_dst["y1"]
+    assert float(nwp_dataset_reproj_dataarray.y.isel(y=0).values) == float(
+        precip_dataarray.y.isel(y=0).values
     ), "y-value lower left corner is not equal to radar composite"
-    assert (
-        metadata_reproj["y2"] == metadata_dst["y2"]
+    assert float(nwp_dataset_reproj_dataarray.y.isel(y=-1).values) == float(
+        precip_dataarray.y.isel(y=-1).values
     ), "y-value upper right corner is not equal to radar composite"
 
     assert (
-        metadata_reproj["projection"] == metadata_dst["projection"]
+        nwp_dataset_reproj.attrs["projection"] == precip_dataset.attrs["projection"]
     ), "projection is different than destination projection"
diff --git a/pysteps/tests/test_verification_probscores.py b/pysteps/tests/test_verification_probscores.py
index c7f9990b8..4fcb8451a 100644
--- a/pysteps/tests/test_verification_probscores.py
+++ b/pysteps/tests/test_verification_probscores.py
@@ -8,10 +8,21 @@
 from pysteps.tests.helpers import get_precipitation_fields
 from pysteps.verification import probscores
 
-precip = get_precipitation_fields(num_next_files=10, return_raw=True)
+precip_dataset = get_precipitation_fields(num_next_files=10, return_raw=True)
+
+precip_var = precip_dataset.attrs["precip_var"]
+precip_dataarray = precip_dataset[precip_var]
+
+# XR: the scorring code has not been made xarray compatible, so we need to convert to numpy arrays. Once changed we can properly test these scores with xarray DataArrays
 
 # CRPS
-test_data = [(precip[:10], precip[-1], 0.01470871)]
+test_data = [
+    (
+        precip_dataarray.isel(time=slice(0, 10)).values,
+        precip_dataarray.isel(time=-1).values,
+        0.01470871,
+    )
+]
 
 
 @pytest.mark.parametrize("X_f, X_o, expected", test_data)
@@ -21,7 +32,16 @@ def test_CRPS(X_f, X_o, expected):
 
 
 # reldiag
-test_data = [(precip[:10], precip[-1], 1.0, 10, 10, 3.38751492)]
+test_data = [
+    (
+        precip_dataarray.isel(time=slice(0, 10)).values,
+        precip_dataarray.isel(time=-1).values,
+        1.0,
+        10,
+        10,
+        3.38751492,
+    )
+]
 
 
 @pytest.mark.parametrize("X_f, X_o, X_min, n_bins, min_count, expected", test_data)
@@ -34,7 +54,16 @@ def test_reldiag_sum(X_f, X_o, X_min, n_bins, min_count, expected):
 
 
 # ROC_curve
-test_data = [(precip[:10], precip[-1], 1.0, 10, True, 0.79557329)]
+test_data = [
+    (
+        precip_dataarray.isel(time=slice(0, 10)).values,
+        precip_dataarray.isel(time=-1).values,
+        1.0,
+        10,
+        True,
+        0.79557329,
+    )
+]
 
 
 @pytest.mark.parametrize(
diff --git a/pysteps/tests/test_verification_spatialscores.py b/pysteps/tests/test_verification_spatialscores.py
index a02bc0773..c7ca1a70c 100644
--- a/pysteps/tests/test_verification_spatialscores.py
+++ b/pysteps/tests/test_verification_spatialscores.py
@@ -6,10 +6,32 @@
 from pysteps.tests.helpers import get_precipitation_fields
 from pysteps.verification import spatialscores
 
-R = get_precipitation_fields(num_prev_files=1, return_raw=True)
+precip_dataset = get_precipitation_fields(num_prev_files=1, return_raw=True)
+
+precip_var = precip_dataset.attrs["precip_var"]
+precip_dataarray = precip_dataset[precip_var]
+
+# XR: the scorring code has not been made xarray compatible, so we need to convert to numpy arrays. Once changed we can properly test these scores with xarray DataArrays
+# BUG: the tests for BMSE below reverse the arrays with [::-1], this should be fixed in the scoring code
 test_data = [
-    (R[0], R[1], "FSS", [1], [10], None, 0.85161531),
-    (R[0], R[1], "BMSE", [1], None, "Haar", 0.99989651),
+    (
+        precip_dataarray.isel(time=0).values,
+        precip_dataarray.isel(time=1).values,
+        "FSS",
+        [1],
+        [10],
+        None,
+        0.85161531,
+    ),
+    (
+        precip_dataarray.isel(time=0).values[::-1],
+        precip_dataarray.isel(time=1).values[::-1],
+        "BMSE",
+        [1],
+        None,
+        "Haar",
+        0.99989651,
+    ),
 ]
 
 
@@ -25,10 +47,36 @@ def test_intensity_scale(X_f, X_o, name, thrs, scales, wavelet, expected):
     )
 
 
-R = get_precipitation_fields(num_prev_files=3, return_raw=True)
+precip_dataset = get_precipitation_fields(num_next_files=3, return_raw=True)
+
+precip_var = precip_dataset.attrs["precip_var"]
+precip_dataarray = precip_dataset[precip_var]
+
 test_data = [
-    (R[:2], R[2:], "FSS", [1], [10], None),
-    (R[:2], R[2:], "BMSE", [1], None, "Haar"),
+    (
+        precip_dataarray.isel(time=slice(0, 2)).values,
+        precip_dataarray.isel(
+            time=slice(
+                2,
+            )
+        ).values,
+        "FSS",
+        [1],
+        [10],
+        None,
+    ),
+    (
+        precip_dataarray.isel(time=slice(0, 2)).values,
+        precip_dataarray.isel(
+            time=slice(
+                2,
+            )
+        ).values,
+        "BMSE",
+        [1],
+        None,
+        "Haar",
+    ),
 ]
 
 
diff --git a/pysteps/utils/conversion.py b/pysteps/utils/conversion.py
index 2ea6a3a12..92e08a450 100644
--- a/pysteps/utils/conversion.py
+++ b/pysteps/utils/conversion.py
@@ -28,12 +28,12 @@
 def cf_parameters_from_unit(unit: str) -> tuple[str, dict[str, str | None]]:
     if unit == "mm/h":
         var_name = "precip_intensity"
-        var_standard_name = None
+        var_standard_name = "instantaneous_precipitation_rate"
         var_long_name = "instantaneous precipitation rate"
         var_unit = "mm/h"
     elif unit == "mm":
         var_name = "precip_accum"
-        var_standard_name = None
+        var_standard_name = "accumulated_precipitation"
         var_long_name = "accumulated precipitation"
         var_unit = "mm"
     elif unit == "dBZ":
diff --git a/pysteps/utils/reprojection.py b/pysteps/utils/reprojection.py
index 6144a52c0..10cf0dc0c 100644
--- a/pysteps/utils/reprojection.py
+++ b/pysteps/utils/reprojection.py
@@ -14,107 +14,69 @@
 from pysteps.exceptions import MissingOptionalDependency
 
 import numpy as np
+import xarray as xr
 
 try:
-    from rasterio import Affine as A
-    from rasterio.warp import reproject, Resampling
+    import pyproj
 
-    RASTERIO_IMPORTED = True
+    PYPROJ_IMPORTED = True
 except ImportError:
-    RASTERIO_IMPORTED = False
+    PYPROJ_IMPORTED = False
 
 
-def reproject_grids(src_array, dst_array, metadata_src, metadata_dst):
+def reproject_grids(src_dataset, dst_dataset):
     """
     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 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.
+    src_dataset: xr.Dataset
+        xr.Dataset containing a precipitation variable which needs to be reprojected
+    dst_dataset: xr.Dataset
+        xr.Dataset containing a precipitation variable which is used to project the provided src_dataset
 
     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.
+    reprojected_dataset: xr.Dataset
+        xr.Dataset containing the reprojected precipitation variable
     """
 
-    if not RASTERIO_IMPORTED:
+    if not PYPROJ_IMPORTED:
         raise MissingOptionalDependency(
-            "rasterio package is required for the reprojection module, but it is "
+            "pyproj 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)
+    x_r = dst_dataset.x.values
+    y_r = dst_dataset.y.values
+    x_2d, y_2d = np.meshgrid(x_r, y_r)
+    # Calculate match  between the two projections
+    transfomer = pyproj.Transformer.from_proj(
+        src_dataset.attrs["projection"], dst_dataset.attrs["projection"]
     )
-
-    # 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)
+    dest_src_x, dest_src_y = transfomer.transform(
+        x_2d.flatten(), y_2d.flatten(), direction="INVERSE"
     )
-
-    # 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
+    dest_src_x, dest_src_y = dest_src_x.reshape(x_2d.shape), dest_src_y.reshape(
+        y_2d.shape
+    )
+    dest_src_x_dataarray = xr.DataArray(
+        dest_src_x, dims=("y_src", "x_src"), coords={"y_src": y_r, "x_src": x_r}
+    )
+    dest_src_y_dataarray = xr.DataArray(
+        dest_src_y, dims=("y_src", "x_src"), coords={"y_src": y_r, "x_src": x_r}
+    )
+    # Select the nearest neighbour in the source dataset for each point in the destination dataset
+    reproj_dataset = src_dataset.sel(
+        x=dest_src_x_dataarray, y=dest_src_y_dataarray, method="nearest"
+    )
+    # Clean up the dataset
+    reproj_dataset = reproj_dataset.drop_vars(["x", "y"])
+    reproj_dataset = reproj_dataset.rename({"x_src": "x", "y_src": "y"})
+    # Fill attributes from dst_dataset to reproj_dataset
+    reproj_dataset.attrs = dst_dataset.attrs
+    reproj_dataset[reproj_dataset.attrs["precip_var"]].attrs = dst_dataset[
+        dst_dataset.attrs["precip_var"]
+    ].attrs
+
+    return reproj_dataset