pysteps.blending.steps.forecast(precip, precip_models, velocity, velocity_models, timesteps, timestep, issuetime, n_ens_members, n_cascade_levels=8, blend_nwp_members=False, precip_thr=None, kmperpixel=None, extrap_method='semilagrangian', decomp_method='fft', bandpass_filter_method='gaussian', noise_method='nonparametric', noise_stddev_adj=None, ar_order=2, vel_pert_method='bps', weights_method='bps', conditional=False, probmatching_method='cdf', mask_method='incremental', callback=None, return_output=True, seed=None, num_workers=1, fft_method='numpy', domain='spatial', outdir_path_skill='./tmp/', extrap_kwargs=None, filter_kwargs=None, noise_kwargs=None, vel_pert_kwargs=None, clim_kwargs=None, mask_kwargs=None, measure_time=False)#

Generate a blended nowcast ensemble by using the Short-Term Ensemble Prediction System (STEPS) method.

  • precip (array-like) – Array of shape (ar_order+1,m,n) containing the input precipitation fields ordered by timestamp from oldest to newest. The time steps between the inputs are assumed to be regular.

  • precip_models (array-like) – Array of shape (n_models,timesteps+1) containing, per timestep (t=0 to lead time here) and per (NWP) model or model ensemble member, a dictionary with a list of cascades obtained by calling a method implemented in pysteps.cascade.decomposition. In case of one (deterministic) model as input, add an extra dimension to make sure precip_models is five dimensional prior to calling this function. It is also possible to supply the original (NWP) model forecasts containing rainfall fields as an array of shape (n_models,timestep+1,m,n), which will then be decomposed in this function. Note that for an operational application or for testing with multiple model runs, it is recommended to decompose the model forecasts outside beforehand, as this reduces calculation times. This is possible with pysteps.blending.utils.decompose_NWP(), pysteps.blending.utils.compute_store_nwp_motion(), and pysteps.blending.utils.load_NWP().

  • velocity (array-like) – Array of shape (2,m,n) containing the x- and y-components of the advection field. The velocities are assumed to represent one time step between the inputs. All values are required to be finite.

  • velocity_models (array-like) – Array of shape (n_models,timestep,2,m,n) containing the x- and y-components of the advection field for the (NWP) model field per forecast lead time. All values are required to be finite.

  • timesteps (int or list of floats) – Number of time steps to forecast or a list of time steps for which the forecasts are computed (relative to the input time step). The elements of the list are required to be in ascending order.

  • timestep (float) – Time step of the motion vectors (minutes). Required if vel_pert_method is not None or mask_method is ‘incremental’.

  • issuetime (datetime) – Datetime object containing the date and time for which the forecast is issued.

  • n_ens_members (int) – The number of ensemble members to generate. This number should always be equal to or larger than the number of NWP ensemble members / number of NWP models.

  • n_cascade_levels (int, optional) – The number of cascade levels to use. Default set to 8 due to default climatological skill values on 8 levels.

  • blend_nwp_members (bool) – Check if NWP models/members should be used individually, or if all of them are blended together per nowcast ensemble member. Standard set to false.

  • precip_thr (float, optional) – Specifies the threshold value for minimum observable precipitation intensity. Required if mask_method is not None or conditional is True.

  • kmperpixel (float, optional) – Spatial resolution of the input data (kilometers/pixel). Required if vel_pert_method is not None or mask_method is ‘incremental’.

  • extrap_method (str, optional) – Name of the extrapolation method to use. See the documentation of pysteps.extrapolation.interface.

  • decomp_method ({'fft'}, optional) – Name of the cascade decomposition method to use. See the documentation of pysteps.cascade.interface.

  • bandpass_filter_method ({'gaussian', 'uniform'}, optional) – Name of the bandpass filter method to use with the cascade decomposition. See the documentation of pysteps.cascade.interface.

  • noise_method ({'parametric','nonparametric','ssft','nested',None}, optional) – Name of the noise generator to use for perturbating the precipitation field. See the documentation of pysteps.noise.interface. If set to None, no noise is generated.

  • noise_stddev_adj ({'auto','fixed',None}, optional) – Optional adjustment for the standard deviations of the noise fields added to each cascade level. This is done to compensate incorrect std. dev. estimates of casace levels due to presence of no-rain areas. ‘auto’=use the method implemented in pysteps.noise.utils.compute_noise_stddev_adjs(). ‘fixed’= use the formula given in [BPS06] (eq. 6), None=disable noise std. dev adjustment.

  • ar_order (int, optional) – The order of the autoregressive model to use. Must be >= 1.

  • vel_pert_method ({'bps',None}, optional) – Name of the noise generator to use for perturbing the advection field. See the documentation of pysteps.noise.interface. If set to None, the advection field is not perturbed.

  • weights_method ({'bps','spn'}, optional) – The calculation method of the blending weights. Options are the method by [BPS06] and the covariance-based method by [SPN13]. Defaults to bps.

  • conditional (bool, optional) – If set to True, compute the statistics of the precipitation field conditionally by excluding pixels where the values are below the threshold precip_thr.

  • mask_method ({'obs','incremental',None}, optional) – The method to use for masking no precipitation areas in the forecast field. The masked pixels are set to the minimum value of the observations. ‘obs’ = apply precip_thr to the most recently observed precipitation intensity field, ‘incremental’ = iteratively buffer the mask with a certain rate (currently it is 1 km/min), None=no masking.

  • probmatching_method ({'cdf','mean',None}, optional) – Method for matching the statistics of the forecast field with those of the most recently observed one. ‘cdf’=map the forecast CDF to the observed one, ‘mean’=adjust only the conditional mean value of the forecast field in precipitation areas, None=no matching applied. Using ‘mean’ requires that mask_method is not None.

  • callback (function, optional) – Optional function that is called after computation of each time step of the nowcast. The function takes one argument: a three-dimensional array of shape (n_ens_members,h,w), where h and w are the height and width of the input field precip, respectively. This can be used, for instance, writing the outputs into files.

  • return_output (bool, optional) – Set to False to disable returning the outputs as numpy arrays. This can save memory if the intermediate results are written to output files using the callback function.

  • seed (int, optional) – Optional seed number for the random generators.

  • num_workers (int, optional) – The number of workers to use for parallel computation. Applicable if dask is enabled or pyFFTW is used for computing the FFT. When num_workers>1, it is advisable to disable OpenMP by setting the environment variable OMP_NUM_THREADS to 1. This avoids slowdown caused by too many simultaneous threads.

  • fft_method (str, optional) – A string defining the FFT method to use (see FFT methods in pysteps.utils.interface.get_method()). Defaults to ‘numpy’ for compatibility reasons. If pyFFTW is installed, the recommended method is ‘pyfftw’.

  • domain ({"spatial", "spectral"}) – If “spatial”, all computations are done in the spatial domain (the classical STEPS model). If “spectral”, the AR(2) models and stochastic perturbations are applied directly in the spectral domain to reduce memory footprint and improve performance [PCH19].

  • outdir_path_skill (string, optional) – Path to folder where the historical skill are stored. Defaults to path_workdir from rcparams. If no path is given, ‘./tmp’ will be used.

  • extrap_kwargs (dict, optional) – Optional dictionary containing keyword arguments for the extrapolation method. See the documentation of pysteps.extrapolation.interface().

  • filter_kwargs (dict, optional) – Optional dictionary containing keyword arguments for the filter method. See the documentation of pysteps.cascade.bandpass_filters.

  • noise_kwargs (dict, optional) – Optional dictionary containing keyword arguments for the initializer of the noise generator. See the documentation of pysteps.noise.fftgenerators.

  • vel_pert_kwargs (dict, optional) –

    Optional dictionary containing keyword arguments ‘p_par’ and ‘p_perp’ for the initializer of the velocity perturbator. The choice of the optimal parameters depends on the domain and the used optical flow method.

    Default parameters from [BPS06]: p_par = [10.88, 0.23, -7.68] p_perp = [5.76, 0.31, -2.72]

    Parameters fitted to the data (optical flow/domain):

    darts/fmi: p_par = [13.71259667, 0.15658963, -16.24368207] p_perp = [8.26550355, 0.17820458, -9.54107834]

    darts/mch: p_par = [24.27562298, 0.11297186, -27.30087471] p_perp = [-7.80797846e+01, -3.38641048e-02, 7.56715304e+01]

    darts/fmi+mch: p_par = [16.55447057, 0.14160448, -19.24613059] p_perp = [14.75343395, 0.11785398, -16.26151612]

    lucaskanade/fmi: p_par = [2.20837526, 0.33887032, -2.48995355] p_perp = [2.21722634, 0.32359621, -2.57402761]

    lucaskanade/mch: p_par = [2.56338484, 0.3330941, -2.99714349] p_perp = [1.31204508, 0.3578426, -1.02499891]

    lucaskanade/fmi+mch: p_par = [2.31970635, 0.33734287, -2.64972861] p_perp = [1.90769947, 0.33446594, -2.06603662]

    vet/fmi: p_par = [0.25337388, 0.67542291, 11.04895538] p_perp = [0.02432118, 0.99613295, 7.40146505]

    vet/mch: p_par = [0.5075159, 0.53895212, 7.90331791] p_perp = [0.68025501, 0.41761289, 4.73793581]

    vet/fmi+mch: p_par = [0.29495222, 0.62429207, 8.6804131 ] p_perp = [0.23127377, 0.59010281, 5.98180004]

    fmi=Finland, mch=Switzerland, fmi+mch=both pooled into the same data set

    The above parameters have been fitten by using and located in the scripts directory.

    See pysteps.noise.motion for additional documentation.

  • clim_kwargs (dict, optional) – Optional dictionary containing keyword arguments for the climatological skill file. Arguments can consist of: ‘outdir_path’, ‘n_models’ (the number of NWP models) and ‘window_length’ (the minimum number of days the clim file should have, otherwise the default is used).

  • mask_kwargs (dict) – Optional dictionary containing mask keyword arguments ‘mask_f’ and ‘mask_rim’, the factor defining the the mask increment and the rim size, respectively. The mask increment is defined as mask_f*timestep/kmperpixel.

  • measure_time (bool) – If set to True, measure, print and return the computation time.


out – If return_output is True, a four-dimensional array of shape (n_ens_members,num_timesteps,m,n) containing a time series of forecast precipitation fields for each ensemble member. Otherwise, a None value is returned. The time series starts from t0+timestep, where timestep is taken from the input precipitation fields precip. If measure_time is True, the return value is a three-element tuple containing the nowcast array, the initialization time of the nowcast generator and the time used in the main loop (seconds).

Return type



[See03], [BPS04], [BPS06], [SPN13], [PCH19]


1. The blending currently does not blend the beta-parameters in the parametric noise method. It is recommended to use the non-parameteric noise method.

2. If blend_nwp_members is True, the BPS2006 method for the weights is suboptimal. It is recommended to use the SPN2013 method instead.

3. Not yet implemented (and neither in the steps nowcasting module): The regression of the lag-1 and lag-2 parameters to their climatological values. See also eq. 12 - 19 in :cite: BPS2004. By doing so, the Phi parameters change over time, which enhances the AR process. This can become a future development if this turns out to be a warranted functionality.