create regression test for example scripts

.github/workflows/test_pysteps.yml

@@ -100,7 +100,17 @@ jobs:
         working-directory: ${{github.workspace}}/pysteps_data
         env:
           PYSTEPSRC: ${{github.workspace}}/pysteps_data/pystepsrc
-        run: pytest --pyargs pysteps --cov=pysteps --cov-report=xml --cov-report=term -ra
+          PYSTEPS_REPO_ROOT: ${{github.workspace}}
+        run: pytest --pyargs pysteps --cov=pysteps --cov-report=xml --cov-report=term -ra --basetemp=${{github.workspace}}/pytest-tmp
+
+      - name: Upload example test images
+        if: always()
+        uses: actions/upload-artifact@v4
+        with:
+          name: example-test-images-${{ matrix.os }}-${{ matrix.python-version }}
+          path: ${{github.workspace}}/pytest-tmp/**/fig_*.png
+          if-no-files-found: ignore
+          retention-days: 7
 
       - name: Upload coverage to Codecov (Linux only)
         if: matrix.os == 'ubuntu-latest'

ci/ci_test_env.yml

@@ -34,3 +34,4 @@ dependencies:
   - pytest-cov
   - pip:
       - cookiecutter
+      - git+https://github.com/pySTEPS/pysteps-nwp-importers

examples/ens_kalman_filter_blended_forecast.py

@@ -22,13 +22,13 @@
 from datetime import datetime, timedelta
 
 import numpy as np
+import pysteps_nwp_importers
 from matplotlib import pyplot as plt
 
 import pysteps
-from pysteps import io, rcparams, blending
+from pysteps import blending, io, rcparams
 from pysteps.utils import aggregate_fields_space
 from pysteps.visualization import plot_precip_field
-import pysteps_nwp_importers
 
 ################################################################################
 # Read the radar images and the NWP forecast
@@ -87,11 +87,9 @@
     + nwp_data_source["fn_ext"],
 )
 nwp_importer = io.get_method("dwd_nwp", "importer")
-kwargs = nwp_data_source["importer_kwargs"]
-# Resolve grid_file_path relative to PYSTEPS_DATA_PATH
-kwargs["grid_file_path"] = os.path.join(
-    os.environ["PYSTEPS_DATA_PATH"], kwargs["grid_file_path"]
-)
+# grid_file_path in importer_kwargs is already resolvable from the repo root,
+# like every other data_sources path in pystepsrc.
+kwargs = dict(nwp_data_source["importer_kwargs"])
 nwp_precip, _, nwp_metadata = nwp_importer(filename, **kwargs)
 # We lower the number of ens members to 10 to reduce the memory needs in the
 # example here. However, it is advised to have a minimum of 20 members for the
@@ -250,14 +248,18 @@
     n_cascade_levels=6,  # No. of cascade levels
     precip_thr=log_thr_prec,  # Precip threshold
     norain_thr=0.0005,  # Minimum of 0.5% precip needed, otherwise 'zero rainfall'
-    num_workers=4,  # No. of parallel threads
+    # No. of parallel threads. Kept at 1 so that the forecast (and the plots
+    # in this example) are exactly reproducible: multi-threaded FFT/dask
+    # reductions are not guaranteed to be bit-for-bit deterministic.
+    num_workers=1,
     noise_stddev_adj="auto",  # Standard deviation adjustment
     noise_method="ssft",  # SSFT as noise method
     enable_combination=True,  # Enable combination
     noise_kwargs={"win_size": (512, 512), "win_fun": "hann", "overlap": 0.5},
     extrap_kwargs={"interp_order": 3, "map_coordinates_mode": "nearest"},
     combination_kwargs=combination_kwargs,
     filter_kwargs={"include_mean": True},
+    seed=42,  # Fixed seed for reproducible ensemble members
 )
 
 # Transform the data back into mm/h
@@ -304,6 +306,7 @@
         colorbar=False,
     )
 
+plt.show()
 
 ################################################################################
 # References

