def forecast(
    precip,
    velocity,
    timesteps,
    feature_method="blob",
    max_num_features=25,
    feature_kwargs=None,
    ari_order=1,
    kernel_type="anisotropic",
    localization_window_radius=None,
    errdist_window_radius=None,
    acf_window_radius=None,
    extrap_method="semilagrangian",
    extrap_kwargs=None,
    add_perturbations=True,
    pert_thrs=(0.5, 1.0),
    n_ens_members=10,
    vel_pert_method="bps",
    vel_pert_kwargs=None,
    kmperpixel=None,
    timestep=None,
    seed=None,
    num_workers=1,
    use_multiprocessing=False,
    measure_time=False,
    callback=None,
    return_output=True,
):
    """
    Generate a deterministic or ensemble nowcast by using the Lagrangian
    INtegro-Difference equation model with Autoregression (LINDA) model.

    Parameters
    ----------
    precip: array_like
        Array of shape (ari_order + 2, m, n) containing the input rain rate
        or reflectivity fields (in linear scale) ordered by timestamp from
        oldest to newest. The time steps between the inputs are assumed to be
        regular.
    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.
    timesteps: int
        Number of time steps to forecast.
    feature_method: {'blob', 'domain' 'shitomasi'}
        Feature detection method:

        +-------------------+-----------------------------------------------------+
        |    Method name    |                  Description                        |
        +===================+=====================================================+
        |  blob             | Laplacian of Gaussian (LoG) blob detector           |
        |                   | implemented in scikit-image                         |
        +-------------------+-----------------------------------------------------+
        |  domain           | no feature detection, the model is applied over the |
        |                   | whole domain without localization                   |
        +-------------------+-----------------------------------------------------+
        |  shitomasi        | Shi-Tomasi corner detector implemented in OpenCV    |
        +-------------------+-----------------------------------------------------+

        Default: 'blob'
    max_num_features: int, optional
        Maximum number of features to use. It is recommended to set this between
        20 and 50, which gives a good tradeoff between localization and
        computation time. Default: 25
    feature_kwargs: dict, optional
        Keyword arguments that are passed as ``**kwargs`` for the feature detector.
        See :py:mod:`pysteps.feature.blob` and :py:mod:`pysteps.feature.shitomasi`.
    ari_order: {1, 2}, optional
        The order of the ARI(p, 1) model. Default: 1
    kernel_type: {"anisotropic", "isotropic"}, optional
        The type of the kernel. Default: 'anisotropic'
    localization_window_radius: float, optional
        The standard deviation of the Gaussian localization window.
        Default: 0.2 * min(m, n)
    errdist_window_radius: float, optional
        The standard deviation of the Gaussian window for estimating the
        forecast error distribution. Default: 0.15 * min(m, n)
    acf_window_radius: float, optional
        The standard deviation of the Gaussian window for estimating the
        forecast error ACF. Default: 0.25 * min(m, n)
    extrap_method: str, optional
        The extrapolation method to use. See the documentation of
        :py:mod:`pysteps.extrapolation.interface`. Default: 'semilagrangian'
    extrap_kwargs: dict, optional
        Optional dictionary containing keyword arguments for the extrapolation
        method. See :py:mod:`pysteps.extrapolation.interface`.
    add_perturbations: bool
        Set to False to disable perturbations and generate a single
        deterministic nowcast. Default: True
    pert_thrs: float
        Two-element tuple containing the threshold values for estimating the
        perturbation parameters (mm/h). Default: (0.5, 1.0)
    n_ens_members: int, optional
        The number of ensemble members to generate. Default: 10
    vel_pert_method: {'bps', None}, optional
        Name of the generator to use for perturbing the advection field. See
        :py:mod:`pysteps.noise.interface`. Default: 'bps'
    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. For the default values and parameters optimized for different
        domains, see :py:func:`pysteps.nowcasts.steps.forecast`.
    kmperpixel: float, optional
        Spatial resolution of the input data (kilometers/pixel). Required if
        vel_pert_method is not None.
    timestep: float, optional
        Time step of the motion vectors (minutes). Required if vel_pert_method
        is not None.
    seed: int, optional
        Optional seed for the random generators.
    num_workers: int, optional
        The number of workers to use for parallel computations. Applicable if
        dask is installed. Default: 1
    use_multiprocessing: bool, optional
        Set to True to improve the performance of certain parallelized parts of
        the code. If set to True, the main script calling linda.forecast must
        be enclosed within the 'if __name__ == "__main__":' block.
        Default: False
    measure_time: bool, optional
        If set to True, measure, print and return the computation time.
        Default: False
    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 precipitation fields, respectively. This can be used, for
        instance, writing the outputs into files. Default: None
    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. Default: True

    Returns
    -------
    out: numpy.ndarray
        A four-dimensional array of shape (n_ens_members, timesteps, m, n)
        containing a time series of forecast precipitation fields for each
        ensemble member. If add_perturbations is False, the first dimension is
        dropped. The time series starts from t0 + timestep, where timestep is
        taken from the input fields. 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). If return_output is set to False, a single None value is
        returned instead.

    Notes
    -----
    It is recommended to choose the feature detector parameters so that the
    number of features is around 20-40. This gives a good tradeoff between
    localization and computation time.

    It is highly recommented to set num_workers>1 to reduce computation time.
    In this case, 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.
    """
    _check_inputs(precip, velocity, timesteps, ari_order)

    if feature_kwargs is None:
        feature_kwargs = dict()
    if extrap_kwargs is None:
        extrap_kwargs = dict()
    else:
        extrap_kwargs = extrap_kwargs.copy()

    if localization_window_radius is None:
        localization_window_radius = 0.2 * np.min(precip.shape[1:])
    if add_perturbations:
        if errdist_window_radius is None:
            errdist_window_radius = 0.15 * min(precip.shape[1], precip.shape[2])
        if acf_window_radius is None:
            acf_window_radius = 0.25 * min(precip.shape[1], precip.shape[2])
        if vel_pert_method is not None:
            if kmperpixel is None:
                raise ValueError("vel_pert_method is set but kmperpixel is None")
            if timestep is None:
                raise ValueError("vel_pert_method is set but timestep is None")
        if vel_pert_kwargs is None:
            vel_pert_kwargs = dict()

    print("Computing LINDA nowcast")
    print("-----------------------")
    print("")

    print("Inputs")
    print("------")
    print(f"dimensions:           {precip.shape[1]}x{precip.shape[2]}")
    print(f"number of time steps: {precip.shape[0]}")
    print("")

    print("Methods")
    print("-------")
    nowcast_type = "ensemble" if add_perturbations else "deterministic"
    print(f"nowcast type:         {nowcast_type}")
    print(f"feature detector:     {feature_method}")
    print(f"extrapolator:         {extrap_method}")
    print(f"kernel type:          {kernel_type}")
    if add_perturbations and vel_pert_method is not None:
        print(f"velocity perturbator: {vel_pert_method}")
    print("")

    print("Parameters")
    print("----------")
    print(f"number of time steps:       {timesteps}")
    print(f"ARI model order:            {ari_order}")
    print(f"localization window radius: {localization_window_radius}")
    if add_perturbations:
        print(f"error dist. window radius:  {errdist_window_radius}")
        print(f"error ACF window radius:    {acf_window_radius}")
        print(f"ensemble size:              {n_ens_members}")
        print(f"parallel workers:           {num_workers}")
        print(f"seed:                       {seed}")
        if vel_pert_method == "bps":
            vp_par = vel_pert_kwargs.get(
                "p_par", noise.motion.get_default_params_bps_par()
            )
            vp_perp = vel_pert_kwargs.get(
                "p_perp", noise.motion.get_default_params_bps_perp()
            )
            print(
                f"velocity perturbations, parallel:      {vp_par[0]:.2f}, {vp_par[1]:.2f}, {vp_par[2]:.2f}"
            )
            print(
                f"velocity perturbations, perpendicular: {vp_perp[0]:.2f}, {vp_perp[1]:.2f}, {vp_perp[2]:.2f}"
            )
            vel_pert_kwargs = vel_pert_kwargs.copy()
            vel_pert_kwargs["vp_par"] = vp_par
            vel_pert_kwargs["vp_perp"] = vp_perp

    extrap_kwargs["allow_nonfinite_values"] = (
        True if np.any(~np.isfinite(precip)) else False
    )

    starttime_init = time.time()

    if check_norain(precip, 0.0, 0.0, None):
        return zero_precipitation_forecast(
            n_ens_members if nowcast_type == "ensemble" else None,
            timesteps,
            precip,
            callback,
            return_output,
            measure_time,
            starttime_init,
        )

    forecast_gen = _linda_deterministic_init(
        precip,
        velocity,
        feature_method,
        max_num_features,
        feature_kwargs,
        ari_order,
        kernel_type,
        localization_window_radius,
        extrap_method,
        extrap_kwargs,
        add_perturbations,
        num_workers,
        measure_time,
    )
    if measure_time:
        forecast_gen, precip_lagr_diff, init_time = forecast_gen
    else:
        forecast_gen, precip_lagr_diff = forecast_gen

    if add_perturbations:
        pert_gen = _linda_perturbation_init(
            precip,
            precip_lagr_diff,
            velocity,
            forecast_gen,
            pert_thrs,
            localization_window_radius,
            errdist_window_radius,
            acf_window_radius,
            vel_pert_method,
            vel_pert_kwargs,
            kmperpixel,
            timestep,
            num_workers,
            use_multiprocessing,
            measure_time,
        )
        if measure_time:
            precip_pert_gen, velocity_pert_gen, pert_init_time = pert_gen
            init_time += pert_init_time
        else:
            precip_pert_gen, velocity_pert_gen = pert_gen
    else:
        precip_pert_gen = None
        velocity_pert_gen = None

    precip_forecast = _linda_forecast(
        precip,
        precip_lagr_diff[1:],
        timesteps,
        forecast_gen,
        precip_pert_gen,
        velocity_pert_gen,
        n_ens_members,
        seed,
        measure_time,
        True,
        return_output,
        callback,
    )

    if return_output:
        if measure_time:
            return precip_forecast[0], init_time, precip_forecast[1]
        else:
            return precip_forecast
    else:
        return None
