How to use the pysteps.io function in pysteps

# Load data source config
root_path = rcparams.data_sources[data_source]["root_path"]
path_fmt = rcparams.data_sources[data_source]["path_fmt"]
fn_pattern = rcparams.data_sources[data_source]["fn_pattern"]
fn_ext = rcparams.data_sources[data_source]["fn_ext"]
importer_name = rcparams.data_sources[data_source]["importer"]
importer_kwargs = rcparams.data_sources[data_source]["importer_kwargs"]
timestep = rcparams.data_sources[data_source]["timestep"]

# Find the radar files in the archive
fns = io.find_by_date(
    date, root_path, path_fmt, fn_pattern, fn_ext, timestep, num_prev_files=2
)

# 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

#
# First, we will show how a netcdf file is currenty imported in the pysteps
# workflow and how this can be done using xarray.

# The path to an example BOM netcdf file
root_path = sp.rcparams.data_sources["bom"]["root_path"]
filename = os.path.join(
    root_path, "prcp-cscn/2/2018/06/16/2_20180616_150000.prcp-cscn.nc"
)

###############################################################################
# pysteps original
# ~~~~~~~~~~~~~~~~

# Read data
R, __, metadata = sp.io.import_bom_rf3(filename)
metadata

# Visualization
sp.visualization.plot_precip_field(
    R,
    map=None,# "cartopy",
    type="depth",
    units="mm",
    geodata=metadata,
    drawlonlatlines=True,
)

###############################################################################
# pysteps using xarray
# ~~~~~~~~~~~~~~~~~~~~

from pysteps.utils import conversion, rapsd, transformation
from pysteps.visualization import plot_precip_field, plot_spectrum1d

###############################################################################
# Read precipitation field
# ------------------------
#
# First thing,  the radar composite is imported and transformed in units
# of dB.
# This image will be used to train the Fourier filters that are necessary to
# produce the fields of spatially correlated noise.

# 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"]

## parameters
r_threshold         = 0.1 # [mm/h]
num_cascade_levels  = 6
unit                = "mm/h"    # mm/h or dBZ
transformation      = "dB"      # None or dB 
adjust_domain       = None      # None or "square"

# Read-in the data
print('Read the data...')
startdate  = datetime.datetime.strptime(startdate_str, "%Y%m%d%H%M")

## import data specifications
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, 0, 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)
R = R.squeeze() # since this contains just one frame
Rmask = np.isnan(R)

# Prepare input files
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