examples/linda_nowcasts.py

@@ -137,6 +137,7 @@
     n_ens_members=40,
     num_workers=8,
     measure_time=True,
+    seed=42,  # Fixed seed for reproducible ensemble members
 )[0]
 
 # Compute 40-member STEPS nowcast for comparison
@@ -151,6 +152,7 @@
     kmperpixel=2.0,
     timestep=datasource_params["timestep"],
     vel_pert_method=None,
+    seed=42,  # Fixed seed for reproducible ensemble members
 )
 
 # Convert reflectivity nowcast to rain rate

examples/optical_flow_methods_convergence.py

@@ -18,16 +18,16 @@
 Let's first load the libraries that we will use.
 """
 
-from datetime import datetime
 import time
+from datetime import datetime
 
 import matplotlib.pyplot as plt
 import numpy as np
 from matplotlib.pyplot import get_cmap
 from scipy.ndimage import uniform_filter
 
 import pysteps as stp
-from pysteps import motion, io, rcparams
+from pysteps import io, motion, rcparams
 from pysteps.motion.vet import morph
 from pysteps.visualization import plot_precip_field, quiver
 
@@ -237,12 +237,14 @@ def create_observations(input_precip, motion_type, num_times=9):
 
     synthetic_observations = np.ma.masked_invalid(synthetic_observations)
 
-    synthetic_observations.data[np.ma.getmaskarray(synthetic_observations)] = 0
+    synthetic_observations.data[np.ma.getmaskarray(synthetic_observations)] = np.nan
 
     return ideal_motion, synthetic_observations
 
 
-def plot_optflow_method_convergence(input_precip, optflow_method_name, motion_type):
+def plot_optflow_method_convergence(
+    input_precip, optflow_method_name, motion_type, remove_nans=False
+):
     """
     Test the convergence to the actual solution of the optical flow method used.
 
