diff --git a/.github/workflows/test_pysteps.yml b/.github/workflows/test_pysteps.yml
index 08a5e5061..ecb0db7b9 100644
--- a/.github/workflows/test_pysteps.yml
+++ b/.github/workflows/test_pysteps.yml
@@ -2,7 +2,7 @@ name: Test Pysteps
 
 env:
   MINIMAL_DEPENDENCIES: cython numpy jsmin jsonschema matplotlib netCDF4 opencv pillow pyproj scipy dask
-  OPTIONAL_DEPENDENCIES: pyfftw cartopy h5py PyWavelets pandas scikit-image rasterio
+  OPTIONAL_DEPENDENCIES: pyfftw cartopy h5py PyWavelets pandas scikit-image
   TEST_DEPENDENCIES: pytest pytest-cov
 
 on:
diff --git a/doc/requirements.txt b/doc/requirements.txt
index cd7856d7e..42eda9fad 100644
--- a/doc/requirements.txt
+++ b/doc/requirements.txt
@@ -1,4 +1,4 @@
-# Additional requirements related to the documentation build only
+# Additional requeriments related to the documentation build only
 sphinx
 numpydoc
 sphinxcontrib.bibtex
diff --git a/doc/source/pysteps_reference/utils.rst b/doc/source/pysteps_reference/utils.rst
index 0622e6bbe..7e5474145 100644
--- a/doc/source/pysteps_reference/utils.rst
+++ b/doc/source/pysteps_reference/utils.rst
@@ -16,4 +16,3 @@ Implementation of miscellaneous utility functions.
 .. automodule:: pysteps.utils.spectral
 .. automodule:: pysteps.utils.tapering
 .. automodule:: pysteps.utils.transformation
-.. automodule:: pysteps.utils.reprojection
diff --git a/environment_dev.yml b/environment_dev.yml
index 337a567d3..34b894cd4 100644
--- a/environment_dev.yml
+++ b/environment_dev.yml
@@ -31,4 +31,3 @@ dependencies:
   - scikit-image
   - pandas
   - xarray
-  - rasterio
diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py
index 8c19c1990..4b00ac9b9 100644
--- a/pysteps/io/importers.py
+++ b/pysteps/io/importers.py
@@ -917,7 +917,7 @@ def _import_fmi_pgm_metadata(filename, gzipped=False):
 
 
 @postprocess_import()
