How to use the pysteps.utils.transformation function in pysteps

# Import the example radar composite
root_path = rcparams.data_sources["mch"]["root_path"]
filename = os.path.join(root_path, "20160711", "AQC161932100V_00005.801.gif")
R, _, metadata = io.import_mch_gif(filename, product="AQC", unit="mm", accutime=5.0)

# Convert to mm/h
R, metadata = conversion.to_rainrate(R, metadata)

# Nicely print the metadata
pprint(metadata)

# Plot the rainfall field
plot_precip_field(R, geodata=metadata)

# Log-transform the data
R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)

# Assign the fill value to all the Nans
R[~np.isfinite(R)] = metadata["zerovalue"]

###############################################################################
# Parametric filter
# -----------------
#
# In the parametric approach, a power-law model is used to approximate the power
# spectral density (PSD) of a given rainfall field.
#
# The parametric model uses  a  piece-wise  linear  function  with  two  spectral
# slopes (beta1 and beta2) and one breaking point

# Fit the parametric PSD to the observation
Fp = initialize_param_2d_fft_filter(R)

n_leadtimes,
    n_ens_members,
    n_cascade_levels=6,
    R_thr=-10.0,
    kmperpixel=2.0,
    timestep=timestep,
    decomp_method="fft",
    bandpass_filter_method="gaussian",
    noise_method="nonparametric",
    vel_pert_method="bps",
    mask_method="incremental",
    seed=seed,
)

# Back-transform to rain rates
R_f = transformation.dB_transform(R_f, threshold=-10.0, inverse=True)[0]

# Plot some of the realizations
fig = figure()
for i in range(4):
    ax = fig.add_subplot(221 + i)
    ax.set_title("Member %02d" % i)
    plot_precip_field(R_f[i, -1, :, :], geodata=metadata, colorbar=False, axis="off")
tight_layout()

###############################################################################
# Verification
# ------------
#
# Pysteps includes a number of verification metrics to help users to analyze 
# the general characteristics of the nowcasts in terms of consistency and 
# quality (or goodness). 

# Read the data from the archive
importer = io.get_method(importer_name, "importer")
R, _, metadata = io.read_timeseries(fns, importer, **importer_kwargs)

# Convert to rain rate
R, metadata = conversion.to_rainrate(R, metadata)

# Upscale data to 2 km to limit memory usage
R, metadata = dimension.aggregate_fields_space(R, metadata, 2000)

# Plot the rainfall field
plot_precip_field(R[-1, :, :], geodata=metadata)

# Log-transform the data to unit of dBR, set the threshold to 0.1 mm/h,
# set the fill value to -15 dBR
R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)

# Set missing values with the fill value
R[~np.isfinite(R)] = -15.0

# Nicely print the metadata
pprint(metadata)

###############################################################################
# Deterministic nowcast with S-PROG
# ---------------------------------
#
# First, the motiong field is estimated using a local tracking approach based
# on the Lucas-Kanade optical flow.
# The motion field can then be used to generate a deterministic nowcast with
# the S-PROG model, which implements a scale filtering appraoch in order to
# progressively remove the unpredictable spatial scales during the forecast.

