a) Reader

The reader module (reader.py) contains functions to read the configuration file and generate file lists. The function wrainfo.reader.create_flist() loads FURUNO raw data for a defined period of time into a list. The file list is compared with an error file list. The error file list contains manually defined erroneous files and these are removed from the created file lists in order to exclude these data from the processing. The two individual file lists are merged into one. Per default, scn/scnx format is supported as input data format, optionally data can be provided in HDF5 format. A file list with all raw data of the defined period is returned.

Scn/scnx files are processed by the function wrainfo.reader.read_single_file(). A raw dataset is read and returned as an xarray dataset. Figure 2 shows the most important polarimetric variables of a FURUNO raw dataset, which are read by the function.

Raw reflectivity (ZH) provides information on the spatial distribution of a recorded precipitation event extending from north to south-east (Figure 2).

By means of the polarimetric variable RHOHV, non-meteorological (ground clutter) and meteorological echoes can be distinguished in the raw reflectivity (ZH). If the cross correlation coefficient (RHOHV) is greater than 0.9, it is assumed that mostly precipitation occurs [17]. Thus, Figure 2 shows that the individual pixels with high reflectivity, such as in the west and south, are ground clutter. KDP allows conclusions to be drawn about the size and concentration of raindrops. Thus, high positive values in KDP correlate with high reflectivity, which can be seen in the north (Figure 2).

The polarimetric variable PHIDP is a measure of the attenuation. Strong attenuation is visible in the north (Figure 2). Even if the attenuation is not strong, an attenuation correction must be made.

b) Clutter detection

The clutter detection module (clutter.py) contains functions to detect clutter and then remove it from the raw data. Furthermore, a clutter map can be created with this module.

A clutter map can be generated using the wrainfo.clutter.create_clutter_map_sequential() function. To do this, a file list must first be created that contains all raw data over a defined time period. This can be done using the wrainfo.reader.create_filelist() function.

The created file list is one of the function parameters. The parameter threshold specifies that only positive values of the uncorrected reflectivity, which are in the boolean dataset of the clutter map, are considered. If the status parameter is set to true, then the processing progress of the clutter map is obtained as a percentage. Only the reflectivity values are extracted from the raw data and added. Ground clutter is located at the points where high reflectivities occur in each dataset. Afterwards a boolean dataset is created where the added reflectivity values must be greater than the threshold. The clutter maps generated with WRaINfo contain a dataset with the summed reflectivities and another boolean dataset (1 = clutter, 0 = no clutter).

Figure 3 shows a clutter map for the elevation angle of 0.5°. Pixels which are identified as clutter, are marked in yellow. These are wind parks and high buildings in the town, where the weather radar is located.

For clutter detection the fuzzy echo classification from wradlib is used [18, 19]. The function wrainfo.clutter.fuzzy_echo_classification() needs as input a two dimensional (azimuth and range) raw dataset of the weather radar and the clutter map as xarray dataset. The raw dataset may additionally contain time as a dimension. Then the parameter “dims” must be adjusted accordingly.

For the fuzzy echo classification all polarimetric variables (RHOHV, PHIDP, ZDR, VRAD) as well as the generated clutter map (CMAP) are used as described in Vulpiani et al. [19].

The default values for the weighting of each variable correspond to the study of Crislogo et. al [18] and can be found in Table 3.

The new variable “FUZZ” is appended to the raw dataset. This is the clutter map resulting from the fuzzy echo classification and is needed for the function wrainfo.clutter.dbzh_no_clutter(). Figure 4 shows the clutter map and results of the fuzzy echo classification. All radar bins that have been classified as clutter are now marked in yellow.

To remove the pixels masked as ground clutter from the raw reflectivity, the wrainfo.clutter.dbzh_no_clutter() function is used. The output dataset of the fuzzy echo classification is needed for this function. In addition, the static clutter map has to be loaded. The ground clutter is excluded from the raw reflectivity using the clutter map and the previously created variable “FUZZ”. All pixels that have been identified as clutter are flagged as NaN values. Figure 5 shows the raw reflectivity and the clutter corrected reflectivity. The returned xarray dataset gets a new variable “DBZH_no_clutter”.

c) Attenuation correction

After the clutter detection is done, the attenuation correction follows. This is accomplished by the wrainfo.attenuation_corr.phase_zphi() function and processes the phase PHIDP. The function needs the data array of PHIDP from the clutter corrected xarray dataset. The parameters “start_range” and “range” of the moving window have to be defined. Within the function a binary array (True = rain, False = not-rain) rolling a range-long sum from PHIDP is generated. Afterwards the first occurrence of the maximum (from front and back) of PHIDP binary array is searched for. This gives the indices of the centres of the rolling window. Then min, max, median and mean of these are calculated. Finally, the true start and end indices of the phase PHIDP are calculated (+/- half window length). The returned dataset additionally contains the first and last values (min, max, median, mean) of the phase PHIDP (Figure 6).

In the Z-PHI method, the rain profile along the radar beam is derived from the reflectivity profile (ZH). The required boundary condition is the total attenuation of the radar beam, which is derived from PHIDP.

