def decompose_NWP(
    R_NWP,
    NWP_model,
    analysis_time,
    timestep,
    valid_times,
    output_path,
    num_cascade_levels=8,
    num_workers=1,
    decomp_method="fft",
    fft_method="numpy",
    domain="spatial",
    normalize=True,
    compute_stats=True,
    compact_output=True,
):
    """Decomposes the NWP forecast data into cascades and saves it in
    a netCDF file

    Parameters
    ----------
    R_NWP: array-like
      Array of dimension (n_timesteps, x, y) containing the precipitation forecast
      from some NWP model.
    NWP_model: str
      The name of the NWP model
    analysis_time: numpy.datetime64
      The analysis time of the NWP forecast. The analysis time is assumed to be a
      numpy.datetime64 type as imported by the pysteps importer
    timestep: int
      Timestep in minutes between subsequent NWP forecast fields
    valid_times: array_like
      Array containing the valid times of the NWP forecast fields. The times are
      assumed to be numpy.datetime64 types as imported by the pysteps importer.
    output_path: str
      The location where to save the file with the NWP cascade. Defaults to the
      path_workdir specified in the rcparams file.
    num_cascade_levels: int, optional
      The number of frequency bands to use. Must be greater than 2. Defaults to 8.
    num_workers: int, optional
      The number of workers to use for parallel computation. Applicable if dask
      is enabled or pyFFTW is used for computing the FFT. When num_workers>1, it
      is advisable to disable OpenMP by setting the environment variable
      OMP_NUM_THREADS to 1. This avoids slowdown caused by too many simultaneous
      threads.

    Other Parameters
    ----------------
    decomp_method: str, optional
      A string defining the decomposition method to use. Defaults to "fft".
    fft_method: str or tuple, optional
      A string or a (function,kwargs) tuple defining the FFT method to use
      (see :py:func:`pysteps.utils.interface.get_method`).
      Defaults to "numpy". This option is not used if input_domain and
      output_domain are both set to "spectral".
    domain: {"spatial", "spectral"}, optional
      If "spatial", the output cascade levels are transformed back to the
      spatial domain by using the inverse FFT. If "spectral", the cascade is
      kept in the spectral domain. Defaults to "spatial".
    normalize: bool, optional
      If True, normalize the cascade levels to zero mean and unit variance.
      Requires that compute_stats is True. Implies that compute_stats is True.
      Defaults to False.
    compute_stats: bool, optional
      If True, the output dictionary contains the keys "means" and "stds"
      for the mean and standard deviation of each output cascade level.
      Defaults to False.
    compact_output: bool, optional
      Applicable if output_domain is "spectral". If set to True, only the
      parts of the Fourier spectrum with non-negligible filter weights are
      stored. Defaults to False.


    Returns
    -------
    None
    """

    if not NETCDF4_IMPORTED:
        raise MissingOptionalDependency(
            "netCDF4 package is required to save the decomposed NWP data, "
            "but it is not installed"
        )

    # Make a NetCDF file
    output_date = f"{analysis_time.astype('datetime64[us]').astype(datetime.datetime):%Y%m%d%H%M%S}"
    outfn = Path(output_path) / f"cascade_{NWP_model}_{output_date}.nc"
    ncf = netCDF4.Dataset(outfn, "w", format="NETCDF4")

    # Express times relative to the zero time
    zero_time = np.datetime64("1970-01-01T00:00:00", "ns")
    valid_times = np.array(valid_times) - zero_time
    analysis_time = analysis_time - zero_time

    # Set attributes of decomposition method
    ncf.domain = domain
    ncf.normalized = int(normalize)
    ncf.compact_output = int(compact_output)
    ncf.analysis_time = int(analysis_time)
    ncf.timestep = int(timestep)

    # Create dimensions
    ncf.createDimension("time", R_NWP.shape[0])
    ncf.createDimension("cascade_levels", num_cascade_levels)
    ncf.createDimension("x", R_NWP.shape[2])
    ncf.createDimension("y", R_NWP.shape[1])

    # Create variables (decomposed cascade, means and standard deviations)
    R_d = ncf.createVariable(
        "pr_decomposed",
        np.float32,
        ("time", "cascade_levels", "y", "x"),
        zlib=True,
        complevel=4,
    )
    means = ncf.createVariable("means", np.float64, ("time", "cascade_levels"))
    stds = ncf.createVariable("stds", np.float64, ("time", "cascade_levels"))
    v_times = ncf.createVariable("valid_times", np.float64, ("time",))
    v_times.units = "nanoseconds since 1970-01-01 00:00:00"

    # The valid times are saved as an array of floats, because netCDF files can't handle datetime types
    v_times[:] = np.array([np.float64(valid_times[i]) for i in range(len(valid_times))])

    # Decompose the NWP data
    filter_g = filter_gaussian(R_NWP.shape[1:], num_cascade_levels)
    fft = utils_get_method(fft_method, shape=R_NWP.shape[1:], n_threads=num_workers)
    decomp_method, _ = cascade_get_method(decomp_method)

    for i in range(R_NWP.shape[0]):
        R_ = decomp_method(
            field=R_NWP[i, :, :],
            bp_filter=filter_g,
            fft_method=fft,
            input_domain=domain,
            output_domain=domain,
            normalize=normalize,
            compute_stats=compute_stats,
            compact_output=compact_output,
        )

        # Save data to netCDF file
        # print(R_["cascade_levels"])
        R_d[i, :, :, :] = R_["cascade_levels"]
        means[i, :] = R_["means"]
        stds[i, :] = R_["stds"]

    # Close the file
    ncf.close()