-------
    R : array-like
        Array of any shape containing the converted units.
    metadata : dict
        The metadata with updated attributes.

    """

    R = R.copy()
    metadata = metadata.copy()

    if metadata["transform"] is not None:

        if metadata["transform"] == "dB":

            R, metadata = transformation.dB_transform(R, metadata,
                                                      inverse=True)

        elif metadata["transform"] in ["BoxCox", "log"]:

            R, metadata = transformation.boxcox_transform(R, metadata,
                                                          inverse=True)

        elif metadata["transform"] == "NQT":

            R, metadata = transformation.NQ_transform(R, metadata,
                                                      inverse=True)

        elif metadata["transform"] == "sqrt":

            R, metadata = transformation.sqrt_transform(R, metadata,
                                                        inverse=True)

del quality  # Not used

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

# Convert to mm/h
R, metadata = conversion.to_rainrate(R, metadata)

# Keep the reference frame in mm/h and its mask (for plotting purposes)
ref_mm = R[0, :, :].copy()
mask = np.ones(ref_mm.shape)
mask[~np.isnan(ref_mm)] = np.nan

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

# Keep the reference frame in dBR (for plotting purposes)
ref_dbr = R[0].copy()
ref_dbr[ref_dbr < -10] = np.nan

# Plot the reference field
plot_precip_field(ref_mm, title="Reference field")
circle = plt.Circle((620, 400), 100, color="b", clip_on=False, fill=False)
plt.gca().add_artist(circle)
plt.show()

###############################################################################
# Notice the "half-in, half-out" precipitation area within the blue circle.
# As we are going to show next, the tracking algorithm can erroneously interpret

# Read the radar composites
importer = io.get_method(importer_name, "importer")
Z, _, metadata = io.read_timeseries(fns, importer, **importer_kwargs)

# Convert to rain rate using the finnish Z-R relationship
R, metadata = conversion.to_rainrate(Z, metadata, 223.0, 1.53)

# Plot the rainfall field
plot_precip_field(R[-1, :, :], geodata=metadata)

# Store the last frame for plotting it later later
R_ = R[-1, :, :].copy()

# Log-transform the data to unit of dBR, set the threshold to 0.1 mm/h,
# set the fill value to -15 dBR
R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)

# Nicely print the metadata
pprint(metadata)

###############################################################################
# Compute the nowcast
# -------------------
#
# The extrapolation nowcast is based on the estimation of the motion field,
# which is here performed using a local tacking approach (Lucas-Kanade).
# The most recent radar rainfall field is then simply advected along this motion
# field in oder to produce an extrapolation forecast.

# Estimate the motion field with Lucas-Kanade
oflow_method = motion.get_method("LK")
V = oflow_method(R[-3:, :, :])

if metadata["transform"] is not None:

        if metadata["transform"] == "dB":

            R, metadata = transformation.dB_transform(R, metadata,
                                                      inverse=True)

        elif metadata["transform"] in ["BoxCox", "log"]:

            R, metadata = transformation.boxcox_transform(R, metadata,
                                                          inverse=True)

        elif metadata["transform"] == "NQT":

            R, metadata = transformation.NQ_transform(R, metadata,
                                                      inverse=True)

        elif metadata["transform"] == "sqrt":

            R, metadata = transformation.sqrt_transform(R, metadata,
                                                        inverse=True)

        else:

            raise ValueError("Unknown transformation %s"
                             % metadata["transform"])

    if metadata["unit"] == "mm/h":

        pass

filename = os.path.join(
    root_path, "20160928", "201609281600_fmi.radar.composite.lowest_FIN_SUOMI1.pgm.gz"
)
R, _, metadata = io.import_fmi_pgm(filename, gzipped=True)

# Convert to rain rate
R, metadata = conversion.to_rainrate(R, metadata)

# Nicely print the metadata
pprint(metadata)

# Plot the rainfall field
plot_precip_field(R, geodata=metadata)

# Log-transform the data
R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)

###############################################################################
# 2D Fourier spectrum
# --------------------
#
# Compute and plot the 2D Fourier power spectrum of the precipitaton field.

# Set Nans as the fill value
R[~np.isfinite(R)] = metadata["zerovalue"]

# Compute the Fourier transform of the input field
F = abs(np.fft.fftshift(np.fft.fft2(R)))

# Plot the power spectrum
M, N = F.shape
fig, ax = pyplot.subplots()

#
# To find a suitable lambda, we will experiment with a range of values
# and select the one that produces the most symmetric distribution, i.e., the
# lambda associated with a value of skewness closest to zero.
# To visually compare the results, the transformed data are standardized.
#
# .. _`Box and Cox (1964)`:https://doi.org/10.1111/j.2517-6161.1964.tb00553.x

data = []
labels = []
skw = []

# Test a range of values for the transformation parameter Lambda
Lambdas = np.linspace(-0.4, 0.4, 11)
for i, Lambda in enumerate(Lambdas):
    R_, _ = transformation.boxcox_transform(R, metadata, Lambda)
    R_ = (R_ - np.mean(R_)) / np.std(R_)
    data.append(R_)
    labels.append("{0:.2f}".format(Lambda))
    skw.append(skew(R_))  # skewness

# Plot the transformed data distribution as a function of lambda
plot_distribution(data, labels, skw)
plt.title("Box-Cox transform")
plt.tight_layout()
plt.show()

# Best lambda
idx_best = np.argmin(np.abs(skw))
Lambda = Lambdas[idx_best]

print("Best parameter lambda: %.2f\n(skewness = %.2f)" % (Lambda, skw[idx_best]))

methods_objects["square"] = dimension.square_domain
    methods_objects["upscale"] = dimension.aggregate_fields_space

    # image processing methods
    methods_objects["shitomasi"] = images.ShiTomasi_detection
    methods_objects["morph_opening"] = images.morph_opening

    # interpolation methods
    methods_objects["rbfinterp2d"] = interpolate.rbfinterp2d

    # spectral methods
    methods_objects["rapsd"] = spectral.rapsd
    methods_objects["rm_rdisc"] = spectral.remove_rain_norain_discontinuity

    # transformation methods
    methods_objects["boxcox"] = transformation.boxcox_transform
    methods_objects["box-cox"] = transformation.boxcox_transform
    methods_objects["db"] = transformation.dB_transform
    methods_objects["decibel"] = transformation.dB_transform
    methods_objects["log"] = transformation.boxcox_transform
    methods_objects["nqt"] = transformation.NQ_transform
    methods_objects["sqrt"] = transformation.sqrt_transform

    # FFT methods
    if name in ["numpy", "pyfftw", "scipy"]:
        if "shape" not in kwargs.keys():
            raise KeyError("mandatory keyword argument shape not given")
        return _get_fft_method(name, **kwargs)
    else:
        try:
            return methods_objects[name]
        except KeyError as e: