Notebook

Run it live:

Download this notebook: 1a_basics.ipynb

Running a Hydrometeorological Ensemble Verification Pipeline with veriflow

Forecast verification helps us evaluate the quality and behaviour of forecasts. This can include comparing forecasts with observations, assessing skill across lead times, examining forecast uncertainty, or checking whether forecast distributions are realistic under different hydrological conditions.

veriflow is a Python-based framework for running forecast verification workflows in a consistent and reproducible way. It uses a configuration-driven approach to load data, define verification experiments, compute verification metrics and statistical criteria, and return results. It can be used to verify both deterministic and probabilistic ensemble forecasts.

1. What this notebook does

This notebook focuses on the workflow, with an initial visualization of verification results:

1. Define a verification configuration

2. Run the veriflow pipeline

3. Inspect the returned results

4. Explore the available data and diagnostics

5. Visualize the results with the interactive app

The aim of this notebook is to demonstrate usage of veriflow, not to provide a complete overview of forecast verification theory.

A second notebook focuses on the interpretation of the data, verification results, and forecast quality.

---

2. Dataset and case study

This example uses hydrometeorological data for the River Rhine basin. We focus on discharge and compare three discharge forecast datasets. Each dataset consists of a 5-member ensemble reforecast. The discharge forecasts were generated by forcing the HBV hydrological model with raw or post-processed ECMWF ensemble reforecasts of precipitation and temperature.

The forecast datasets are:

• ``raw_raw``: baseline streamflow forecast driven by raw meteorological ensemble forecasts

• ``lin_log``: streamflow forecast driven by meteorological ensembles post-processed using a linear-log transformation

• ``qqt_qqt``: streamflow forecast driven by meteorological ensembles post-processed using a quantile-to-quantile transformation

In this notebook, all forecast datasets are evaluated against the same observation dataset. We focus on a key location, Lobith.

The data are read directly from a remote Zarr store, so no local download is required. For more information, see the References section.

3. Imports and environment setup

We start by importing the required Python libraries and veriflow components.

The imports include:

• Required Python tools

• veriflow components, including configuration objects, Zarr datasource for reading the example data remotely, and pipeline runner

• A local app for interactive visualization

# Add automatic reloading of modules in case of changes
%load_ext autoreload
%autoreload 2

from datetime import datetime, timezone
import logging
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt
from pathlib import Path
from veriflow.configuration import Config, GeneralInfoConfig
from veriflow.configuration.utils import (
    LeadTimes,
    Range,
    TimeUnits,
    VerificationPair,
    VerificationPeriod,
)
from veriflow.constants import DataType
from veriflow.datasources import ZarrConfig
from veriflow.pipeline import run_pipeline
from veriflow.scores import CrpsForEnsembleConfig, RankHistogramConfig

import sys

sys.path.append(str(Path("..").resolve()))
from verification_plots import (
    crps_plot,
    find_crps_variable,
    rank_histogram_3d_plot,
    rank_histogram_plot,
    scatter_plot,
)

---

4. Defining the verification configuration

veriflow is a configuration-driven workflow. The configuration tells the pipeline:

• which data sources to read

• which forecasts and observations to compare

• which lead times to evaluate

• which verification metrics and statistical criteria to compute

We define these parts below and then run the full workflow with a single pipeline call.

4.1. General verification settings

The general settings define the structure of the verification experiment:

• the verification period

• the forecast lead times

• the verification pairs

Each verification pair represents a single comparison between a forecast dataset and the observations.

general = GeneralInfoConfig(
    # Evaluate forecasts issued during this period.
    verification_period=VerificationPeriod(
        start=datetime(1988, 1, 1, tzinfo=timezone.utc),
        end=datetime(2008, 12, 31, tzinfo=timezone.utc),
        dimension="forecast_reference_time",
    ),
    # Define which forecast datasets are compared against observations.
    verification_pairs=[
        VerificationPair(id="raw_raw", obs="obs", sim="raw_raw", variable="Q"),
        VerificationPair(id="lin_log", obs="obs", sim="lin_log", variable="Q"),
        VerificationPair(id="qqt_qqt", obs="obs", sim="qqt_qqt", variable="Q"),
    ],
    # Evaluate lead times from 0 to 10 days.
    lead_times=LeadTimes(
        unit=TimeUnits.day,
        values=Range(start=0, end=10, step=1),
    ),
)

4.2. Data sources

The example data is stored in remote Zarr datasets. Zarr is well suited for cloud-based, multidimensional data access.

Here we define one observation source and three forecast sources. The source names must match the names used in the verification pairs above.

datasources = [
    ZarrConfig(
        general=general,
        import_adapter="zarr",
        source="obs",
        data_type=DataType.observed_historical,
        consolidated=True,
        path="https://s3.deltares.nl/deltares-verification-assets/rhine_dataset_verkade_2013/optimized/obs_Q.zarr",
    ),
    ZarrConfig(
        general=general,
        import_adapter="zarr",
        source="raw_raw",
        data_type=DataType.simulated_forecast_ensemble,
        consolidated=True,
        path="https://s3.deltares.nl/deltares-verification-assets/rhine_dataset_verkade_2013/optimized/case-raw-raw_Q.zarr",
    ),
    ZarrConfig(
        general=general,
        import_adapter="zarr",
        source="lin_log",
        data_type=DataType.simulated_forecast_ensemble,
        consolidated=True,
        path="https://s3.deltares.nl/deltares-verification-assets/rhine_dataset_verkade_2013/optimized/case-lin-log_Q.zarr",
    ),
    ZarrConfig(
        general=general,
        import_adapter="zarr",
        source="qqt_qqt",
        data_type=DataType.simulated_forecast_ensemble,
        consolidated=True,
        path="https://s3.deltares.nl/deltares-verification-assets/rhine_dataset_verkade_2013/optimized/case-qqt-qqt_Q.zarr",
    ),
]

4.3. Verification metrics and statistical criteria

In this notebook, we compute two commonly used criteria for ensemble forecasts:

Continuous Ranked Probability Score (CRPS) assesses how well an ensemble forecast matches the observed outcome.

• A CRPS of 0 indicates a perfect forecast.

• Lower values indicate better probabilistic forecast quality, indicating the ensemble is centred around the observation and does not unnecessarily spread out.

• CRPS increases when the forecast distribution is biased, too wide, or too narrow relative to the observation.

Rank Histogram assesses statistical reliability by showing where observations fall within the ensemble distribution.

• A roughly uniform histogram indicates a well-calibrated ensemble.

• Systematic patterns can reveal biases or issues with ensemble spread.

The interactive app later in this notebook can be used to explore both reuslting CRPS and rank histogram.

scores = [
    # CRPS
    CrpsForEnsembleConfig(
        score_adapter="crps_for_ensemble",
        general=general,
        reduce_dims=[],
    ),
    # Rank histogram
    RankHistogramConfig(
        score_adapter="rank_histogram",
        general=general,
        reduce_dims=["forecast_reference_time"],
    ),
]

config = Config(
    fileversion="0.1.0",
    general=general,
    datasources=datasources,
    scores=scores,
)

---

5. Running the verification pipeline

With the configuration in place, we can run the veriflow pipeline, which combines the following steps into a workflow:

1. Load the input datasets

2. Pair forecasts and observations

3. Compute the requested verification metrics and statistical criteria

4. Return everything in an OutputDataset

import logging

logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s %(levelname)s %(name)s: %(message)s",
    handlers=[logging.StreamHandler()],
)

output_dataset = run_pipeline(config)

2026-07-24 15:52:25,183 INFO veriflow.pipeline: Successfully initialized the configuration.
         verification_period_start = 1988-01-01 00:00:00
         verification_period_end = 2008-12-31 00:00:00
