Initializing noise field with default arguments returns array of zeros. Is this expected behavior?

[bwalraven]

I am running a probabilistic nowcast for a rectangular precipitation field (199x113) using the steps function with the following settings:

domain = 'spatial'
extrap_method = 'semilagrangian'
decomp_method = 'fft'
fft_method = 'numpy'
noise_method = 'nonparametric'
noise_stddev_adj = 'auto'

and run into the following error before the noise adjustment coefficients are computed:
ValueError: field contains non-finite values.

This points me to the following statement in the cascade decomposition:

if np.any(~np.isfinite(field)):
    raise ValueError("field contains non-finite values") 

which in turn points me to decomp_N = decomp_method(N, F, mask=MASK_) in the following function:

pysteps/pysteps/noise/utils.py

Line 24 in e332585

def compute_noise_stddev_adjs(

This stems from the fact that the noise field F returned in the function:

pysteps/pysteps/noise/fftgenerators.py

Line 214 in e332585

def initialize_nonparam_2d_fft_filter(field, **kwargs):

is an array of zeroes, and thus N in the following lines in noise generation returns a nan.

It appears the issue lies in the following lines in the function

pysteps/pysteps/noise/fftgenerators.py

Line 214 in e332585

def initialize_nonparam_2d_fft_filter(field, **kwargs):

F = np.zeros(fft_shape, dtype=complex)
 for i in range(nr_fields):
     if use_full_fft:
         F += fft.fft2(field[i, :, :] * tapering)
     else:
         F += fft.rfft2(field[i, :, :] * tapering)
 F /= nr_fields

Using the default arguments, with window_function="tukey" and a field size of 199x113, the tapering returns an oval shape (black for values >0):

This essentially means that if we have a precipitation field, with precipitation only at the edges, and use the default steps settings, this would return a noise field of zeros and eventually lead to an error in the cascade decomposition.

Though the function appears to numerically behave as expected, is this an error or is this desired behavior? Does this mean it's not possible to compute a nowcast if we have precipitation only around the edges of our domain?

[bwalraven]

Adding a small piece of code as example:

# Initialize zero array
rainy_edges = np.zeros((3, 199, 113))

# Disturb the top 20 and bottom 20 rows for each field with values between 1-10
for i in np.arange(0, 3):
    rainy_edges[i, :20, :] = np.random.randint(1, 10, (20, 113))
    rainy_edges[i, -20:, :] = np.random.randint(1, 10, (20, 113))


# Set default arguments
field = rainy_edges
use_full_fft = 'False' # default
win_fun = 'tukey' # default

# Copied code from the non-parametric noise initializer
# dims
if len(field.shape) == 2:
    field = field[None, :, :]
nr_fields = field.shape[0]
field_shape = field.shape[1:]
if use_full_fft:
    fft_shape = (field.shape[1], field.shape[2])
else:
    fft_shape = (field.shape[1], int(field.shape[2] / 2) + 1)

# make sure non-rainy pixels are set to zero
field -= field.min(axis=(1, 2))[:, None, None]

if win_fun is not None:
    tapering = utils.tapering.compute_window_function(
        field_shape[0], field_shape[1], win_fun
    )
else:
    tapering = np.ones(field_shape)

F = np.zeros(fft_shape, dtype=complex)
for i in range(nr_fields):
    if use_full_fft:
        F += fft.fft2(field[i, :, :] * tapering)
    else:
        F += fft.rfft2(field[i, :, :] * tapering)
F /= nr_fields

# returns a field with only zeros
print(F.real.sum())

# plotting the tapering
plt.figure(figsize=(6, 10))
plt.title('Tapering')
plt.pcolormesh(tapering>0, cmap='Greys')

[dnerini]

Hi @bwalraven thanks for the detailed issue report.

I'll try to have a close look, but in essence the tapering is probably behaving as expected as we exclude rainy pixels at the edges when computing the power spectra to avoid artefacts.

However, the plot of the tapering function does look a bit weird, and, more importantly, it should handle zero fields more graciously instead of raising an error.

[dnerini]

I also wonder if this is not an edge case that was not covered when addressing #309 @ladc

[bwalraven]

Hi @dnerini thanks for looking into this. Curious to what you can find out.

And understandable to exclude some pixels at the edges, though what strikes me in the example I showed is the difference in excluded pixels in x and y direction (approx. 1 pixel vs. almost 50 pixels). Note that the image I included is a little bit distorted as there are more y-pixels (199) than x-pixels (113).