-def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1000.0, **kwargs):
+def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1.0, **kwargs):
     """Import a precipitation or reflectivity field (and optionally the quality
     field) from a HDF5 file conforming to the KNMI Data Centre specification.
 
@@ -934,8 +934,8 @@ def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1000.0, **kwa
         is used as default, but hourly, daily and monthly accumulations
         are also available.
     pixelsize: float
-        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
+        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
         size are also available.
 
     {extra_kwargs_doc}
@@ -1013,7 +1013,7 @@ def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1000.0, **kwa
     # The 'where' group of mch- and Opera-data, is called 'geographic' in the
     # KNMI data.
     geographic = f["geographic"]
-    proj4str = "+proj=stere +lat_0=90 +lon_0=0.0 +lat_ts=60.0 +a=6378137 +b=6356752 +x_0=0 +y_0=0"
+    proj4str = geographic["map_projection"].attrs["projection_proj4_params"].decode()
     pr = pyproj.Proj(proj4str)
     metadata["projection"] = proj4str
 
@@ -1044,7 +1044,7 @@ def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1000.0, **kwa
     metadata["y2"] = y2
     metadata["xpixelsize"] = pixelsize
     metadata["ypixelsize"] = pixelsize
-    metadata["cartesian_unit"] = "m"
+    metadata["cartesian_unit"] = "km"
     metadata["yorigin"] = "upper"
     metadata["institution"] = "KNMI - Royal Netherlands Meteorological Institute"
     metadata["accutime"] = accutime
diff --git a/pysteps/io/nwp_importers.py b/pysteps/io/nwp_importers.py
index a507dae1e..f27e8dbb4 100644
--- a/pysteps/io/nwp_importers.py
+++ b/pysteps/io/nwp_importers.py
@@ -371,20 +371,8 @@ def _import_rmi_nwp_geodata_xr(
     xpixelsize = abs(ds_in.x[1] - ds_in.x[0])
     ypixelsize = abs(ds_in.y[1] - ds_in.y[0])
 
-    x_coords = (
-        np.arange(xmin, xmin + float(xpixelsize) * ds_in.x.shape[0], float(xpixelsize))
-        + float(xpixelsize) / 2
-    )
-    y_coords = (
-        np.arange(ymin, ymin + float(ypixelsize) * ds_in.y.shape[0], float(ypixelsize))
-        + float(ypixelsize) / 2
-    )
-
     cartesian_unit = ds_in.x.units
 
-    # Assign new coordinates that are in the middle of the cell
-    ds_in = ds_in.assign_coords(x=x_coords, y=y_coords)
-
     # Add metadata needed by pySTEPS as attrs in X and Y variables
 
     ds_in.x.attrs.update(
@@ -529,28 +517,16 @@ def _import_knmi_nwp_geodata_xr(
     xpixelsize = abs(ds_in.x[1] - ds_in.x[0])
     ypixelsize = abs(ds_in.y[1] - ds_in.y[0])
 
-    x_coords = (
-        np.arange(xmin, xmin + float(xpixelsize) * ds_in.x.shape[0], float(xpixelsize))
-        + float(xpixelsize) / 2
-    )
-    y_coords = (
-        np.arange(ymin, ymin + float(ypixelsize) * ds_in.y.shape[0], float(ypixelsize))
-        + float(ypixelsize) / 2
-    )
-
-    cartesian_unit_x = ds_in.x.units
-    cartesian_unit_y = ds_in.y.units
-
-    # Assign new coordinates that are in the middle of the cell
-    ds_in = ds_in.assign_coords(x=x_coords, y=y_coords)
+    cartesian_unit = ds_in.x.units
 
     # Add metadata needed by pySTEPS as attrs in X and Y variables
+
     ds_in.x.attrs.update(
         {
             # TODO: Remove before final 2.0 version
             "x1": xmin,
             "x2": xmax,
-            "cartesian_unit": cartesian_unit_x,
+            "cartesian_unit": cartesian_unit,
         }
     )
 
@@ -559,14 +535,12 @@ def _import_knmi_nwp_geodata_xr(
             # TODO: Remove before final 2.0 version
             "y1": ymin,
             "y2": ymax,
-            "cartesian_unit": cartesian_unit_y,
+            "cartesian_unit": cartesian_unit,
         }
     )
 
     # Add metadata needed by pySTEPS as attrs in rainfall variable
     da_rainfall = ds_in[varname].isel({varname_time: 0})
-    # Set values below 0.0 to 0.0
-    da_rainfall = da_rainfall.where(da_rainfall >= 0.0, 0.0)
 
     ds_in[varname].attrs.update(
         {
diff --git a/pysteps/pystepsrc b/pysteps/pystepsrc
index 24b133c88..2f91010c6 100644
--- a/pysteps/pystepsrc
+++ b/pysteps/pystepsrc
@@ -4,8 +4,7 @@
     "silent_import": false,
     "outputs": {
         // path_outputs : path where to save results (figures, forecasts, etc)
-        "path_outputs": "./",
-        "path_workdir": "./tmp/"
+        "path_outputs": "./"
     },
     "plot": {
         // "motion_plot" : "streamplot" or "quiver"
diff --git a/pysteps/tests/test_io_bom_nwp.py b/pysteps/tests/test_io_bom_nwp.py
index d7deae63f..a3ad71301 100644
--- a/pysteps/tests/test_io_bom_nwp.py
+++ b/pysteps/tests/test_io_bom_nwp.py
@@ -53,7 +53,7 @@ def test_io_import_bom_nwp_xarray_attrs(variable, expected, tolerance):
 
 
 test_attrs_xr_coord_x = [
-    ("cartesian_unit", "m", None),
+    ("units", "m", None),
     ("x1", -127750.0, 0.1),
     ("x2", 127750.0, 0.1),
 ]
@@ -66,7 +66,7 @@ def test_io_import_bom_nwp_xarray_attrs_coordx(variable, expected, tolerance):
 
 
 test_attrs_xr_coord_y = [
-    ("cartesian_unit", "m", None),
+    ("units", "m", None),
     ("y1", -127750.0, 0.1),
     ("y2", 127750.0, 0.1),
 ]
diff --git a/pysteps/tests/test_io_knmi_hdf5.py b/pysteps/tests/test_io_knmi_hdf5.py
index aee3076d0..657e522bf 100644
--- a/pysteps/tests/test_io_knmi_hdf5.py
+++ b/pysteps/tests/test_io_knmi_hdf5.py
@@ -25,7 +25,9 @@ def test_io_import_knmi_hdf5_shape():
 # test_metadata: list of (variable,expected, tolerance) tuples
 
 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"
+    "+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"
 )
 
 # list of (variable,expected,tolerance) tuples
diff --git a/pysteps/tests/test_io_knmi_nwp.py b/pysteps/tests/test_io_knmi_nwp.py
index 95c332dcd..2bb61b54b 100644
--- a/pysteps/tests/test_io_knmi_nwp.py
+++ b/pysteps/tests/test_io_knmi_nwp.py
@@ -53,7 +53,7 @@ def test_io_import_knmi_nwp_xarray_attrs(variable, expected, tolerance):
 
 
 test_attrs_xr_coord_x = [
-    ("cartesian_unit", "degrees_east", None),
+    ("units", "degrees_east", None),
     ("x1", 0.0, 0.0001),
     ("x2", 11.063, 0.0001),
 ]
@@ -66,7 +66,7 @@ def test_io_import_knmi_nwp_xarray_attrs_coordx(variable, expected, tolerance):
 
 
 test_attrs_xr_coord_y = [
-    ("cartesian_unit", "degrees_north", None),
+    ("units", "degrees_north", None),
     ("y1", 49.0, 0.0001),
     ("y2", 55.877, 0.0001),
 ]
diff --git a/pysteps/tests/test_io_rmi_nwp.py b/pysteps/tests/test_io_rmi_nwp.py
index e4ab07bed..ef1689dae 100644
--- a/pysteps/tests/test_io_rmi_nwp.py
+++ b/pysteps/tests/test_io_rmi_nwp.py
@@ -56,7 +56,7 @@ def test_io_import_rmi_nwp_xarray_attrs(variable, expected, tolerance):
 
 
 test_attrs_xr_coord_x = [
-    ("cartesian_unit", "m", None),
+    ("units", "m", None),
     ("x1", 0, 0.1),
     ("x2", 731900.0, 0.1),
 ]
@@ -69,7 +69,7 @@ def test_io_import_rmi_nwp_xarray_attrs_coordx(variable, expected, tolerance):
 
 
 test_attrs_xr_coord_y = [
-    ("cartesian_unit", "m", None),
+    ("units", "m", None),
     ("y1", -731900.0, 0.1),
     ("y2", 0.0, 0.1),
 ]
diff --git a/pysteps/tests/test_utils_reprojection.py b/pysteps/tests/test_utils_reprojection.py
deleted file mode 100644
index 7bc5c9b10..000000000
--- a/pysteps/tests/test_utils_reprojection.py
+++ /dev/null
@@ -1,92 +0,0 @@
-# -*- coding: utf-8 -*-
-
-import os
-import xarray as xr
-import pytest
-import pysteps
-from pysteps.utils import reprojection
-
-pytest.importorskip("rasterio")
-
-root_path_radar = pysteps.rcparams.data_sources["knmi"]["root_path"]
-root_path_nwp = pysteps.rcparams.data_sources["knmi_nwp"]["root_path"]
-
-rel_path_nwp = os.path.join("2018", "09", "05")
-rel_path_radar = os.path.join(
-    "2010", "08"
-)  # Different day, but that does not matter for the tester
-
-filename_nwp = os.path.join(
-    root_path_nwp, rel_path_nwp, "20180905_0600_Pforecast_Harmonie.nc"
-)
-filename_radar = os.path.join(
-    root_path_radar, rel_path_radar, "RAD_NL25_RAP_5min_201008260230.h5"
-)
-
-# Open the radar and NWP data
-radar_array_xr = pysteps.io.importers.import_knmi_hdf5(filename_radar)
-nwp_array_xr = pysteps.io.import_knmi_nwp_xr(filename_nwp)
-
-steps_arg_names = (
-    "radar_array_xr",
-    "nwp_array_xr",
-    "t",
-)
-
-steps_arg_values = [
-    (radar_array_xr, nwp_array_xr, 1),
-    (radar_array_xr, nwp_array_xr, 4),
-    (radar_array_xr, nwp_array_xr, 8),
-]
-
-
-@pytest.mark.parametrize(steps_arg_names, steps_arg_values)
-def test_utils_reprojection(
-    radar_array_xr,
-    nwp_array_xr,
-    t,
-):
-
-    # Reproject
-    nwp_array_xr_reproj = reprojection(nwp_array_xr, radar_array_xr)
-
-    # The tests
-    assert (
-        nwp_array_xr_reproj["t"].shape[0] == nwp_array_xr["t"].shape[0]
-    ), "Time dimension has not the same length as source"
-    assert (
-        nwp_array_xr_reproj["y"].shape[0] == radar_array_xr["y"].shape[0]
-    ), "y dimension has not the same length as radar composite"
-    assert (
-        nwp_array_xr_reproj["x"].shape[0] == radar_array_xr["x"].shape[0]
-    ), "x dimension has not the same length as radar composite"
-
-    assert (
-        nwp_array_xr_reproj.x.attrs["x1"] == radar_array_xr.x.attrs["x1"]
-    ), "x-value lower left corner is not equal to radar composite"
-    assert (
-        nwp_array_xr_reproj.x.attrs["x2"] == radar_array_xr.x.attrs["x2"]
-    ), "x-value upper right corner is not equal to radar composite"
-    assert (
-        nwp_array_xr_reproj.y.attrs["y1"] == radar_array_xr.y.attrs["y1"]
-    ), "y-value lower left corner is not equal to radar composite"
-    assert (
-        nwp_array_xr_reproj.y.attrs["y2"] == radar_array_xr.y.attrs["y2"]
-    ), "y-value upper right corner is not equal to radar composite"
-
-    assert (
-        nwp_array_xr_reproj.x.min().values == radar_array_xr.x.min().values
-    ), "First x-coordinate does not equal first x-coordinate radar composite"
-    assert (
-        nwp_array_xr_reproj.y.min().values == radar_array_xr.y.min().values
-    ), "First y-coordinate does not equal first y-coordinate radar composite"
-    assert (
-        nwp_array_xr_reproj.x.max().values == radar_array_xr.x.max().values
-    ), "Last x-coordinate does not equal last x-coordinate radar composite"
-    assert (
-        nwp_array_xr_reproj.y.max().values == radar_array_xr.y.max().values
-    ), "Last y-coordinate does not equal last y-coordinate radar composite"
-
-    assert (
-        nwp_array_xr_reproj.attrs["projection"] == radar_array_xr.attrs["projection"]
-    ), "projection is different than destionation projection"
diff --git a/pysteps/utils/__init__.py b/pysteps/utils/__init__.py
index 9594a75ae..214d9a309 100644
--- a/pysteps/utils/__init__.py
+++ b/pysteps/utils/__init__.py
@@ -11,4 +11,3 @@
 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 eebf84b1d..ee62dd48e 100644
--- a/pysteps/utils/interface.py
+++ b/pysteps/utils/interface.py
@@ -18,7 +18,6 @@
 from . import fft
 from . import images
 from . import interpolate
-from . import reprojection
 from . import spectral
 from . import tapering
 from . import transformation
@@ -188,9 +187,6 @@ def donothing(R, metadata=None, *args, **kwargs):
     methods_objects["rbfinterp2d"] = interpolate.rbfinterp2d
     methods_objects["idwinterp2d"] = interpolate.idwinterp2d
 
-    # reprojection methods
-    methods_objects["reprojection"] = reprojection.reprojection
-
     # 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
deleted file mode 100644
index a702aef29..000000000
--- a/pysteps/utils/reprojection.py
+++ /dev/null
@@ -1,113 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-pysteps.utils.reprojection
-==================
-
-Reprojection tool to reproject the grid and adjust the grid cell size of an 
-input field to a destination field.
-
-.. autosummary::
-    :toctree: ../generated/
-
-    reprojection
-"""
-from pysteps.exceptions import MissingOptionalDependency
-import xarray as xr
-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 reprojection(src_array, dst_array):
-    """Reprojects precipitation fields to the domain of another precipitation
-    field.
-
-    Parameters
-    ----------
-    src_array: xr.DataArray
-        Three-dimensional xarray DataArray with dimensions (t, x, y) containing
-        a time series of precipitation fields. These precipitation fields
-        will be reprojected.
-    dst_array: xr.DataArray
-        Xarray DataArray containing a precipitation field or a time series of
-        precipitation fields. The xarray src_array will be reprojected to the
-        domain of dst_array.
-
-    Returns
-    -------
-    r_rprj: xr.DataArray
-        Three-dimensional xarray DataArray with dimensions (t, x, y) containing
-        the precipitation fields of src_array, but reprojected to the domain of
-        dst_array.
-    """
-
-    if not RASTERIO_IMPORTED:
-        raise MissingOptionalDependency(
-            "rasterio package is required for the reprojection tool, but it is "
-            "not installed"
-        )
-
-    # Extract the grid info from src_array
-    src_crs = src_array.attrs["projection"]
-    x1_src = src_array.x.attrs["x1"]
-    y2_src = src_array.y.attrs["y2"]
-    xpixelsize_src = src_array.attrs["xpixelsize"]
-    ypixelsize_src = src_array.attrs["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 = dst_array.attrs["projection"]
-    x1_dst = dst_array.x.attrs["x1"]
-    y2_dst = dst_array.y.attrs["y2"]
-    xpixelsize_dst = dst_array.attrs["xpixelsize"]
-    ypixelsize_dst = dst_array.attrs["ypixelsize"]
-    dst_transform = A.translation(float(x1_dst), float(y2_dst)) * A.scale(
-        float(xpixelsize_dst), float(-ypixelsize_dst)
-    )
-
-    # Initialise the reprojected (x)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 src_array.attrs["yorigin"] != dst_array.attrs["yorigin"]:
-        src_array = src_array[:, ::-1, :]
-
-    for i in range(src_array["t"].shape[0]):
-        reproject(
-            src_array.isel(t=i).values,
-            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,
-        )
-
-    # Assign the necessary attributes from src_array and dst_array to R_rprj
-    r_rprj = xr.DataArray(
-        data=r_rprj,
-        dims=("t", "y", "x"),
-        coords=dict(
-            t=("t", src_array.coords["t"].data),
-            x=("x", dst_array.coords["x"].data),
-            y=("y", dst_array.coords["y"].data),
-        ),
-    )
-
-    r_rprj.attrs.update(src_array.attrs)
-    r_rprj.x.attrs.update(dst_array.x.attrs)
-    r_rprj.y.attrs.update(dst_array.y.attrs)
-    for key in ["projection", "yorigin", "xpixelsize", "ypixelsize"]:
-        r_rprj.attrs[key] = dst_array.attrs[key]
-
-    return r_rprj
diff --git a/requirements_dev.txt b/requirements_dev.txt
index e1e4a7d44..89ea53ca7 100644
--- a/requirements_dev.txt
+++ b/requirements_dev.txt
@@ -19,7 +19,6 @@ cartopy>=0.18
 h5py
 scikit-image
 pandas
-rasterio
 
 # Testing
 pytest
diff --git a/tox.ini b/tox.ini
index 880d41939..207f92e46 100644
--- a/tox.ini
+++ b/tox.ini
@@ -36,7 +36,6 @@ conda_deps =
     pyproj
     cartopy
     pygrib
-	rasterio
 conda_channels = conda-forge
 
 setenv =