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).

[mats-knmi]

We get the same errors semi-regularly, when there is only a small amount of precipitation in the model. There are 2 issues I see here:

1. The tapering method
   For us the tapering does seem to look alright (for a domain of 700x765 pixels). It looks to me like this "tukey" tapering method always procuces a perfect circle that fits the domain. In the case of @bwalraven (for a very rectangular domain) this is not a good tapering method, and I think this should maybe not be the default. Maybe we can stretch this tapering to match the input field shape (and produce an oval tapering).
   
2. The error handling
   In the solution of issue How to handle no rain in pysteps? #309 some checks were implemented for no rain in NWP and/or QPE. However, these checks are done before the tapering is applied. This results in the issue described here. In my opinion we should check if applying the tapering to self.__precip results in zero precip and then handle that the same as self.__params.zero_precip_radar. Next we should check if applying the tapering to self.__state.precip_noise_input results in zero precip and then handle that the same as self.__params.zero_precip_model_fields (i.e. returning self.__zero_precipitation_forecast()). I will create a PR with this proposed solution. Let me know what you think @dnerini. (but also @ladc and @RubenImhoff, since you worked on How to handle no rain in pysteps? #309).

[RubenImhoff]

Thanks, @mats-knmi! I think combining both points that you mention is a good way to find a solution for this issue. :) Happy to give the proposed solution a review.

[mats-knmi]

I created a PR: #460

[SKB-7]

I'm also facing the same issue as mentioned in #451, " is an array of zeroes, and thus N in the following lines in [noise generation] returns a nan."

I have seen PR: #460 created by @mats-knmi. It is not completely merged yet, right? Does something need to be fixed? I'm happy to get involved.

Many Thanks.

[mats-knmi]

Hi @SKB-7, thanks for your enthusiasm. Technically #460 is ready to be merged. I am just waiting for @dnerini to answer, since he first picked up this issue.

It would be helpful for us to know if that PR solves your issue though. Can you try it out and see if it works? You could try to do the following to test that PR out.

pip install git+https://github.com/pySTEPS/pysteps.git@fix-noise-tapering-and-error-handling

[SKB-7]

Hi @mats-knmi, After installing the branch pysteps.git@fix-noise-tapering-and-error-handling, I have tested my code, but still, noise generation using initialize_nonparam_2d_ssft_filter and generate_noise_2d_ssft_filter returns a nan values.

[mats-knmi]

@SKB-7 Interesting, which part of pysteps do you run? My PR should fix this issue for the steps blending. Do you maybe have an example of how you run pysteps, so that I can reproduce this?

[SKB-7]

Below is the function I'm using to generate noise for n ensembles.

def getCorrNoise(n, h, dt, obsFile, U_windFile, V_windFile, seednum):

ds1 = xr.open_dataset(U_windFile)
ds2 = xr.open_dataset(V_windFile)

u_values = ds1['u'][0:(0+n),:,:] # eastward wind [m/s]
v_values = ds2['v'][0:(0+n),:,:] # northward wind [m/s]

# ----  open precip data  ----
ds = xr.open_dataset(obsFile)

har = ds['prcp'][0:(0+n),:,:].data
har[har<0.1] = 0 
ds=None
xsize=har.shape[2]
ysize=har.shape[1]

yind=np.repeat(np.arange(0,ysize),xsize).reshape(ysize,xsize).astype('int16')
xind=np.tile(np.arange(0,xsize),ysize).reshape(ysize,xsize).astype('int16')

dx = dy = 10000

# create array to hold simulated noise
s = np.empty((n,ysize,xsize),dtype=np.float32)

# --- SIMULATE FIRST FIELD ---
if len(har[0,:,:][har[0,:,:]>=0.1]) > 1133:          
    # print("entered 1st if")
    # only use pysteps smoothing filter if there is actually rain in the study area
    Fnp = initialize_nonparam_2d_ssft_filter(har[0,:,:])#, win_size=(64,64)) 
    s[0,:,:] = generate_noise_2d_ssft_filter(Fnp, seed=seednum)
    firstRain=True
    
