How to use the pysteps.cascade.get_method function in pysteps

print("number of time steps:     %d" % n_timesteps)
    print("parallel threads:         %d" % num_workers)
    print("number of cascade levels: %d" % n_cascade_levels)
    print("order of the AR(p) model: %d" % ar_order)
    print("precip. intensity threshold: %g" % R_thr)

    if measure_time:
        starttime_init = time.time()

    fft = utils.get_method(fft_method, shape=R.shape[1:],
                           n_threads=num_workers)

    M, N = R.shape[1:]

    # initialize the band-pass filter
    filter_method = cascade.get_method(bandpass_filter_method)
    filter = filter_method((M, N), n_cascade_levels, **filter_kwargs)

    decomp_method = cascade.get_method(decomp_method)

    extrapolator_method = extrapolation.get_method(extrap_method)

    R = R[-(ar_order + 1):, :, :].copy()
    R_min = R.min()

    if conditional:
        MASK_thr = np.logical_and.reduce([R[i, :, :] >= R_thr for i in range(R.shape[0])])
    else:
        MASK_thr = None

    # initialize the extrapolator
    x_values, y_values = np.meshgrid(np.arange(R.shape[2]),

if conditional or mask_method is not None:
        print("precip. intensity threshold: %g" % R_thr)

    num_ensemble_workers = n_ens_members if num_workers > n_ens_members \
        else num_workers

    if measure_time:
        starttime_init = time.time()

    fft = utils.get_method(fft_method, shape=R.shape[1:], n_threads=num_workers)

    M, N = R.shape[1:]

    # initialize the band-pass filter
    filter_method = cascade.get_method(bandpass_filter_method)
    filter = filter_method((M, N), n_cascade_levels, **filter_kwargs)

    decomp_method = cascade.get_method(decomp_method)

    extrapolator_method = extrapolation.get_method(extrap_method)

    x_values, y_values = np.meshgrid(np.arange(R.shape[2]),
                                     np.arange(R.shape[1]))

    xy_coords = np.stack([x_values, y_values])

    R = R[-(ar_order + 1):, :, :].copy()

    if conditional:
        MASK_thr = np.logical_and.reduce([R[i, :, :] >= R_thr for i in range(R.shape[0])])
    else:

R_min = metadata["zerovalue"]

    num_ensemble_workers = n_ens_members if num_workers > n_ens_members else num_workers

    if measure_time:
        starttime_init = time.time()

    # get methods
    extrapolator_method = extrapolation.get_method(extrap_method)

    x_values, y_values = np.meshgrid(np.arange(R.shape[2]), np.arange(R.shape[1]))

    xy_coords = np.stack([x_values, y_values])

    decomp_method = cascade.get_method(decomp_method)
    filter_method = cascade.get_method(bandpass_filter_method)
    if noise_method is not None:
        init_noise, generate_noise = noise.get_method(noise_method)

    # advect the previous precipitation fields to the same position with the
    # most recent one (i.e. transform them into the Lagrangian coordinates)
    R = R[-(ar_order + 1) :, :, :].copy()
    extrap_kwargs = extrap_kwargs.copy()
    extrap_kwargs["xy_coords"] = xy_coords
    res = []
    f = lambda R, i: extrapolator_method(
        R[i, :, :], V, ar_order - i, "min", **extrap_kwargs
    )[-1]
    for i in range(ar_order):
        if not dask_imported:
            R[i, :, :] = f(R, i)
        else:

num_ensemble_workers = n_ens_members if num_workers > n_ens_members \
        else num_workers

    if measure_time:
        starttime_init = time.time()

    fft = utils.get_method(fft_method, shape=R.shape[1:], n_threads=num_workers)

    M, N = R.shape[1:]

    # initialize the band-pass filter
    filter_method = cascade.get_method(bandpass_filter_method)
    filter = filter_method((M, N), n_cascade_levels, **filter_kwargs)

    decomp_method = cascade.get_method(decomp_method)

    extrapolator_method = extrapolation.get_method(extrap_method)

    x_values, y_values = np.meshgrid(np.arange(R.shape[2]),
                                     np.arange(R.shape[1]))

    xy_coords = np.stack([x_values, y_values])

    R = R[-(ar_order + 1):, :, :].copy()

    if conditional:
        MASK_thr = np.logical_and.reduce([R[i, :, :] >= R_thr for i in range(R.shape[0])])
    else:
        MASK_thr = None

    # advect the previous precipitation fields to the same position with the

F = abs(np.fft.fftshift(np.fft.fft2(R_)))
fig = plt.figure()
M,N = F.shape
im = plt.imshow(np.log(F**2), vmin=4, vmax=24, cmap=cm.jet, 
                extent=(-N/2, N/2, -M/2, M/2))
cb = fig.colorbar(im)
plt.xlabel("Wavenumber $k_x$")
plt.ylabel("Wavenumber $k_y$")
plt.title("Log-power spectrum of R")
plt.show()

# Cascade decomposition

## construct the Gaussian bandpass filter
bandapass_filter = stp.cascade.get_method("gaussian")
filter = bandapass_filter(R_.shape, num_cascade_levels, gauss_scale=0.5, 
                          gauss_scale_0=0.5)

## plot the bandpass filter weights
fig = plt.figure()
ax = fig.gca()
L = max(N, M)
for k in range(num_cascade_levels):
    ax.semilogx(np.linspace(0, L/2, len(filter["weights_1d"][k, :])), 
                filter["weights_1d"][k, :], "k-", 
                basex=pow(0.5*L/3, 1.0/(num_cascade_levels-2)))
                
ax.set_xlim(1, L/2)
ax.set_ylim(0, 1)
ax.get_xaxis().set_major_formatter(ticker.ScalarFormatter())
xt = np.hstack([[1.0], filter["central_wavenumbers"][1:]])

filter["weights_1d"][k, :], "k-", 
                basex=pow(0.5*L/3, 1.0/(num_cascade_levels-2)))
                
