Add importer for DWD radar data products and reprojection onto a regular grid as part of #464

pysteps/io/interface.py

@@ -29,6 +29,8 @@
     odim_hdf5=importers.import_odim_hdf5,
     opera_hdf5=importers.import_opera_hdf5,
     knmi_hdf5=importers.import_knmi_hdf5,
+    dwd_hdf5=importers.import_dwd_hdf5,
+    dwd_radolan=importers.import_dwd_radolan,
     saf_crri=importers.import_saf_crri,
 )
 

pysteps/pystepsrc

@@ -110,6 +110,17 @@
                 "gzipped": true
             }
         },
+        "dwd": {
+            "root_path": "./radar/dwd",
+            "path_fmt": "%Y/%m/%d",
+            "fn_pattern": "%Y%m%d_%H%M_RY",
+            "fn_ext": "h5",
+            "importer": "dwd_hdf5",
+            "timestep": 5,
+            "importer_kwargs": {
+                "qty": "RATE"
+            }
+        },
         "bom_nwp": {
             "root_path": "./nwp/bom",
             "path_fmt": "%Y/%m/%d",
@@ -142,6 +153,19 @@
             "importer_kwargs": {
                 "gzipped": true
             }
+        },
+        "dwd_nwp": {
+            "root_path": "./nwp/dwd",
+            "path_fmt": "%Y/%m/%d",
+            "fn_ext": ".grib2",
+            "fn_pattern": "%Y%m%d_%H%M_PR_GSP_060_120",
+            "importer": "dwd_nwp",
+            "timestep": 5,
+            "importer_kwargs": {
+                "varname": "PR_GSP",
+                "grid_file_path": "./aux/grid_files/dwd/icon/R19B07/icon_grid_0047_R19B07_L.nc"
+            }
+
         }
     }
 }

pysteps/tests/test_io_dwd_hdf5.py

@@ -0,0 +1,61 @@
+# -*- coding: utf-8 -*-
+
+import os
+
+import pytest
+
+import pysteps
+from pysteps.tests.helpers import smart_assert
+
+pytest.importorskip("h5py")
+
+# Test for RADOLAN RY product
+
+root_path = pysteps.rcparams.data_sources["dwd"]["root_path"]
+
+filename = os.path.join(root_path, "RY", "2025", "06", "04", "20250604_1700_RY.h5")
+precip_ry, _, metadata_ry = pysteps.io.import_dwd_hdf5(filename, qty="RATE")
+
+
+def test_io_import_dwd_hdf5_ry_shape():
+    """Test the importer DWD HDF5."""
+    assert precip_ry.shape == (1200, 1100)
+
+
+# Test_metadata
+# Expected projection definition
+expected_proj = (
+    "+proj=stere +lat_0=90 +lat_ts=60 "
+    "'+lon_0=10 +a=6378137 +b=6356752.3142451802"
+    "+no_defs +x_0=543196.83521776402 "
+    "+y_0=3622588.8619310018 +units=m"
+)
+
+# 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-10),
+    ("x2", 1099500.0, 1e-10),
+    ("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-10),
+    ("threshold", 0.12, 1e-10),
+]
+
+
+@pytest.mark.parametrize("variable, expected, tolerance", test_ry_attrs)
+def test_io_import_opera_hdf5_odyssey_dataset_attrs(variable, expected, tolerance):
+    """Test the importer OPERA HDF5."""
+    smart_assert(metadata_ry[variable], expected, tolerance)

pysteps/tests/test_utils_reprojection.py

@@ -5,8 +5,11 @@
 import pytest
 import pysteps
 from pysteps.utils import reprojection as rpj
+import pysteps_nwp_importers
 
 pytest.importorskip("rasterio")
+pytest.importorskip("pyproj")
+# pytest.importorskip("pysteps_nwp_importers")
 
 root_path_radar = pysteps.rcparams.data_sources["rmi"]["root_path"]
 
@@ -109,4 +112,67 @@ def test_utils_reproject_grids(
 
     assert (
         metadata_reproj["projection"] == metadata_dst["projection"]
-    ), "projection is different than destionation projection"
+    ), "projection is different than destination projection"
+
+
+root_path_radar = pysteps.rcparams.data_sources["dwd"]["root_path"]
+root_path_nwp = pysteps.rcparams.data_sources["dwd_nwp"]["root_path"]
+
+filename_radar = os.path.join(
+    root_path_radar, "RY", "2025", "06", "04", "20250604_1700_RY.h5"
+)
+_, _, metadata_dst = pysteps.io.import_dwd_hdf5(filename_radar, qty="RATE")
+
+filename_nwp = os.path.join(
+    root_path_nwp, "2025", "06", "04", "20250604_1600_PR_GSP_060_120.grib2"
+)
+kwargs = {
+    "varname": "PR_GSP",
+    "grid_file_path": "./aux/grid_files/dwd/icon/R19B07/icon_grid_0047_R19B07_L.nc",
+}
+array_src, _, metadata_src = pysteps_nwp_importers.importer_dwd_nwp.import_dwd_nwp(
+    filename_nwp, **kwargs
+)
+
+# Since output of NWP importer is based on xarray and the restructure function is
+# rewritten for np.ndarray, there's is some improvisation
+metadata_src["clon"] = array_src["longitude"].values
+metadata_src["clat"] = array_src["latitude"].values
+array_src = array_src.values
+
+restructure_arg_names = ("array_src", "metadata_src", "metadata_dst")
+restructure_arg_values = [(array_src, metadata_src, metadata_dst)]
+
+
+@pytest.mark.parametrize(restructure_arg_names, restructure_arg_values)
+def test_utils_unstructured2regular(array_src, metadata_src, metadata_dst):
+    # Run unstructured2regular
+    array_rprj, metadata_rprj = rpj.unstructured2regular(
+        array_src, metadata_src, metadata_dst
+    )
+
+    # The tests
+    assert (
+        array_rprj.shape[0] == array_src.shape[0]
+    ), "Time dimension has not the same length as source"
+    assert (
+        array_rprj.shape[1] == array_src.shape[0],
+        "Ensemble member dimension has not the same length as source",
+    )
+
+    assert (
+        metadata_rprj["x1"] == metadata_dst["x1"]
+    ), "x-value lower left corner is not equal to radar composite"
+    assert (
+        metadata_rprj["x2"] == metadata_dst["x2"]
+    ), "x-value upper right corner is not equal to radar composite"
+    assert (
+        metadata_rprj["y1"] == metadata_dst["y1"]
+    ), "y-value lower left corner is not equal to radar composite"
+    assert (
+        metadata_rprj["y2"] == metadata_dst["y2"]
+    ), "y-value upper right corner is not equal to radar composite"
+
+    assert (
+        metadata_rprj["projection"] == metadata_dst["projection"]
+    ), "projection is different than destination projection"

pysteps/utils/reprojection.py

@@ -15,7 +15,6 @@
 from scipy.interpolate import griddata
 
 import numpy as np
-import xarray as xr
 
 try:
     from rasterio import Affine as A