else:
    # print("entered 1st else")
    # otherwise just use smoothed white noise
    s[0,:,:] = np.random.normal(size=(ysize,xsize)) 
    #**seed kwargs are not there in normal fnc so commented,seed=seednum)
    firstRain=False

last_rn = s[0,:,:] # record this instance of random noise as most recent instance of random noise


# set initial value of alpha parameter - determines degree to which noise is incorporated into field at each timestep
alpha = 0.95

# --- ADVECT FIELD FORWARD IN TIME ---
for hr in range(1,n):
    # print(f"hr: {hr}")
    
    # get wind values from PREVIOUS DAY and calculate dx, dy that correspond to this time step
    u = u_values[hr-1,:,:] # eastward wind [m/s]
    deltax = u*60.*60.*24. # distance wind travels each time step in positive x dir [m]
    dix = deltax//dx # number of grid steps to shift in semi-lagrangian scheme in units of [0.1 deg pixels]
    
    v = v_values[hr-1,:,:] # northward wind [m/s]
    deltay = -v*60.*60.*24. # distance wind travels each time step in positive y dir [m]
    diy = deltay//dy # number of grid steps to shift in semi-lagrangian scheme [0.1 deg pixels]
    
    harh = har[hr,:,:] 
    
    # use pySTEPS to generate noise with power spectrum of Precip field
    # only use pysteps smoothing filter if there is actually rain in the study area
    
    if len(harh[harh>=0.1]) > 1133: # atleast 5% of study area must be rainy
        # print("entered 2nd if")
        
        Fnp = initialize_nonparam_2d_ssft_filter(harh)#, win_size=(64,64))
        rn = generate_noise_2d_ssft_filter(Fnp, seed=seednum+hr)
        # print(f'seednum {seednum+hr}')
        # get average temporal correlation of field in previous timestep
        alpha = getAlpha(harh,har[hr-1,:,:])
       
        firstRain=True            
    else:
        # print("entered 2nd else")
        
        # if field does not have at least 5% rainy pixels, generate random field 
        # using spatial correlation structure of last instance of rainfall 
        
        if firstRain==True:
            print("entered 3rd if")                
            rn = generate_noise_2d_ssft_filter(Fnp, seed=seednum+hr)
        
        else:
            # print("entered 3rd else")
            # if first instance of a "rainy field" hasn't occurred yet in study period, use white noise
            rn = np.random.normal(size=(ysize,xsize))
            # **seed kwargs are not there in normal fnc so commented,seed=seednum+hr) # white noise
            # print(f'random noise {hr}')
                    
    
    if np.isnan(rn[0,0])==True: # if pysteps smoothing process doesn't work for some reason, use last realization of smoothed noise
        print("NaN in noise field at timestep %d, using last realization of noise field"%hr)
        rn = last_rn
        # print(f'error ssft_filter {hr}')
    else:
        last_rn = rn
    

    # -- advect noise values from previous time step & perturb with noise field rn --
    ybefore=np.array((yind-diy)%ysize,dtype='int16')
    xbefore=np.array((xind-dix)%xsize,dtype='int16')
    s[hr,:] = (alpha*s[hr-1,ybefore,xbefore] + np.sqrt(1.-alpha**2)*rn).astype('float32')
    # print(f"hr: {hr} ,Size of s: {s.shape}")

# return noise field
return(s)

[mats-knmi]

@SKB-7 I think I understand the issue now. I have fixed this for the steps blending, not for any of the nowcasts (and also not if you call these noise generation functions directly, like you are). Re-reading the issue I see that @bwalraven probably also had the issue with the nowcasts and not the blending.

In these cases there is currently no check for no rain built into the methods themselves, so you have to write them yourself before calling the method. I see that you have already implemented this in your code if len(harh[harh>=0.1]) > 1133, but this does not cover the case where there is only precipitation at the edges.

Basically the way I solved this for the blending is by improving on that check. I suppose a solution would be to add this improved norain check as a util function to the pysteps code, so that you can use this to check if you can call pysteps noise generation, so that you can decide on a fallback.

Also, just out of curiosity, what are you using pysteps for? If you are generating an ensemble nowcast, why are you not using the steps nowcast (or any of the other nowcast modules)? And if you are creating your own improved nowcast method, would you maybe consider contributing your alternative method to the pysteps project?