@@ -276,7 +278,11 @@ def plot_optflow_method_convergence(input_precip, optflow_method_name, motion_ty
 
     elapsed_time = time.perf_counter()
 
-    computed_motion = oflow_method(precip_obs, verbose=False)
+    precip_for_motion = precip_obs.data
+    if remove_nans:
+        precip_for_motion[precip_obs.mask] = np.nanmin(precip_for_motion)
+
+    computed_motion = oflow_method(precip_for_motion, verbose=False)
 
     print(
         f"{optflow_method_name} computation time: "
@@ -299,11 +305,13 @@ def plot_optflow_method_convergence(input_precip, optflow_method_name, motion_ty
     # Compare retrieved motion field with the ideal one
     plt.figure(figsize=(9, 4))
     plt.subplot(1, 2, 1)
-    ax = plot_precip_field(precip_obs[0], title="Reference motion")
+    precip_plot_data = precip_obs[0].data
+    precip_plot_data[precip_obs[0].mask] = np.nan
+    ax = plot_precip_field(precip_plot_data, title="Reference motion")
     quiver(ideal_motion, step=25, ax=ax)
 
     plt.subplot(1, 2, 2)
-    ax = plot_precip_field(precip_obs[0], title="Retrieved motion")
+    ax = plot_precip_field(precip_plot_data, title="Retrieved motion")
     quiver(computed_motion, step=25, ax=ax)
 
     # To evaluate the accuracy of the computed_motion vectors, we will use
@@ -314,9 +322,7 @@ def plot_optflow_method_convergence(input_precip, optflow_method_name, motion_ty
     mse = ((ideal_motion - computed_motion)[:, precip_mask] ** 2).mean()
 
     rel_mse = mse / (ideal_motion[:, precip_mask] ** 2).mean()
-    plt.suptitle(
-        f"{optflow_method_name} " f"Relative RMSE: {np.sqrt(rel_mse) * 100:.2f}%"
-    )
+    plt.suptitle(f"{optflow_method_name} Relative RMSE: {np.sqrt(rel_mse) * 100:.2f}%")
     plt.show()
 
 
@@ -362,17 +368,17 @@ def plot_optflow_method_convergence(input_precip, optflow_method_name, motion_ty
 #
 # Constant motion x-direction
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-plot_optflow_method_convergence(reference_field, "DARTS", "linear_x")
+plot_optflow_method_convergence(reference_field, "DARTS", "linear_x", remove_nans=True)
 
 ################################################################################
 # Constant motion y-direction
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-plot_optflow_method_convergence(reference_field, "DARTS", "linear_y")
+plot_optflow_method_convergence(reference_field, "DARTS", "linear_y", remove_nans=True)
 
 ################################################################################
 # Rotational motion
 # ~~~~~~~~~~~~~~~~~
-plot_optflow_method_convergence(reference_field, "DARTS", "rotor")
+plot_optflow_method_convergence(reference_field, "DARTS", "rotor", remove_nans=True)
 
 ################################################################################
 # Farneback

examples/rainfarm_downscale.py

@@ -114,6 +114,9 @@ def upscale_field(precip, metadata, scale_factor):
 #
 # We can now use RainFARM to downscale the precipitation field.
 
+# Fixed seed for reproducible stochastic noise fields
+np.random.seed(42)
+
 # Basic downscaling
 precip_hr = rainfarm.downscale(precip_lr, ds_factor=scale_factor)
 

examples/steps_blended_forecast.py

@@ -16,7 +16,7 @@
 from matplotlib import pyplot as plt
 
 import pysteps
-from pysteps import io, rcparams, blending, nowcasts
+from pysteps import blending, io, nowcasts, rcparams
 from pysteps.visualization import plot_precip_field
 
 ################################################################################
@@ -169,6 +169,7 @@
     kmperpixel=radar_metadata["xpixelsize"] / 1000.0,
     noise_stddev_adj="auto",
     vel_pert_method=None,
+    seed=42,  # Fixed seed for reproducible ensemble members
 )
 
 # Transform the data back into mm/h
@@ -349,6 +350,7 @@
     )
     ax3.axis("off")
 
+plt.show()
 
 ################################################################################
 # Blend the external nowcast with NWP - ensemble mode
@@ -375,6 +377,7 @@
     probmatching_method="cdf",
     mask_method="incremental",
     weights_method="bps",
+    seed=42,  # Fixed seed for reproducible ensemble members
 )
 
 # Transform the data back into mm/h

pysteps/blending/ens_kalman_filter_methods.py

@@ -433,6 +433,16 @@ def __init__(self, config, params):
         self.__degradation_timestep = 0.2
         self.__inflation_factor_obs_tmp = 1.0
 
+        # Dedicated, per-member random generators for the iterative probability
+        # matching resampling below, so that "seed" fully controls the
+        # reproducibility of the forecast instead of relying on the global,
+        # unseeded numpy random state.
+        seed = self._config.seed
+        self.__randgen_probmatch = [
+            np.random.RandomState(None if seed is None else seed + j)
+            for j in range(self._config.n_ens_members)
+        ]
+
         print("Initialize masked ensemble Kalman filter")
         print("========================================")
         print("")
@@ -603,6 +613,7 @@ def worker(j):
                     first_array=background_ensemble[j],
                     second_array=observation_ensemble[j],
                     probability_first_array=1 - self.__sampling_probability,
+                    randgen=self.__randgen_probmatch[j],
                 ).reshape(self.__params.len_y, self.__params.len_x)
 
             dask_worker_collection = []

pysteps/blending/pca_ens_kalman_filter.py

@@ -666,6 +666,15 @@ def __init__(
 
         self.__ens_member = ens_member
 
+        # Dedicated, per-member random generator for selecting noise fields from
+        # the pool in __iterate, so that the "seed" config option fully controls
+        # the reproducibility of the forecast instead of relying on the global,
+        # unseeded numpy random state.
+        seed = self.__forecast_state.config.seed
+        self.__randgen_select = np.random.RandomState(
+            None if seed is None else seed + ens_member
+        )
+
     # Bundle single steps of the forecast.
     def run_forecast_step(self, nwp, is_correction_timestep=False):
         # Decompose precipitation field.
@@ -813,7 +822,9 @@ def __iterate(self):
         # precipitation mask and the standard deviation coefficients.
         epsilon = (
             self.__forecast_state.noise_field_pool[
-                np.random.randint(self.__forecast_state.config.n_noise_fields)
+                self.__randgen_select.randint(
+                    self.__forecast_state.config.n_noise_fields
+                )
             ]
             * self.__forecast_state.precip_mask[self.__ens_member][None, :, :]
             * self.__forecast_state.params.noise_std_coeffs[:, None, None]

pysteps/io/importers.py

@@ -519,10 +519,10 @@ def _import_bom_rf3_geodata(ds_rainfall):
         if getattr(ds_rainfall.variables["x"], "units") == "km":
             factor_scale = 1000.0
 
-    geodata["x1"] = xmin * factor_scale
-    geodata["y1"] = ymin * factor_scale
-    geodata["x2"] = xmax * factor_scale
-    geodata["y2"] = ymax * factor_scale
+    geodata["x1"] = (xmin - xpixelsize * 0.5) * factor_scale
+    geodata["y1"] = (ymin - ypixelsize * 0.5) * factor_scale
+    geodata["x2"] = (xmax + xpixelsize * 0.5) * factor_scale
+    geodata["y2"] = (ymax + ypixelsize * 0.5) * factor_scale
     geodata["xpixelsize"] = xpixelsize * factor_scale
     geodata["ypixelsize"] = ypixelsize * factor_scale
     geodata["cartesian_unit"] = "m"
@@ -1807,37 +1807,6 @@ def import_dwd_hdf5(filename, qty="RATE", **kwargs):
     if precip is None:
         raise IOError("requested quantity %s not found" % qty)
 
-    # Get the projection and grid information from the HDF5 file
-    pr = pyproj.Proj(file_content["where"]["projdef"])
-    ll_x, ll_y = pr(
-        file_content["where"]["LL_lon"],
-        file_content["where"]["LL_lat"],
-    )
-    ur_x, ur_y = pr(
-        file_content["where"]["UR_lon"],
-        file_content["where"]["UR_lat"],
-    )
-
-    # Determine domain corners in geographic and carthesian coordinates
-    if len([k for k in file_content["where"].keys() if "_lat" in k]) == 4:
-        lr_x, lr_y = pr(
-            file_content["where"]["LR_lon"],
-            file_content["where"]["LR_lat"],
-        )
-        ul_x, ul_y = pr(
-            file_content["where"]["UL_lon"],
-            file_content["where"]["UL_lat"],
-        )
-        x1 = min(ll_x, ul_x)
-        y1 = min(ll_y, lr_y)
-        x2 = max(lr_x, ur_x)
-        y2 = max(ul_y, ur_y)
-    else:
-        x1 = ll_x
-        y1 = ll_y
-        x2 = ur_x
-        y2 = ur_y
-
     # Get the grid cell size
     if (
         "where" in file_content["dataset1"].keys()
@@ -1849,8 +1818,12 @@ def import_dwd_hdf5(filename, qty="RATE", **kwargs):
         xpixelsize = file_content["where"]["xscale"]
         ypixelsize = file_content["where"]["yscale"]
     else:
-        xpixelsize = None
-        ypixelsize = None
+        raise ValueError("No xscale found in ODIM hdf5 file")
+
+    x1 = -xpixelsize * 0.5
+    x2 = -xpixelsize * 0.5 + precip.shape[1] * xpixelsize
+    y1 = xpixelsize * 0.5 - precip.shape[0] * ypixelsize
+    y2 = ypixelsize * 0.5
 
     # Get the unit and transform
     if qty == "ACRR":

pysteps/motion/darts.py

@@ -11,8 +11,9 @@
     DARTS
 """
 
-import numpy as np
 import time
+
+import numpy as np
 from numpy.linalg import lstsq, svd
 
 from pysteps import utils