Decomposes the NWP forecast data into cascades and saves it in
a netCDF file

Parameters
----------
R_NWP: array-like
Array of dimension (n_timesteps, x, y) containing the precipitation forecast
from some NWP model.
NWP_model: str
The name of the NWP model
analysis_time: numpy.datetime64
The analysis time of the NWP forecast. The analysis time is assumed to be a
numpy.datetime64 type as imported by the pysteps importer
timestep: int
Timestep in minutes between subsequent NWP forecast fields
valid_times: array_like
Array containing the valid times of the NWP forecast fields. The times are
assumed to be numpy.datetime64 types as imported by the pysteps importer.
output_path: str
The location where to save the file with the NWP cascade. Defaults to the
path_workdir specified in the rcparams file.
num_cascade_levels: int, optional
The number of frequency bands to use. Must be greater than 2. Defaults to 8.
num_workers: int, optional
The number of workers to use for parallel computation. Applicable if dask
is enabled or pyFFTW is used for computing the FFT. When num_workers>1, it
is advisable to disable OpenMP by setting the environment variable
OMP_NUM_THREADS to 1. This avoids slowdown caused by too many simultaneous
threads.

Other Parameters
----------------
decomp_method: str, optional
A string defining the decomposition method to use. Defaults to "fft".
fft_method: str or tuple, optional
A string or a (function,kwargs) tuple defining the FFT method to use
(see :py:func:`pysteps.utils.interface.get_method`).
Defaults to "numpy". This option is not used if input_domain and
output_domain are both set to "spectral".
domain: {"spatial", "spectral"}, optional
If "spatial", the output cascade levels are transformed back to the
spatial domain by using the inverse FFT. If "spectral", the cascade is
kept in the spectral domain. Defaults to "spatial".
normalize: bool, optional
If True, normalize the cascade levels to zero mean and unit variance.
Requires that compute_stats is True. Implies that compute_stats is True.
Defaults to False.
compute_stats: bool, optional
If True, the output dictionary contains the keys "means" and "stds"
for the mean and standard deviation of each output cascade level.
Defaults to False.
compact_output: bool, optional
Applicable if output_domain is "spectral". If set to True, only the
parts of the Fourier spectrum with non-negligible filter weights are
stored. Defaults to False.


Returns
-------
None
pysteps/blending/utils.py