How to use the pysteps.utils function in pysteps

ds = stp.rcparams.data_sources[data_source]

## find radar field filenames
input_files = stp.io.find_by_date(startdate, ds.root_path, ds.path_fmt, ds.fn_pattern, 
                                  ds.fn_ext, ds.timestep, n_prvs_times, 0)

## read radar field files
importer = stp.io.get_method(ds.importer, type="importer")
R, _, metadata = stp.io.read_timeseries(input_files, importer, **ds.importer_kwargs)
Rmask = np.isnan(R)

# Prepare input files
print("Prepare the data...")

## if necessary, convert to rain rates [mm/h]    
converter = stp.utils.get_method("mm/h")
R, metadata = converter(R, metadata)

## threshold the data
R[R

Lucas, B. D. and Kanade, T.: An iterative image registration technique with
    an application to stereo vision, in: Proceedings of the 1981 DARPA Imaging
    Understanding Workshop, pp. 121–130, 1981.
    """

    input_images = input_images.copy()

    if verbose:
        print("Computing the motion field with the Lucas-Kanade method.")
        t0 = time.time()

    nr_fields = input_images.shape[0]
    domain_size = (input_images.shape[1], input_images.shape[2])

    feature_detection_method = utils.get_method(fd_method)
    interpolation_method = utils.get_method(interp_method)

    if fd_kwargs is None:
        fd_kwargs = dict()

    if lk_kwargs is None:
        lk_kwargs = dict()

    if interp_kwargs is None:
        interp_kwargs = dict()

    xy = np.empty(shape=(0, 2))
    uv = np.empty(shape=(0, 2))
    for n in range(nr_fields - 1):

        # extract consecutive images
        prvs_img = input_images[n, :, :].copy()

print("probability matching:   %s" % probmatching_method)
    print("FFT method:             %s" % fft_method)
    print("")

    print("Parameters:")
    print("-----------")
    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:

print("Prepare the data...")

## if necessary, convert to rain rates [mm/h]    
converter = stp.utils.get_method(unit)
R, metadata = converter(R, metadata)

## threshold the data
R[R

ds = []
    for fname, tstep in zip(*fns):
        ds.append(load_mch(fname, tstep, **importer_kwargs))

    ds = xr.concat(ds, dim="time")
    return ds


ds = create_timeseries(fns, **importer_kwargs)
print(ds)

###############################################################################
# Convert and transform the data
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

ds["RR"] = sp.utils.to_rainrate(ds.precipitation)
ds["dBZ"] = sp.utils.dB_transform(ds.precipitation)
print(ds)

# Read the reference radar composite
importer = io.get_method(importer_name, "importer")
reference_field, quality, metadata = io.read_timeseries(fns, importer,
                                                        **importer_kwargs)

del quality  # Not used

reference_field = np.squeeze(reference_field)  # Remove time dimension

###############################################################################
# Preprocess the data
# ~~~~~~~~~~~~~~~~~~~

# Convert to mm/h
reference_field, metadata = stp.utils.to_rainrate(reference_field, metadata)

# Mask invalid values
reference_field = np.ma.masked_invalid(reference_field)

# Plot the reference precipitation
plot_precip_field(reference_field, title="Reference field")
plt.show()

# Log-transform the data [dBR]
reference_field, metadata = stp.utils.dB_transform(reference_field,
                                                   metadata,
                                                   threshold=0.1,
                                                   zerovalue=-15.0)

print("Precip. pattern shape: " + str(reference_field.shape))

print("Verifying the nowcast (%02d) ..." % countnwc)
        
        # Read observations
        
        ## find radar field filenames
        input_files = stp.io.find_by_date(startdate, ds.root_path, ds.path_fmt, ds.fn_pattern,
                                          ds.fn_ext, ds.timestep, 0, p["n_lead_times"])
                                          
        ## read radar field files
        importer = stp.io.get_method(ds.importer, type="importer")
        R_obs, _, metadata_obs = stp.io.read_timeseries(input_files, importer, **ds.importer_kwargs)
        R_obs = R_obs[1:,:,:]
        metadata_obs["timestamps"] = metadata_obs["timestamps"][1:]
        
        ## if necessary, convert to rain rates [mm/h]   
        converter = stp.utils.get_method("mm/h")        
        R_obs, metadata_obs = converter(R_obs, metadata_obs)  
        
        ## threshold the data
        R_obs[R_obs < p["r_threshold"]] = 0.0
        metadata_obs["threshold"] = p["r_threshold"]
            
        # Load the nowcast
        
        ## filename of the nowcast netcdf
        infn = os.path.join(path_to_nwc, "%s_nowcast.netcdf" % startdate.strftime("%Y%m%d%H%M"))
        
        print("     read: %s" % infn)
            
        ## read netcdf
        R_fct, metadata_fct = stp.io.import_netcdf_pysteps(infn)
        timestamps = metadata_fct["timestamps"]

ds.fn_ext, ds.timestep, n_prvs_times, 0)

## read radar field files
importer = stp.io.get_method(ds.importer, "importer")
R, _, metadata = stp.io.read_timeseries(input_files, importer, **ds.importer_kwargs)
Rmask = np.isnan(R)

# Prepare input files
print("Prepare the data...")

## if requested, make sure we work with a square domain
reshaper = stp.utils.get_method(adjust_domain)
R, metadata = reshaper(R, metadata, method="pad")

## if necessary, convert to rain rates [mm/h]
converter = stp.utils.get_method("mm/h")
R, metadata = converter(R, metadata)

## threshold the data
R[R

compact_output : bool
        Applicable if output_domain is "spectral". If set to True, only the
        parts of the Fourier spectrum with nonzero filter weights are stored.
        Defaults to False.

    Returns
    -------
    out : ndarray
        A dictionary described in the module documentation.
        The number of cascade levels is determined from the filter
        (see :py:mod:`pysteps.cascade.bandpass_filters`).

    """
    fft = kwargs.get("fft_method", "numpy")
    if isinstance(fft, str):
        fft = utils.get_method(fft, shape=field.shape)
    normalize = kwargs.get("normalize", False)
    mask = kwargs.get("MASK", None)
    input_domain = kwargs.get("input_domain", "spatial")
    output_domain = kwargs.get("output_domain", "spatial")
    compute_stats = kwargs.get("compute_stats", False)
    compact_output = kwargs.get("compact_output", False)

    if normalize and not compute_stats:
        raise ValueError("incorrect input arguments: normalize=True but compute_stats=False")

    if len(field.shape) != 2:
        raise ValueError("The input is not two-dimensional array")

    if mask is not None and mask.shape != field.shape:
        raise ValueError("Dimension mismatch between field and MASK:"
                         + "field.shape=" + str(field.shape)

if output_domain == "spatial" or (compute_stats and mask is not None):
            field__ = fft.irfft2(field_)
        else:
            field__ = field_

        if compute_stats:
            if output_domain == "spatial" or (compute_stats and mask is not None):
                if mask is not None:
                    masked_field = field__[mask]
                else:
                    masked_field = field__
                mean = np.mean(masked_field)
                std = np.std(masked_field)
            else:
                mean = utils.spectral.mean(field_, bp_filter["shape"])
                std = utils.spectral.std(field_, bp_filter["shape"])
            
            means.append(mean)
            stds.append(std)

        if output_domain == "spatial":
            field_ = field__
        if normalize:
            field_ = (field_ - mean) / std
        if output_domain == "spectral" and compact_output:
            weight_mask = bp_filter["weights_2d"][k, :, :] > 1e-12
            field_ = field_[weight_mask]
            weight_masks.append(weight_mask)
        field_decomp.append(field_)

    result["domain"] = output_domain