2026-07-24 15:52:25,184 INFO veriflow.datasources.base: Starting data fetch for obs from Zarr.
2026-07-24 15:52:27,538 INFO veriflow.datasources.base: Starting data fetch for raw_raw from Zarr.
2026-07-24 15:52:33,328 INFO veriflow.datasources.base: Starting data fetch for lin_log from Zarr.
2026-07-24 15:52:36,195 INFO veriflow.datasources.base: Starting data fetch for qqt_qqt from Zarr.
2026-07-24 15:52:39,172 INFO veriflow.pipeline: Successfully loaded all data from sources.
2026-07-24 15:52:40,887 INFO veriflow.pipeline: Successfully computed CrpsForEnsemble for verification pair raw_raw.
2026-07-24 15:52:42,485 INFO veriflow.pipeline: Successfully computed CrpsForEnsemble for verification pair lin_log.
2026-07-24 15:52:44,099 INFO veriflow.pipeline: Successfully computed CrpsForEnsemble for verification pair qqt_qqt.
2026-07-24 15:52:45,073 INFO veriflow.pipeline: Successfully computed RankHistogram for verification pair raw_raw.
2026-07-24 15:52:46,039 INFO veriflow.pipeline: Successfully computed RankHistogram for verification pair lin_log.
2026-07-24 15:52:47,131 INFO veriflow.pipeline: Successfully computed RankHistogram for verification pair qqt_qqt.
2026-07-24 15:52:47,131 INFO veriflow.pipeline: Verification pipeline completed successfully.

---

6. Inspecting the OutputDataset

The pipeline returns an ``OutputDataset``, which provides access to the verification results. Results are organized by verification pair, where each pair defines a comparison between an observation dataset and a forecast dataset (for example, raw_raw versus observations). Each verification pair contains the input data together with the computed verification metrics.

We start by inspecting the output structure and then extract the results for a single verification pair.

type(output_dataset)

veriflow.datamodel.main.OutputDataset

The output contains three verification pairs, corresponding to the three forecast datasets defined in the configuration.

We can now select one verification pair and inspect its dataset. The returned object is an xarray.Dataset containing:

• the forecast data

• the corresponding observations

• the computed verification metrics and statistical criteria (in this case, CRPS and rank histogram)

• metadata and coordinates (e.g., station, lead time, and forecast reference time)

output_dataset.get(output_dataset.verification_pairs[0])

<xarray.Dataset> Size: 83MB
Dimensions:                  (station: 88, forecast_reference_time: 2924,
                              lead_time: 10, realization: 5, rank: 6)
Coordinates:
  * station                  (station) <U11 4kB 'H-RN-0001' ... 'H-RN-WURZ'
    lat                      (station) float64 704B dask.array<chunksize=(12,), meta=np.ndarray>
    lon                      (station) float64 704B dask.array<chunksize=(12,), meta=np.ndarray>
  * forecast_reference_time  (forecast_reference_time) datetime64[ns] 23kB 19...
  * lead_time                (lead_time) timedelta64[ns] 80B 1 days ... 10 days
    time                     (forecast_reference_time, lead_time) datetime64[ns] 234kB dask.array<chunksize=(2924, 10), meta=np.ndarray>
  * realization              (realization) int64 40B 0 1 2 3 4
  * rank                     (rank) float64 48B 1.0 2.0 3.0 4.0 5.0 6.0
Data variables:
    obs                      (station, forecast_reference_time, lead_time) float32 10MB dask.array<chunksize=(88, 730, 10), meta=np.ndarray>
    raw_raw                  (forecast_reference_time, lead_time, station, realization) float32 51MB dask.array<chunksize=(2924, 10, 12, 5), meta=np.ndarray>
    crps_for_ensemble        (forecast_reference_time, lead_time, station) float64 21MB dask.array<chunksize=(730, 10, 12), meta=np.ndarray>
    histogram_rank           (station, lead_time, rank) float64 42kB dask.array<chunksize=(12, 10, 6), meta=np.ndarray>
Attributes:
    data_type:  observed_historical
    units:      m^3/s

xarray.Dataset

• Dimensions:
  • station: 88
  • forecast_reference_time: 2924
  • lead_time: 10
  • realization: 5
  • rank: 6
