Implement no-rain check on radar fields/set default autocorrelation parameters

[mpvginde]

Task related to this pySTEPS issue 309

From the A. Seed comments as quoted by @ladc:

• Set threshold at around 3% raining fraction.
• Default auto-correlations [code from STEPS1] below this fraction. ▶ Default (lag2) auto-correlations use a power law for the correlation length as a function of scale + black magic

[RubenImhoff]

With no rain in the precip radar fields, the noise cannot be initialized. The proposed change is to initialize the noise in step 3 instead of step 2 (hence, after the cascade decompositions). This allows for a check if there is no rain in both the radar field and NWP fields. If that is the case, the resulting forecast will directly consist of only zeroes and thus the noise initialization and all subsequent steps will be bypassed. Otherwise, if the radar field is zero and the NWP field is not, the initialization of noise will take place with the NWP field at timestep zero (which is at that time step zero in 'Lagrangian coordinates').

[mpvginde]

Link with issue #6

[RubenImhoff]

Changed the auto-correlation derivation function to make sure it gives the climatological auto-correlation values when the radar rainfall fields are zero:

def _estimate_ar_parameters_radar(
    R_c, ar_order, n_cascade_levels, MASK_thr, zero_precip_radar
):
    """Estimate AR parameters for the radar rainfall field."""
    # If there are values in the radar fields, compute the autocorrelations
    GAMMA = np.empty((n_cascade_levels, ar_order))
    if not zero_precip_radar:
        # compute lag-l temporal autocorrelation coefficients for each cascade level
        for i in range(n_cascade_levels):
            GAMMA[i, :] = correlation.temporal_autocorrelation(R_c[i], mask=MASK_thr)

    # Else, use standard values for the autocorrelations
    else:
        # Get the climatological lag-1 and lag-2 autocorrelation values from Table 2
        # in `BPS2004`.
        # Hard coded, change to own (climatological) values when present.
        GAMMA = np.array(
            [
                [0.99805, 0.9925, 0.9776, 0.9297, 0.796, 0.482, 0.079, 0.0006],
                [0.9933, 0.9752, 0.923, 0.750, 0.367, 0.069, 0.0018, 0.0014],
            ]
        )

        # Check whether the number of cascade_levels is correct
        if GAMMA.shape[1] > n_cascade_levels:
            GAMMA = GAMMA[:, 0:n_cascade_levels]
        elif GAMMA.shape[1] < n_cascade_levels:
            # Get the number of cascade levels that is missing
            n_extra_lev = n_cascade_levels - GAMMA.shape[1]
            # Append the array with correlation values of 10e-4
            GAMMA = np.append(
                GAMMA,
                [np.repeat(0.0006, n_extra_lev), np.repeat(0.0014, n_extra_lev)],
                axis=1,
            )

        # Finally base GAMMA.shape[0] on the AR-level
        if ar_order == 1:
            GAMMA = GAMMA[0, :]
        if ar_order > 2:
            for repeat_index in range(ar_order - 2):
                GAMMA = np.vstack((GAMMA, GAMMA[1, :]))

        # Finally, transpose GAMMA to ensure that the shape is the same as np.empty((n_cascade_levels, ar_order))
        GAMMA = GAMMA.transpose()
        assert GAMMA.shape == (n_cascade_levels, ar_order)

    # Print the GAMMA value
    nowcast_utils.print_corrcoefs(GAMMA)

    if ar_order == 2:
        # adjust the lag-2 correlation coefficient to ensure that the AR(p)
        # process is stationary
        for i in range(n_cascade_levels):
            GAMMA[i, 1] = autoregression.adjust_lag2_corrcoef2(GAMMA[i, 0], GAMMA[i, 1])

    # estimate the parameters of the AR(p) model from the autocorrelation
    # coefficients
    PHI = np.empty((n_cascade_levels, ar_order + 1))
    for i in range(n_cascade_levels):
        PHI[i, :] = autoregression.estimate_ar_params_yw(GAMMA[i, :])

    nowcast_utils.print_ar_params(PHI)
    return PHI

[RubenImhoff]

Fixed with commit https://github.com/mpvginde/pysteps/commit/b375e69b242bd71a5e7349fbaf371ff5287c43a5

[ladc]

A possible improvement here is to replace the climatological values from the paper by the ones calculated for our domain.