Open-Source Data Processing Chain for Marche Region X-band Weather Radar

: In the framework of three projects co-funded by the European Union (EU), an X-band (8-12 GHz frequency range) polarimetric radar was installed in the Marche Region territory (East-Central Italy) at Cingoli municipality, province of Macerata. The radar site is located at about 750 m above sea level and about 30 km away from the Adriatic Sea. The radar, managed by the Marche Region Civil Protection Service, is employed for weather monitoring purposes and is in pre-operational stage. It is known that radar measurements are affected by various sources of error, to be addressed in order to improve the accuracy of final products. Among these, the most important are radar calibration, ground and sea clutter, beam blockage, rain attenuation, wet-radome attenuation, beam-broadening, non-uniform beam filling, vertical variability of precipitation and wireless local-area-network (WLAN) interferences. Nowadays quantitative rainfall estimation using X-band weather radar are essential to meet requirements for flood forecasting, water management and many hydro-meteorological applications. Besides higher resolution, X-band radars are cost-effective compared to S-or C-band radars because of smaller antenna size. On the other hand, main disadvantages of such systems are the large influence of attenuation by liquid water and a relatively short range. In this work, we will present the data-processing chain developed ad hoc in order to remove or at least reduce the sources of error affecting Cingoli radar data. The performance of the data-processing chain was evaluated in the light of case studies related to meteorological events that interested Marche region territory in the last two years. The software was developed using open source technology. The current version of the chain does not take into account the echoes from sea clutter and the attenuation due to wet radome that can be significant at X-band; such issues will be addressed in a future work.


Introduction
Radars at different frequencies have been used over many decades for weather purposes. In particular, conventional S-and C-band systems are employed for long-range coverage observations while dualpolarization X-band radars are useful for networked applications [1]. In recent years, because of smaller antenna size, low-cost X-band radars have been employed to observe local weather phenomena and used for flood forecasting [2][3][4]. Compared to large and expensive S-and C-band systems, X-band radars cover a smaller area but with a higher spatio-temporal resolution [5] and are strongly affected by attenuation due to atmospheric liquid water [6]. In the last decade, a number of initiatives towards a more open weather radar science have been proposed [8]; some relevant examples are BALTRAD (written in a mixture of C/C++, Python and Java) [9] and Pythonbased libraries wradlib [10] and Py-ART [11]. For the present work we chose free high-level programming language Python since it is easier to integrate in an existing processing chain than BALTRAD and free wradlib library in particular as it is mostly data agnostic [7], i.e., it can process data in multiple formats or from multiple sources. The software was developed as a "wrapper" of wradlib. The effectiveness of the proposed dataprocessing chain was evaluated in the light of case studies related to meteorological events that recently interested Marche Region territory. Figure 1 shows the territory of Marche Region highlighted in red (on the left) and the position of the radar site (on the right). Cingoli radar covered by radome and its placement on the tower at about 20 m above the ground are shown in Figure 2. While in Figure 3, a picture of the radar on the tower without its radome is reported. Main technical characteristics of Cingoli dual polarization X-band radar are listed in Table 1.    Parameters of the operational volumetric scan for the Cingoli radar are reported in Table 2 and in Figure 4 a representation of the radar beam for 0°elevation and 50°azimuth is shown together with the beam-blocking fraction; Figure 4 was created using a Digital Elevation Model (DEM) of Marche region territory with resolution 100 m. We chose polar coordinates centered on the radar only for Figure 4 in order to highlight the distances relative to Cingoli site and the maximum radar range (120 km) in particular.  The definition of the operational scan has undergone changes over time and the choice of parameter values was made taking into account the fact that Cingoli radar is used by the Civil Protection Service for meteorological monitoring purposes. The first operation carried out for this work was converting volumetric radar data from the proprietary Eldes binary format (called MSM) into HDF5 (Hierarchical Data Format version 5) files which conform to the ODIM (OPERA Data Information Model) standard (http://www.eumetnet.eu/opera). ODIM is an information model for use with weather radar data and products [12]. All the filtering operations described here have been performed on volumetric radar data in polar format. Filtering the data in polar format has made it possible to considerably reduce the calculation times compared to processing the interpolated data. Radar data is interpolated on a regular Cartesian grid only at the end of the whole filtering chain in order to easily display standard products such as Vertical Maximum Intensity (VMI), Constant Altitude Plan Position Indicator (CAPPI), Surface Rainfall Intensity (SRI) and Surface Rainfall Total (SRT).

Cingoli Radar Data Processing Chain
This section details Cingoli radar data processing chain; in order to show the performance of the algorithms used, radar data is presented before and after filtering operations. As an example, we report data measured by the Cingoli radar on 2021-7-19 at 1:30 UTC at the lowest elevation in presence of precipitation; in particular, Figure 5 shows the unfiltered horizontal reflectivity (UZH) on the left and horizontal reflectivity Doppler-filtered by Eldes (ZH) on the right. Note that in July 2019 Cingoli radar scan consisted of seven elevations whose lowest angle was set at 1.5°. In addition to the weather feature (labeled as A in Figure 5), UZH data clearly shows the presence of ground clutter mostly due to Apennine chain (labeled as B), radio interferences (labeled as C), RACON (RAdar beaCON) echoes (labeled as D) and clutter from the Adriatic Sea (labeled as E). From Figure 5 onwards, radar plots for this work have been created by interpolating data on a regular Cartesian grid and in particular using the projected coordinate reference system (CRS) EPSG:3004 (European Petroleum Survey Group, https://epsg.org/home.html), also named "Monte Mario / Italy zone 2 -Datum: Roma 40 -Projection: Gauss-Boaga -Zone: East". The area of use for EPSG:3004 is Italy -East of 12°onshore and offshore. The proposed processing chain was designed to save output radar products in Geographic Tagged Image File Format (GeoTiff, https://www.ogc.org/standard/geotiff) so that any geographical information system (GIS) software can easily ingest them. GeoTiff standard requires metadata providing projection information for the exchange of georeferenced and geocoded imagery so we chose EPSG:3004 as projected CRS for output radar products related to Central Eastern Italy. A scheme of Cingoli radar filtering and processing chain is shown in Figure 6. The processing chain is divided into four main sections: residual ground clutter filtering, attenuation correction, separation of meteorological from non-meteorological targets and removal of Racon echoes. The processing chain receives as input an HDF5 file with the entire polar radar volume and it starts working on Doppler-filtered reflectivity (Z) data. Before describing in detail the individual steps of the chain, we report the final result obtained by applying the algorithms used for the present work to the previously mentioned case; Figure 7 then shows the VMI (on the left) and the SRI (on the right), two of the outputs of the processing chain, related to 2021-7-19 at 1:30 UTC. SRI is calculated starting from VMI data through the power law proposed by Marshall and Palmer, consistent with an exponential drop-size distribution [13]: where a and b are adjustable parameters while Z and R are the reflectivity factor in mm 6 m -3 and the rainfall rate in mm h -1 , respectively; according to Marshall and Palmer study, in this work we used a=200 and b=1.6. It is evident from Figure 7 how the processing chain is able to suppress the various sources of noise while preserving the meteorological signal at the same time.   Figure 7). Furthermore, maximum rainfall intensity data coming from the regional meteorological and hydrological monitoring network (SIRMIP) and measured every 15 minutes, confirm that on 2021-7-19 at 1:30 UTC the only measured values greater than 0 mm per minute are related to rain gauges located in the South of Marche region. SIRMIP is managed by Marche Region Civil Protection Service [14] making data freely available through the SIRMIP on-line (SOL) web application at http://app.protezionecivile.marche.it/sol [15].

Ground Clutter
Weather radars are often installed on the top of a hill or a mountain, with the purpose of surveillance of complex-orography territories; this is also the case of Cingoli radar, installed at about 750 meters on the Apennine chain. Since a radar scan from a high site can imply a large number of pixels affected by groundclutter, its rejection is then essential for the use of radar data, both for quantitative and qualitative purposes. The texture-based technique proposed by Gabella and Notarpietro [16] and used in the present work is based on the fact that non-stationary ground clutter and anomalously propagated echoes decorrelate rapidly in space and are spatially heterogeneous. The signature of such echoes may be recognized in reflectivity data, as their spatial variability is larger than the weather echoes. The filtering technique focuses on the horizontal spatial variability of the radar reflectivity field. Gabella algorithm can be divided into two logical parts: the first part consists of a spatial-proximity filter and the second is a test of compactness. The former is a consequence of the larger spatial continuity of precipitation fields than ground clutter, and the latter of the different spatial characteristics of residual clutter and anomalously propagated echoes with respect to precipitation field. The first filter is applied to each pixel in order to eliminate data that are weakly spatially correlated to the surrounding ones. For this purpose, a square window of side length wsize pixel around the considered polar bin was used. The reflectivity value of the pixel is assumed to be a meteorological echo when the differences between it and np surrounding pixels in the window are below a certain threshold Tr1; otherwise, the pixel is considered to be affected by ground clutter and a flag replaces its value. The second filter identifies a minimum echo area considering adjacent pixels of not null intensity: the pixels are considered to belong to the same group if they touch on any of the eight possible directions (including diagonal directions). The ratio Rp of the total number of pixels in the group and the number of pixels defining its boundary is evaluated for each group and a given group is classified as unwanted clutter when Rp is less than a minimum threshold Tr2. As suggested by Gabella, for this work we used the following value for the filters: wsize = 5, Tr1 = 6 dBZ, np = 8 and Tr2 = 1.3. The choice of Tr2 implies that the smallest compact group accepted as a meteorological feature must have more than 11 pixels.
In order to show the performance of the processing chain in its steps, the polar volume measured by the Cingoli radar on 2022-6-15 at 16 UTC will be considered. During that day, Italy was positioned on the eastern edge of a high-pressure area on which northwesterly currents flow at high altitude. On the ground, the baric gradient is very low and the circulation follows the regime of sea breezes that, combined with daytime heating, contribute to the formation of cumulus clouds on the Apennine ridge in the early afternoon.
Forecasted vertical profile data from European Centre for Medium-Range Weather Forecasts (ECMWF, https://www.ecmwf.int) and reported in Figure 9 for 2022-6-15 at 12 UTC shows a low predisposition to deep convection due the high level of free convection, so only some cumulus clouds evolved in severe storm. In particular, an intense storm system remained stationary upon Umbria region and northwest of Marche region, bringing heavy rains from 14:45 to 17:00 UTC. In Figure 10 Figure 11 shows the reflectivity Doppler-filtered by Eldes (input of the processing chain, on the left) at elevation 0°and the reflectivity filtered according to the Gabella algorithm for residual ground clutter removal (right). It is evident from right side of Figure 11 how the Gabella filter is able to remove the clutter echoes present in the left side of Figure 11 and especially in the northeastern quadrant.

Attenuation Correction
Radar data was gate-by-gate corrected for the attenuation using the iterative approach of Jacobi and Heistermann [17] that represents a modification of the algorithm proposed by Kraemer and Verworn [18]. The algorithm used aims at finding the optimal combination of the parameters c and d that relate the specific attenuation k (in dB/km) and radar reflectivity Z via a power law: d k = c Z (2) The optimization process started with a couple of initial values for linear coefficient c and exponential coefficient d, i.e. cmax and dmax. The coefficients are then successively reduced for contiguous beam sectors to minimal allowed values cmin and dmin until stability is achieved. The algorithm also fixed the number of iterations cn between cmax and cmin and dn between dmax and dmin. As a criterion to detect instability, the attenuation correction along a beam is considered as unstable if the value of attenuation-corrected reflectivity in any gate exceeded a maximum value Zcorr.max or if the path-integrated attenuation (PIA) exceeds a maximum value PIAmax.
For the present work we used the following values: cmax = 1.67 10 -4 , cmin = 2.33 10 -5 , cn = 4, dmax = 0.7, dmin = 0.65, dn = 6, Zcorr.max = 59 dBZ and PIAmax = 20 dB; furthermore, as reported in Table 2, gate resolution of Cingoli radar data is 125 m. The attenuation-corrected reflectivity according to the Jacobi and Heistermann algorithm is shown in Figure 12; as can be observed, reflectivity values related to convective cells are increased with respect to the right side of Figure 11, i.e., after the application of Gabella filter.

Separation and Removal of Non-meteorological Targets
In order to separate meteorological from non-meteorological features in radar data and therefore be able to discard the latter, we decided to implement the methodology proposed by Kilambi et al. [19] in the framework of the processing chain described here. The methodology is based on a commonly observed property of meteorological targets compared to other targets, i.e. they present simultaneously high values of correlation coefficient (ρHV) and values of differential reflectivity (ZDR) close to 1 in linear units. In other words, meteorological targets are basically uniform in shape and close to spherical. Starting from ρHV and ZDR data measured by dual-polarization systems such as Cingoli radar, depolarization ratio (DR) can be derived: where ZDR and DR are in linear units. DR has values between 0 and 1 but we converted into decibel (dB) units for the implementation of the algorithm. According to Kilambi et al., non-meteorological radar echoes are characterized by values of DR > -12 dB and reflectivity Z < 35 dBZ. This simple method has proven effective for identifying non-meteorological targets at C-as well as at X-band, so we used it for Cingoli radar data. Since we measured echoes from light rain having reflectivity values of about 15 dBZ, a less stringent condition for the reflectivity (Z < 15 dBZ) was adopted for this work. Misclassifications between weather and nonmeteorological echoes are typically caused by isolated pixels in the melting layer or at the edge of echo patterns [19] then, as suggested by Kilambi et al., we adopted a despeckling algorithm in order to reduce such misclassifications. A window of width 3 around the pixel has been used for despeckling. In Figure 13 the reflectivity related to radar data measured on 2022-6-15 at 16 UTC (elevation 0°) is reported after removal of non-meteorological targets by means of depolarization ratio and despeckling. The effect of such algorithms is not so evident from the Figure because they are devoted to isolated pixels, operating in cascade to two other filters, Doppler by Eldes and Gabella, all aimed at removing non-meteorological from meteorological echoes. It is nevertheless important to point out that, as expected, applying depolarization ratio and despeckling algorithms preserves all meteorological features present in radar data. Ground clutter filtering, attenuation correction and algorithm for separation of meteorological from nonmeteorological targets were applied to the reflectivity of each elevation of the radar volume.

Racon Filtering
Randomly, Cingoli radar data present a strong echo from the Adriatic Sea most likely due to a Racon. A Racon is a radar transponder commonly used for maritime navigation to mark hazards; when a Racon receives a radar pulse, it responds with a user-programmed signal at the same frequency placing a characteristic signature on the radar display (see label D North-East from radar site on Figure 5). In its normal operation, a Racon is activated by X-band radars present on ships as an aid to navigation.
Echoes from Racon are observable on Cingoli radar data at the first two elevations only and they always appear in the same geographical position at both elevations. In polar coordinates, Racon echoes have been identified in the area bounded by azimuth between 17 and 23 degrees and range bin between 657 and 734 (i.e., between 82.125 and 91.75 km away from the radar site). The echoes from Racon can appear by chance at both elevations or in only one of the two; however, it is eliminated or at least reduced by means of the filter created ad-hoc for the present work.
We have assumed that echoes outside the identified area are not affected by Racon so we focused our attention only on radar bins within the Racon area and belonging to the first three elevations.
Basically, for each pixel with a valid reflectivity value present in the Racon area and in the selected elevations, the corresponding pixels along its vertical are examined; if at least a pixel has a not valid reflectivity value then the pixel under examination is classified as affected by Racon and then discarded. Removal of Racon echoes is evident from VMI related to data measured on 2021-7-19 at 1:30 UTC and shown in Figure 7 compared to unfiltered and Doppler-filtered reflectivity in Figure 5.
In Figures 14 and 15, the SRI and the three-dimensional VMI related to radar data measured on 2022-6-15 at 16 UTC are shown as output of the processing chain. The vertical structure of the precipitating system up to 12.5 km is clearly visible from Figure 15, obtained by interpolating polar coordinates of all elevations into EPSG:3004 Cartesian grid with horizontal and vertical resolution of 125 m and 500 m, respectively.
Real-time VMI and SRI, two of the outputs of the processing chain described here, are available online at Marche Region Civil Protection Service webpage http://console.protezionecivile.marche.it/MonitoraggioMeteo/ observations/radar/html/radar.html.

Conclusion
In this work a data-processing chain developed ad hoc in order to deal with the sources of error affecting Cingoli X-band weather radar data has been presented in detail. The software was written using open source technology and is operationally used by the Civil Protection Service of Marche Region, Italy. The processing chain has been found to be able to greatly reduce or even suppress the unwanted radar echoes while preserving meteorological features at the same time. The performance of the data-processing chain has been evaluated in the light of case studies related to meteorological events that recently interested Marche region territory.
Cingoli radar data present clutter echoes received from the Adriatic Sea at the first two elevations of the scan. Measured Doppler velocity pattern due to sea clutter is similar to that due to the presence of moderate rainfall so Doppler filtering cannot work properly. Furthermore, operational and research weather radars are often covered by a radome that serves several purposes such as protecting the antenna and pedestal from weathering, providing a consistent environment and a direction-independent wind load. In addition to these technical reasons, radomes also mitigate the visual impact of operational weather radars by hiding the moving parts. Despite all this, radomes have also the deleterious effect of disrupting the transmitted and received radar signals; in particular, rain and ice on the radome degrade the signals beyond the nominal dry-radome attenuation, accurately measured by the radar vendor. For frequencies below C-band, these impacts can usually be ignored but at higher frequencies, and in particular at X-band, the impact of radome attenuation can be significant.
Since Cingoli radar is about 30 km away from the Adriatic Sea and is covered by a radome, the processing chain described here will be reviewed and updated in a future work in order to take into account the effects of sea clutter and wet radome attenuation on radar data.