Generate a deterministic or ensemble nowcast by using the Lagrangian
INtegro-Difference equation model with Autoregression (LINDA) model.

Parameters
----------
precip: array_like
Array of shape (ari_order + 2, m, n) containing the input rain rate
or reflectivity fields (in linear scale) ordered by timestamp from
oldest to newest. The time steps between the inputs are assumed to be
regular.
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.
timesteps: int
Number of time steps to forecast.
feature_method: {'blob', 'domain' 'shitomasi'}
Feature detection method:

+-------------------+-----------------------------------------------------+
|    Method name    |                  Description                        |
+===================+=====================================================+
|  blob             | Laplacian of Gaussian (LoG) blob detector           |
|                   | implemented in scikit-image                         |
+-------------------+-----------------------------------------------------+
|  domain           | no feature detection, the model is applied over the |
|                   | whole domain without localization                   |
+-------------------+-----------------------------------------------------+
|  shitomasi        | Shi-Tomasi corner detector implemented in OpenCV    |
+-------------------+-----------------------------------------------------+

Default: 'blob'
max_num_features: int, optional
Maximum number of features to use. It is recommended to set this between
20 and 50, which gives a good tradeoff between localization and
computation time. Default: 25
feature_kwargs: dict, optional
Keyword arguments that are passed as ``**kwargs`` for the feature detector.
See :py:mod:`pysteps.feature.blob` and :py:mod:`pysteps.feature.shitomasi`.
ari_order: {1, 2}, optional
The order of the ARI(p, 1) model. Default: 1
kernel_type: {"anisotropic", "isotropic"}, optional
The type of the kernel. Default: 'anisotropic'
localization_window_radius: float, optional
The standard deviation of the Gaussian localization window.
Default: 0.2 * min(m, n)
errdist_window_radius: float, optional
The standard deviation of the Gaussian window for estimating the
forecast error distribution. Default: 0.15 * min(m, n)
acf_window_radius: float, optional
The standard deviation of the Gaussian window for estimating the
forecast error ACF. Default: 0.25 * min(m, n)
extrap_method: str, optional
The extrapolation method to use. See the documentation of
:py:mod:`pysteps.extrapolation.interface`. Default: 'semilagrangian'
extrap_kwargs: dict, optional
Optional dictionary containing keyword arguments for the extrapolation
method. See :py:mod:`pysteps.extrapolation.interface`.
add_perturbations: bool
Set to False to disable perturbations and generate a single
deterministic nowcast. Default: True
pert_thrs: float
Two-element tuple containing the threshold values for estimating the
perturbation parameters (mm/h). Default: (0.5, 1.0)
n_ens_members: int, optional
The number of ensemble members to generate. Default: 10
vel_pert_method: {'bps', None}, optional
Name of the generator to use for perturbing the advection field. See
:py:mod:`pysteps.noise.interface`. Default: 'bps'
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. For the default values and parameters optimized for different
domains, see :py:func:`pysteps.nowcasts.steps.forecast`.
kmperpixel: float, optional
Spatial resolution of the input data (kilometers/pixel). Required if
vel_pert_method is not None.
timestep: float, optional
Time step of the motion vectors (minutes). Required if vel_pert_method
is not None.
seed: int, optional
Optional seed for the random generators.
num_workers: int, optional
The number of workers to use for parallel computations. Applicable if
dask is installed. Default: 1
use_multiprocessing: bool, optional
Set to True to improve the performance of certain parallelized parts of
the code. If set to True, the main script calling linda.forecast must
be enclosed within the 'if __name__ == "__main__":' block.
Default: False
measure_time: bool, optional
If set to True, measure, print and return the computation time.
Default: False
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 precipitation fields, respectively. This can be used, for
instance, writing the outputs into files. Default: None
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. Default: True

Returns
-------
out: numpy.ndarray
A four-dimensional array of shape (n_ens_members, timesteps, m, n)
containing a time series of forecast precipitation fields for each
ensemble member. If add_perturbations is False, the first dimension is
dropped. The time series starts from t0 + timestep, where timestep is
taken from the input fields. 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). If return_output is set to False, a single None value is
returned instead.

Notes
-----
It is recommended to choose the feature detector parameters so that the
number of features is around 20-40. This gives a good tradeoff between
localization and computation time.

It is highly recommented to set num_workers>1 to reduce computation time.
In this case, 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.
pysteps/nowcasts/linda.py