ax.set_xlim(1, L/2)
ax.set_ylim(0, 1)
ax.get_xaxis().set_major_formatter(ticker.ScalarFormatter())
xt = np.hstack([[1.0], filter["central_wavenumbers"][1:]])
ax.set_xticks(xt)
ax.set_xticklabels(["%.2f" % cf for cf in filter["central_wavenumbers"]])
ax.set_xlabel("Radial wavenumber $|\mathbf{k}|$")
ax.set_ylabel("Normalized weight")
ax.set_title("Bandpass filter weights")
plt.show()

## compute the cascade decomposition
decomposition = stp.cascade.get_method("fft")
cascade = decomposition(R_, filter)

## plot the normalized cascade levels (mean zero and standard deviation one)
grid_res_km = max(metadata["xpixelsize"], metadata["ypixelsize"])/1000.
mu,sigma = cascade["means"],cascade["stds"]
nrows = int(np.ceil((1+num_cascade_levels)/4.))
plt.subplot(nrows,4,1)
for k in range(num_cascade_levels+1):
    if k==0:
        plt.subplot(nrows,4,k+1)
        stp.plt.plot_precip_field(R, units=unit, title="Rainfall field", colorbar=False)
    else:
        R_k = cascade["cascade_levels"][k-1, :, :]
        R_k = (R_k - mu[k-1]) / sigma[k-1]
        plt.subplot(nrows,4,k+1)
        im = plt.imshow(R_k, cmap=cm.jet, vmin=-6, vmax=6)

print("order of the AR(p) model: %d" % ar_order)
    print("precip. intensity threshold: %g" % R_thr)

    if measure_time:
        starttime_init = time.time()

    fft = utils.get_method(fft_method, shape=R.shape[1:],
                           n_threads=num_workers)

    M, N = R.shape[1:]

    # initialize the band-pass filter
    filter_method = cascade.get_method(bandpass_filter_method)
    filter = filter_method((M, N), n_cascade_levels, **filter_kwargs)

    decomp_method = cascade.get_method(decomp_method)

    extrapolator_method = extrapolation.get_method(extrap_method)

    R = R[-(ar_order + 1):, :, :].copy()
    R_min = R.min()

    if conditional:
        MASK_thr = np.logical_and.reduce([R[i, :, :] >= R_thr for i in range(R.shape[0])])
    else:
        MASK_thr = None

    # initialize the extrapolator
    x_values, y_values = np.meshgrid(np.arange(R.shape[2]),
                                     np.arange(R.shape[1]))

    xy_coords = np.stack([x_values, y_values])

R_thr = metadata["threshold"]
    R_min = metadata["zerovalue"]

    num_ensemble_workers = n_ens_members if num_workers > n_ens_members else num_workers

    if measure_time:
        starttime_init = time.time()

    # get methods
    extrapolator_method = extrapolation.get_method(extrap_method)

    x_values, y_values = np.meshgrid(np.arange(R.shape[2]), np.arange(R.shape[1]))

    xy_coords = np.stack([x_values, y_values])

    decomp_method = cascade.get_method(decomp_method)
    filter_method = cascade.get_method(bandpass_filter_method)
    if noise_method is not None:
        init_noise, generate_noise = noise.get_method(noise_method)

    # advect the previous precipitation fields to the same position with the
    # most recent one (i.e. transform them into the Lagrangian coordinates)
    R = R[-(ar_order + 1) :, :, :].copy()
    extrap_kwargs = extrap_kwargs.copy()
    extrap_kwargs["xy_coords"] = xy_coords
    res = []
    f = lambda R, i: extrapolator_method(
        R[i, :, :], V, ar_order - i, "min", **extrap_kwargs
    )[-1]
    for i in range(ar_order):
        if not dask_imported:
            R[i, :, :] = f(R, i)