[nowcasts.steps] vel_pert_method="bps" returns NaNs when motion is zero

[dnerini]

With scipy>=1.6, a call to steps.forecast with arguments vel_pert_method="bps" and V=np.zeros((2, R[0, :].shape*)) returns all NaNs.

Steps to reproduce

Set up test environment:

conda create -n test_py38 python=3.8
conda activate test_py38
conda install -c conda-forge cython numpy jsmin jsonschema matplotlib netCDF4 opencv pillow pyproj pygrib scipy dask pyfftw cartopy h5py PyWavelets pandas scikit-image pytest pytest-cov
pip install git+https://github.com/pySTEPS/pysteps.git@master

In python:

import numpy as np

from pysteps import nowcasts
from pysteps.tests.helpers import get_precipitation_fields

precip_input, metadata = get_precipitation_fields(
    num_prev_files=2,
    num_next_files=0,
    return_raw=False,
    metadata=True,
    upscale=2000,
)
precip_input = precip_input.filled()
field_shape = (precip_input.shape[1], precip_input.shape[2])
motion_field = np.zeros((2, *field_shape))

precip_output = nowcasts.get_method("steps")(
    precip_input,
    motion_field,
    timesteps=3,
    R_thr=metadata["threshold"],
    kmperpixel=2.0,
    timestep=metadata["accutime"],
    seed=42,
    n_ens_members=2,
    vel_pert_method="bps",
)

assert np.isnan(precip_output).all()

Expected output

precip_output should be an array of finite values.

Cause

The output of steps.forecast being all NaNs is caused by the call to extrapolator_method(R_f_ip, V_pert, [t_diff_prev], **extrap_kwargs_) (L#757 of steps.py) (Note that this is calling semilagrangian.extrapolate) with V_pert being all NaNs.

This is caused by the call to generate_vel_noise(vps[j], t_total[j] * timestep) (L#754 in steps.py) at lines N = linalg.norm(V, axis=0) and V_n = V / np.stack([N, N]) (L#127-128 in noise/motion.py) when variable V is all zeros.

I can see two main problems here that would need to be fixed:

1. Call to generate_vel_noise(vps[j], t_total[j] * timestep) should not return NaNs when zero motion is passed, but rather zero or little velocity noise.
2. Call to semilagrangian.extrapolate(R_f_ip, V_pert, [t_diff_prev], **extrap_kwargs_) should not silently fail when V_pert is non-finite, but rather raise an exception.

Originally posted by @dnerini in https://github.com/pySTEPS/pysteps/issues/202#issuecomment-757163327

[dnerini]

@pulkkins do you want to take this over or should I give it a try first?

[pulkkins]

I'll take care of this.

[pulkkins]

This should be fixed in commits 9e7a951, c46a001 and b4198d3. Commit b4198d3 also fixed the problem that the test for the velocity field was inside wrong if clause.

[dnerini]

Very good, @pulkkins , many thanks for taking care of this!