Applying the check_norain function to all timesteps

[bwalraven]

In PR #460 the possibility to add a check_norain to a steps nowcast was added, which is very nice!

Currently, however, the check_norain function checks if the entire precip array is above the no_rain threshold or not. This means that when it is given a 3d precip array (time, lat, lon), it can pass the no_rain check, and still fail later on if one of the timesteps contains a zero precip field but the total of all time steps is above the threshold.

Example: the following dummy precip array passes the no_rain check:

import numpy as np
import pysteps
from pysteps import utils

# Create a dummy 3d precip array of which one slice has zero precip
np.random.seed(24)

slice1 = np.random.rand(100, 100)
slice2 = np.random.rand(100, 100)
slice3 = np.random.rand(100, 100)
slice4 = np.zeros((100, 100))

precip_arr = np.stack([slice1, slice2, slice3, slice4], axis=0)

precip_thr = 0.1
norain_thr = 0.03
win_fun = 'tukey'

norain = utils.check_norain.check_norain(precip_arr[-3:, :, :], precip_thr, norain_thr, win_fun)

print(norain)

However, if we use this array of shape (4, 100, 100) to compute a motion field and then want to perform the nowcast using the following code:

# First compute the motion field
motion_method = 'LK'
oflow_method = pysteps.motion.get_method(motion_method)
velocity = oflow_method(precip_arr)

nwc_method = pysteps.nowcasts.get_method("steps")
precip_forecast = nwc_method(precip_arr[-3:, :, :],
                             velocity=velocity, 
                             timesteps=12,
                             norain_thr=0.03,
                             precip_thr=0.1, 
                             kmperpixel=1.0,
                             timestep=5)

the following error is returned pointing to the temporal autocorrelation function

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 7
      4 velocity = oflow_method(precip_arr)
      6 nwc_method = pysteps.nowcasts.get_method("steps")
----> 7 precip_forecast = nwc_method(precip_arr[-3:, :, :],
      8                              velocity=velocity, 
      9                              timesteps=12,
     10                              norain_thr=0.03,
     11                              precip_thr=0.1, 
     12                              kmperpixel=1.0,
     13                              timestep=5)

File ~/.conda/envs/pysteps/lib/python3.10/site-packages/pysteps/nowcasts/steps.py:1533, in forecast(precip, velocity, timesteps, n_ens_members, n_cascade_levels, precip_thr, norain_thr, kmperpixel, timestep, extrap_method, decomp_method, bandpass_filter_method, noise_method, noise_stddev_adj, ar_order, vel_pert_method, conditional, probmatching_method, mask_method, seed, num_workers, fft_method, domain, extrap_kwargs, filter_kwargs, noise_kwargs, vel_pert_kwargs, mask_kwargs, measure_time, callback, return_output)
   1529 # Create an instance of the new class with all the provided arguments
   1530 nowcaster = StepsNowcaster(
   1531     precip, velocity, timesteps, steps_config=nowcaster_config
   1532 )
-> 1533 forecast_steps_nowcast = nowcaster.compute_forecast()
   1534 nowcaster.reset_states_and_params()
   1535 # Call the appropriate methods within the class

File ~/.conda/envs/pysteps/lib/python3.10/site-packages/pysteps/nowcasts/steps.py:377, in StepsNowcaster.compute_forecast(self)
    366     return zero_precipitation_forecast(
    367         self.__config.n_ens_members,
    368         self.__time_steps,
   (...)
    373         self.__start_time_init,
    374     )
    376 self.__perform_extrapolation()
--> 377 self.__apply_noise_and_ar_model()
    378 self.__initialize_velocity_perturbations()
    379 self.__initialize_precipitation_mask()

File ~/.conda/envs/pysteps/lib/python3.10/site-packages/pysteps/nowcasts/steps.py:832, in StepsNowcaster.__apply_noise_and_ar_model(self)
    827 self.__params.autocorrelation_coefficients = np.empty(
    828     (self.__config.n_cascade_levels, self.__config.ar_order)
    829 )
    830 for i in range(self.__config.n_cascade_levels):
    831     self.__params.autocorrelation_coefficients[i, :] = (
--> 832         correlation.temporal_autocorrelation(
    833             self.__state.precip_cascades[i], mask=self.__state.mask_threshold
    834         )
    835     )
    837 nowcast_utils.print_corrcoefs(self.__params.autocorrelation_coefficients)
    839 # Adjust the lag-2 correlation coefficient if AR(2) model is used