• Coordinates: (8)
  • station
    (station)
    <U11
    'H-RN-0001' ... 'H-RN-WURZ'

    array(['H-RN-0001', 'H-RN-0024', 'H-RN-0026', 'H-RN-0028', 'H-RN-0029',
           'H-RN-0031', 'H-RN-0036', 'H-RN-0038', 'H-RN-0039', 'H-RN-0052',
           'H-RN-0053', 'H-RN-0627', 'H-RN-0659', 'H-RN-0668', 'H-RN-0689',
           'H-RN-0693', 'H-RN-0808', 'H-RN-0847', 'H-RN-0888', 'H-RN-0900',
           'H-RN-0908', 'H-RN-0913', 'H-RN-0943', 'H-RN-0947', 'H-RN-0950',
           'H-RN-0957', 'H-RN-0984', 'H-RN-1025', 'H-RN-1026', 'H-RN-1027',
           'H-RN-2289', 'H-RN-BFG001', 'H-RN-BFG002', 'H-RN-BFG003', 'H-RN-BFG005',
           'H-RN-BFG007', 'H-RN-BFG008', 'H-RN-BFG009', 'H-RN-BFG014',
           'H-RN-BFG015', 'H-RN-BFG017', 'H-RN-BFG018', 'H-RN-BFG019',
           'H-RN-BFG021', 'H-RN-BFG023', 'H-RN-BFG025', 'H-RN-BFG026',
           'H-RN-BFG028', 'H-RN-BFG029', 'H-RN-BFG031', 'H-RN-BFG033',
           'H-RN-BFG035', 'H-RN-BFG036', 'H-RN-BFG037', 'H-RN-BFG039',
           'H-RN-BFG040', 'H-RN-BFG041', 'H-RN-BFG042', 'H-RN-BFG043',
           'H-RN-BFG044', 'H-RN-BFG045', 'H-RN-BFG046', 'H-RN-BFG048',
           'H-RN-BFG049', 'H-RN-BFG050', 'H-RN-BFG052', 'H-RN-BFG058',
           'H-RN-BFG059', 'H-RN-BFG060', 'H-RN-BFG061', 'H-RN-BFG063',
           'H-RN-BFG066', 'H-RN-BFG068', 'H-RN-BFG069', 'H-RN-BFG071',
           'H-RN-BFG072', 'H-RN-BFG073', 'H-RN-BFG074', 'H-RN-BFG075',
           'H-RN-BFG076', 'H-RN-BFG077', 'H-RN-BFG078', 'H-RN-BFG079', 'H-RN-BOLL',
           'H-RN-DIET', 'H-RN-PLOC', 'H-RN-TRIE', 'H-RN-WURZ'], dtype='<U11')

  • lat
    (station)
    float64
    dask.array<chunksize=(12,), meta=np.ndarray>

    Array Chunk Bytes 704 B 96 B Shape (88,) (12,) Dask graph 8 chunks in 17 graph layers Data type float64 numpy.ndarray

  • lon
    (station)
    float64
    dask.array<chunksize=(12,), meta=np.ndarray>

    Array Chunk Bytes 704 B 96 B Shape (88,) (12,) Dask graph 8 chunks in 17 graph layers Data type float64 numpy.ndarray

  • forecast_reference_time
    (forecast_reference_time)
    datetime64[ns]
    1990-03-13 ... 2007-12-31

    array(['1990-03-13T00:00:00.000000000', '1990-03-20T00:00:00.000000000',
           '1990-03-27T00:00:00.000000000', ..., '2007-12-25T00:00:00.000000000',
           '2007-12-30T00:00:00.000000000', '2007-12-31T00:00:00.000000000'],
          shape=(2924,), dtype='datetime64[ns]')

  • lead_time
    (lead_time)
    timedelta64[ns]
    1 days 2 days ... 9 days 10 days

    array([ 86400000000000, 172800000000000, 259200000000000, 345600000000000,
           432000000000000, 518400000000000, 604800000000000, 691200000000000,
           777600000000000, 864000000000000], dtype='timedelta64[ns]')

  • time
    (forecast_reference_time, lead_time)
    datetime64[ns]
    dask.array<chunksize=(2924, 10), meta=np.ndarray>

    Array Chunk Bytes 228.44 kiB 228.44 kiB Shape (2924, 10) (2924, 10) Dask graph 1 chunks in 8 graph layers Data type datetime64[ns] numpy.ndarray

  • realization
    (realization)
    int64
    0 1 2 3 4

    array([0, 1, 2, 3, 4])

  • rank
    (rank)
    float64
    1.0 2.0 3.0 4.0 5.0 6.0

    data_type :

    observed_historical

    units :

    m^3/s

    array([1., 2., 3., 4., 5., 6.])

• Data variables: (4)
  • obs
    (station, forecast_reference_time, lead_time)
    float32
    dask.array<chunksize=(88, 730, 10), meta=np.ndarray>

    data_type :

    observed_historical

    units :

    m^3/s

    Array Chunk Bytes 9.82 MiB 2.45 MiB Shape (88, 2924, 10) (88, 730, 10) Dask graph 7 chunks in 14 graph layers Data type float32 numpy.ndarray

  • raw_raw
    (forecast_reference_time, lead_time, station, realization)
    float32
    dask.array<chunksize=(2924, 10, 12, 5), meta=np.ndarray>

    data_type :

    simulated_forecast_ensemble

    units :

    m^3/s

    Array Chunk Bytes 49.08 MiB 6.69 MiB Shape (2924, 10, 88, 5) (2924, 10, 12, 5) Dask graph 8 chunks in 27 graph layers Data type float32 numpy.ndarray

  • crps_for_ensemble
    (forecast_reference_time, lead_time, station)
    float64
    dask.array<chunksize=(730, 10, 12), meta=np.ndarray>

    units :

    m^3/s

    data_type :

    simulated_forecast_ensemble

    Array Chunk Bytes 19.63 MiB 684.38 kiB Shape (2924, 10, 88) (730, 10, 12) Dask graph 56 chunks in 71 graph layers Data type float64 numpy.ndarray

  • histogram_rank
    (station, lead_time, rank)
    float64
    dask.array<chunksize=(12, 10, 6), meta=np.ndarray>

    data_type :

    simulated_forecast_ensemble

    units :

    m^3/s

    Array Chunk Bytes 41.25 kiB 5.62 kiB Shape (88, 10, 6) (12, 10, 6) Dask graph 8 chunks in 50 graph layers Data type float64 numpy.ndarray

• Attributes: (2)

  data_type :

  observed_historical

  units :

  m^3/s

---

7. Visualization of verification results

Now that we have an overview of the observed and ensemble forecast time series and the verification results through the ``OutputDataset``, we can visualize the results using a small collection of reusable plotting functions in verification_plots.py.

Each function returns a Plotly figure for a pre-sliced dataset. There is no interactive layer: selections that used to be interactive controls - verification pair, station, and lead time - are now made explicitly by slicing the xarray.Dataset before calling a plot function. This makes the plots easy to reuse in notebooks, scripts, or reports.

The available functions are:

• scatter_plot - observation vs. simulation scatter for a single station and lead time

• crps_plot - CRPS across lead times for one or more pre-selected stations

• rank_histogram_plot - rank histogram for a single station and lead time

• rank_histogram_3d_plot - 3D rank histogram showing all ranks across lead times for a single station

The interpretation of these plots and diagnostics is covered in the second notebook.

first_ds = output_dataset.get(output_dataset.verification_pairs[0])
stations = first_ds.coords["station"].values
lead_times = first_ds.coords["lead_time"].values
print("Stations:", stations)