In the following function wrainfo.attenuation_corr.zphi() the specific attenuation (AH) is calculated using an optimized algorithm based on the ZPHI-Method [7, 8]. The function needs the clutter corrected reflectivity and the added dataset of the function wrainfo.attenuation_corr.phase_zphi() as function parameters. To define the boundary conditions, within the function wrainfo.attenuation_corr.zphi() first $\Delta {\Phi }_{DP}$M7 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[\Delta {\Phi _{DP}}\] \end{document} is calculated (eq. 1) from the parameters “last_median” and “first_median” (Figure 7). In eq. 1 r₁ and r₂ represent the integration bounds.

These are contained in the returned xarray dataset of the function wrainfo.attenuation_corr.phase_zphi(). First any $\Delta {\Phi }_{DP}$M8 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[\Delta {\Phi _{\rm{DP}}}\] \end{document} values less than 0 are set to 0. Then the function of $\Delta {\Phi }_{DP}$M9 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[\Delta {\Phi _{\rm{DP}}}\] \end{document} ($\Delta {\Phi }_{DP}$M10 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[{\rm{f\Delta }} {{\rm{\Phi }}_{{\rm{DP}}}}\] \end{document} ) is calculated (eq. 2). According to eq. 3 we prepare attenuated reflectivity to the power of $b={[{\text{Z}}_{\text{a}}(\text{r})]}^{\text{b}}$M11 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[b = {[{{\rm{Z}}_{\rm{a}}}({\rm{r}})]^{\rm{b}}}\] \end{document} , where b is a constant parameter (usually within 0.6–0.9 at microwave frequencies) and the integral of that from each range bin to maximum range [9]. The calculation of that integral is done per ray by cumulative trapezoidal integration over the range dimension and following subtraction from its maximum value. α is calculated as $Z = \text{a} \cdot R^{\text{b}}$M12 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[\alpha {\rm{ }} = \frac{{\rm{A}}}{{{\rm{KDP}}}}\] \end{document} and is assumed as constant or temperature dependent [9, 20]. A stands for the specific attenuation.

For performance reasons eq. 3 is rearranged as

The integral (r₁, r₂) is extracted from the above integral by selection using the first index provided from wrainfo.attenuation_corr.phase_zphi() function. In advance KDP is derived from AH and further PHIDP is recalculated from KDP.

The ZPHI-based AH, KDP, and PHIDP are merged into one xarray dataset which is finally returned.

Figure 7 shows the recalculated polarimetric variables PHIDP and KDP, respectively, and the respective raw dataset as well as the calculated specific attenuation (AH).

Subsequently, the attenuation correction of the reflectivity is performed with the recalculated PHIDP (eq. 4). Figure 8 shows a comparison of the reflectivity before and after attenuation correction. α depends on temperature and drop size distribution. The range of values for X-Band weather radars is between 0.14 – 0.35 dB/deg. α is mostly assumed to be constant (α = 0.28) and corresponds to the climatic mean. Certain uncertainties, such as variability in summer and winter, cannot be excluded. Thus, 0.28 is used for α in eq. 5.

d) Precipitation estimation

The precipitation module (precipitation.py) is used to estimate precipitation rates from clutter and attenuation corrected reflectivity based on the Z-R relationship. This is an empirical relationship between the reflectivity (Z) of the weather radar and the precipitation rate (R). In equation 6 a and b are the relationship parameter and depend on the rainfall type (stratiform or convective). Marshall and Palmer [21] defines default values for a and b are: a = 200 and b = 1.6.

The function wrainfo.precipitation.qpe_zr() is used for the precipitation estimation. The function needs as input parameter the returned dataset of the attenuation correction and the name of the data array, which contains the attenuation corrected reflectivity. The directory to the configuration file also has to be defined as a function parameter. In the configuration file the scan interval in seconds and the dimensions of the dataset have to be defined.

First, the precipitation intensity [mm/h] is derived from the corrected reflectivity in this function. Second, the amount of precipitation [mm] within the specified scan interval is calculated based on the precipitation intensity [mm/h]. Finally, the returned xarray dataset contains a new variable “PREC_ZR”. Figure 9 shows the estimated precipitation amount for the clutter and attenuation corrected reflectivity.

e) Georeferencing and gridding

The geometry module is used to project the polar data into a Cartesian grid and to write georeferenced radar data as netCDF files to the hard disk.

The function wrainfo.geometry.furuno_georeferencing() uses as input parameter an xarray dataset and the directory to the configuration file, because the target coordinate reference system, the spatial resolution in pixels and the moments, which should be georeferenced have to be defined in the configuration file. The function returns a georeference xarray dataset, which shows in Figure 10.

The writing of (georeferenced) radar data as netCDF files to the hard disk is done with the function wrainfo.geometry.furuno_sweep_to_netcdf(). This function has the same parameters as the function for georeferencing. Additionally, the data type can be defined. This refers to the level of data processing. The output directory has to be specified in the configuration file.

f) Error and compression module

The error file list (error_flist.py) and compression (compression.py) module allow the data management of FURUNO data on the hard disk. With the former it is possible to exclude erroneous files from the data processing. With the compression module it is possible to store raw data per day in a compressed tar.gz file.

g) Process chains

WRaINfo contains the module process_chains.py, which is based on the modules reader.py, clutter.py, attenuation_corr.py, precipitation.py and geometry.py (Figure 1).

The module process_chains.py shows how the operational mode can be set up to process the raw data of the FURUNO weather radars and subsequently generate georeferenced precipitation products with any temporal resolution. In this paragraph, the individual process chains are described in detail which allow an operational workflow.