File ~/.conda/envs/pysteps/lib/python3.10/site-packages/pysteps/timeseries/correlation.py:107, in temporal_autocorrelation(x, d, domain, x_shape, mask, use_full_fft, window, window_radius)
    102     raise ValueError(
    103         "dimension mismatch between x and mask: x.shape[1:]=%s, mask.shape=%s"
    104         % (str(x.shape[1:]), str(mask.shape))
    105     )
    106 if np.any(~np.isfinite(x)):
--> 107     raise ValueError("x contains non-finite values")
    109 if d == 1:
    110     x = np.diff(x, axis=0)

ValueError: x contains non-finite values

To fix this we could modify the check_norain function slightly like:

    if win_fun is not None:
        tapering = utils.tapering.compute_window_function(
            precip_arr.shape[-2], precip_arr.shape[-1], win_fun
        )
    else:
        tapering = np.ones((precip_arr.shape[-2], precip_arr.shape[-1]))

    tapering_mask = tapering == 0.0
    masked_precip = precip_arr.copy()
    masked_precip[..., tapering_mask] = np.nanmin(precip_arr)

    if precip_thr is None:
        precip_thr = np.nanmin(masked_precip)
    
    if precip_arr.ndim == 3:
        n_fields = precip_arr.shape[0]
        for i in range(n_fields):
            precip_field = masked_precip[i]
            rain_pixels = precip_field[precip_field > precip_thr]
            norain = rain_pixels.size / precip_field.size <= norain_thr
            print(
                  f"Rain fraction in precip_field {i} is: {str(rain_pixels.size / precip_field.size)}, while minimum fraction is {str(norain_thr)}"
              )
            if norain == True:
            # all fields should be False, if any is True return True
                    return True
    elif precip_arr.ndim == 2:
        rain_pixels = masked_precip[masked_precip > precip_thr]
        norain = rain_pixels.size / masked_precip.size <= norain_thr
        print(
            f"Rain fraction is: {str(rain_pixels.size / masked_precip.size)}, while minimum fraction is {str(norain_thr)}"
        )
        return norain

[dnerini]

hi @bwalraven thanks for the detailed issue, perhaps @mats-knmi could take this up since you actively work on the check_norain method?

@bwalraven would you be willing to contribute with a PR to introduce the changes you suggest?

[mats-knmi]

Hi, just got back from my holiday, so I can take a look.

@bwalraven I am not sure if we need to change the check norain method or the way this method is called in the steps nowcast (and possibly also the blending?). The check norain method is used for both the model data (in steps blending) and for the radar data (in both the blending and the nowcast).
I am not 100% on this but I am fairly certain that for the model data you only want this to return true if all members and all timesteps contain no rain.
For the steps nowcast it does indeed seem to be the case that if just one of the (typically) 4 input frames has norain then the processing will fail, so maybe we need to update the nowcast code to call the check norain method for each input frame?

I am not completely sure how this holds up in the blending though, I think especially it would be a loss there if a single empty frame in the entire model forecast would result in not using the forecast at all. But it also seems like it uses the radar input data in a similar way, so maybe the norain check there needs to be updated as well?

@RubenImhoff do you have any ideas about this?

[RubenImhoff]

Hi all, it indeed seems to be helpful to perform the no_rain check on all the input frames instead of the total at once. I also agree with @mats-knmi that this is not warranted for the NWP data, so maybe we should have a fix in which we feed it one frame at a time from the radar data to perform the check? I can even imagine that if we feed it more than the minimum number of frames required by the optical flow algorithm, that we can give the option to work with a reduced number of provided frames (given that the reduced number of frames are all above the threshold). What do you think?

[bwalraven]

Hi, sorry this got snowed under a bit..
@mats-knmi thanks for taking a look. Personally I have only used pysteps for nowcasting (radar) data, so you probably know best how the no_rain check is used elsewhere for NWP data and blending. Both ways you propose would solve the issue I mentioned, so if you say calling the check_no_rain function for each input frame in the steps nowcast is preferred, rather than changing the check_no_rain function itself, then that seems fine to me.