Stations: ['H-RN-0001' 'H-RN-0024' 'H-RN-0026' 'H-RN-0028' 'H-RN-0029' 'H-RN-0031'
 'H-RN-0036' 'H-RN-0038' 'H-RN-0039' 'H-RN-0052' 'H-RN-0053' 'H-RN-0627'
 'H-RN-0659' 'H-RN-0668' 'H-RN-0689' 'H-RN-0693' 'H-RN-0808' 'H-RN-0847'
 'H-RN-0888' 'H-RN-0900' 'H-RN-0908' 'H-RN-0913' 'H-RN-0943' 'H-RN-0947'
 'H-RN-0950' 'H-RN-0957' 'H-RN-0984' 'H-RN-1025' 'H-RN-1026' 'H-RN-1027'
 'H-RN-2289' 'H-RN-BFG001' 'H-RN-BFG002' 'H-RN-BFG003' 'H-RN-BFG005'
 'H-RN-BFG007' 'H-RN-BFG008' 'H-RN-BFG009' 'H-RN-BFG014' 'H-RN-BFG015'
 'H-RN-BFG017' 'H-RN-BFG018' 'H-RN-BFG019' 'H-RN-BFG021' 'H-RN-BFG023'
 'H-RN-BFG025' 'H-RN-BFG026' 'H-RN-BFG028' 'H-RN-BFG029' 'H-RN-BFG031'
 'H-RN-BFG033' 'H-RN-BFG035' 'H-RN-BFG036' 'H-RN-BFG037' 'H-RN-BFG039'
 'H-RN-BFG040' 'H-RN-BFG041' 'H-RN-BFG042' 'H-RN-BFG043' 'H-RN-BFG044'
 'H-RN-BFG045' 'H-RN-BFG046' 'H-RN-BFG048' 'H-RN-BFG049' 'H-RN-BFG050'
 'H-RN-BFG052' 'H-RN-BFG058' 'H-RN-BFG059' 'H-RN-BFG060' 'H-RN-BFG061'
 'H-RN-BFG063' 'H-RN-BFG066' 'H-RN-BFG068' 'H-RN-BFG069' 'H-RN-BFG071'
 'H-RN-BFG072' 'H-RN-BFG073' 'H-RN-BFG074' 'H-RN-BFG075' 'H-RN-BFG076'
 'H-RN-BFG077' 'H-RN-BFG078' 'H-RN-BFG079' 'H-RN-BOLL' 'H-RN-DIET'
 'H-RN-PLOC' 'H-RN-TRIE' 'H-RN-WURZ']

# Scatter: pre-slice to a single station and a single lead time.
scatter_plot(output_dataset, station=stations[0], lead_time=lead_times[5]).show()

crps_plot(output_dataset, stations=["H-RN-0001"]).show()

# Rank histogram: pre-slice to a single station and a single lead time.
rank_histogram_plot(output_dataset, station="H-RN-0001", lead_time=lead_times[5]).show()

/home/runner/work/veriflow/veriflow/.venv/lib/python3.11/site-packages/xskillscore/core/probabilistic.py:920: RuntimeWarning: invalid value encountered in cast
  ranks = ranks_min + np.int32((ranks_max - ranks_min + 1) * np.random.rand(*xy.shape))

# Advanced and abstract.
rank_histogram_3d_plot(output_dataset, station="H-RN-0001").show()

/home/runner/work/veriflow/veriflow/.venv/lib/python3.11/site-packages/xskillscore/core/probabilistic.py:920: RuntimeWarning: invalid value encountered in cast
  ranks = ranks_min + np.int32((ranks_max - ranks_min + 1) * np.random.rand(*xy.shape))

8. Next step

You have now configured and run a complete veriflow verification workflow, inspected the returned OutputDataset, and opened the interactive visualization app.

The second notebook, Interpreting Hydrological Ensemble Verification Results, builds directly on this output. It focuses on how to inspect the aligned data, interpret scatter plots, CRPS curves, and rank histograms, and draw conclusions about forecast performance, uncertainty, and ensemble reliability.

Optional: YAML-based configuration

This notebook uses Python objects to define the configuration because it is convenient for interactive work.

veriflow can also support YAML-based configuration files, which are useful for reproducible and shareable workflows. A YAML example can be added here once the preferred release configuration format is finalized.

References

The dataset used in the notebook was originally used in: > Verkade et al. (2013), Post-processing ECMWF precipitation and temperature ensemble reforecasts for operational hydrologic forecasting at various spatial scales, Journal of Hydrology. https://doi.org/10.1016/j.jhydrol.2013.07.039

Additional information is available via 4TU.ResearchData:

> https://data.4tu.nl/articles/_/12694658/1