mpvginde / pysteps

Python framework for short-term ensemble prediction systems.
https://pysteps.github.io/
BSD 3-Clause "New" or "Revised" License
1 stars 0 forks source link

Correct noise cascade initialisation #6

Closed mpvginde closed 1 year ago

mpvginde commented 1 year ago

Task related to this pySTEPS issue: 332

Now, the all noise cascade levels are initialized to zero. It is probably better to initialize the noise with the appropriate noise field as provided by the EPS term.

The noise cascade has holds 2 time levels of noise. Both time levels should be initialized with the same noise field. If not the overall power of the noise cascade level will not converge to the correct power.

mpvginde commented 1 year ago

Created a new helper function in the pysteps.blending.steps.py file which initializes each level of the noise cascade with noise having the correct spatial correlation. All ar_order steps have identical noise to ensure that noise cascade immediately has the correct power in the subsequent timesteps (no spinup needed).

def _init_noise_cascade(
    shape,
    n_ens_members,
    n_cascade_levels,
    generate_noise,
    decompositor,
    pp,
    randgen_prec,
    fft_objs,
    bp_filter,
    domain,
    noise_method,
    noise_std_coeffs,
    ar_order
):
    noise_cascade = np.zeros(shape)
    if noise_method:
        for j in range(n_ens_members):        
            EPS = generate_noise(
                pp,
                randstate=randgen_prec[j],
                fft_method=fft_objs[j],
                domain=domain
                )
            EPS = decompositor(
                EPS,
                bp_filter,
                fft_method=fft_objs[j],
                input_domain=domain,
                output_domain=domain,
                compute_stats=True,
                normalize=True,
                compact_output=True
                )
            for i in range(n_cascade_levels):
                EPS_ = EPS["cascade_levels"][i]
                EPS_ *= noise_std_coeffs[i]
                for n in range(ar_order):
                    noise_cascade[j][i][n] = EPS_
            EPS = None
            EPS_ = None
    return(noise_cascade)
mpvginde commented 1 year ago

The noise was previously initialized before the initialization of the random generator. But since you now need them to initialize the noise. The noise_cascade = _init_noise_cascade(...) should be called after the random generator initializations. Also check with developments in issue #1

mpvginde commented 1 year ago

Fixed with commit 2eb125bc89f4cd13205bd88ff763d6a0097c938e

RubenImhoff commented 1 year ago

This is excellent!!