Filtering and retrievals on raw Swiss C-band data

In this exercice we will load raw unfiltered Swiss C-band data during a thunderstorm event and process it to ultimately estimate the precipitation intensities. The following topics will be tackled.

• Ground clutter detection

• Attenuation correction

• KDP estimation

• Hydrometeor classification

• QPE

# Imports

import numpy as np
import cartopy.crs as ccrs
import cartopy
import matplotlib.pyplot as plt

import pyart
pyart.config.load_config('mch_config.py')

Note that you can create your own Py-ART configuration file, which defines default field names, default colormaps, limits, and much more. This is the one we use at MeteoSwiss. You can then either load it at startup in your python code or define the environment variable PYART_CONFIG to point to your file in your work environment.

Reading the data

Let’s start by loading our radar file which is of the standard CFRadial type. It corresponds to the third sweep of the operational radar scans, which is a PPI at 1° elevation. It contains raw radar data (before pre-processing) at a resolution of 83 m. We then add the temperature obtained from the COSMO NWP model to our radar object (note that this temperature was previously interpolated from the model grid to the radar polar grid). Note that the freezing level is quite high in this example (around 4200 m.).

# Open radar file
file_radar = './data/exercice1_swiss_thunderstorm/MHL2217907250U.003.nc'
radar = pyart.io.read_cfradial(file_radar)

# Add temperature
temp = pyart.io.read_cfradial('./data/exercice1_swiss_thunderstorm/20220628073500_savevol_COSMO_LOOKUP_TEMP.nc')
radar.add_field('temperature', temp.fields['temperature'])

Ground-clutter and noise removal

Py-ART uses gatefilters which are a kind of mask to filter out problematic measurements. Most processing routines can take a gatefilter as input and will ignore pixels that were filtered out.

Here we create a gate filter based on the radar moments and their texture to filter out noise and ground clutter. Since we are interested in a strong thunderstorm we also extend this filter to remove all measurements with a SNR ratio of less than 10 dB.

gtfilter = pyart.filters.moment_and_texture_based_gate_filter(radar)
gtfilter.exclude_below('signal_to_noise_ratio', 10)

Let’s compare visually the reflectivity before and after filtering. Note that the plot function of Py-ART take a gatefilter as input.

fig, ax = plt.subplots(1,2, figsize=(10,6), sharex= True, sharey=True)
display = pyart.graph.RadarDisplay(radar)
display.plot_ppi('reflectivity', 0, vmin=0, vmax=60., ax = ax[0], colorbar_label = 'Raw')
display.plot_ppi('reflectivity', 0, vmin=0, vmax=60., gatefilter = gtfilter, 
                 ax = ax[1], colorbar_label = 'Filtered')
ax[0].set_xlim([-50,50])
ax[0].set_ylim([-50,50])
ax[0].set_aspect('equal', 'box')
ax[1].set_aspect('equal', 'box')

Here it is clear that most ground clutter (mostly north west and east of the radar), as well as noise have been filtered out.

Attenuation correction

We can expect strong attenuation behind a thunderstorm like this. So it is a good idea to try to correct for it. Knowledge of the specific attenuation can also be very insightful.

out = pyart.correct.calculate_attenuation_zphi(radar, fzl = 4200,
                           gatefilter=gtfilter,
                           phidp_field = 'uncorrected_differential_phase',
                           temp_field = 'temperature',
                           temp_ref = 'temperature')
spec_at, pia, cor_z, spec_diff_at, pida, cor_zdr = out
radar.add_field('corrected_reflectivity', cor_z)
radar.add_field('corrected_differential_reflectivity', cor_zdr)
radar.add_field('specific_attenuation', spec_at)

Here we use the Z-PHI method, which uses the relation between differential phase shift and specific attenuation. However it works only in the liquid phase. So you need to provide it either with a fixed freezing level height, a field of freezing level heights or a field of temperature. Here we provide the later.

This method provides us with 5 output variables

• specific attenuation dB/km

• path integrated attenuation dB

• corrected reflectivity dBZ

• differential specific attenuation dB

• path integrated differential attenuation dB

• corrected differential reflectivity (ZDR) dB

We will now plot the specific attenuation as well as the raw and corrected reflectivities.

fig, ax = plt.subplots(1,3, figsize=(16,6), sharex= True, sharey=True)
display = pyart.graph.RadarDisplay(radar)
display.plot_ppi('specific_attenuation', 0, vmin=0, vmax=1.5, gatefilter = gtfilter,
                     ax = ax[0])
display.plot_ppi('reflectivity', 0, vmin=0, vmax=60., ax = ax[1],  gatefilter = gtfilter,
                 colorbar_label = 'ZH with attenuation [dBZ]')
display.plot_ppi('corrected_reflectivity', 0, vmin=0, vmax=60., gatefilter = gtfilter,
                     ax = ax[2], colorbar_label = 'ZH attenuation corrected [dBZ]')
ax[0].set_xlim([-50,50])
ax[0].set_ylim([-50,50])
ax[0].set_aspect('equal', 'box')
ax[1].set_aspect('equal', 'box')
ax[2].set_aspect('equal', 'box')

We can clearly observe a strong specific attenuation within the thunderstorm as well as a significant difference in reflectivity before/after correction behind the thunderstorm to the west.

KDP estimation

Another very interesting radar variable is the specific differential phase shift KDP. Large KDP indicates the presence of large oblate drops and is linked to very strong precipitation. KDP is also needed for the hydrometeor classification algorithm. However KDP is not measured directly and needs to be estimated numerically from the raw differential phase shift (PHIDP). Py-ART provides three different retrieval methods. We will use the method by Maesaka et al. (2012) which is fast and robust but assumes KDP to be positive and is therefore limited to rainfall below the melting layer and/or warm clouds.

kdp, _, _ = pyart.retrieve.kdp_maesaka(radar, gatefilter = gtfilter,
                                       psidp_field = 'uncorrected_differential_phase')
radar.add_field('specific_differential_phase', kdp)

fig, ax = plt.subplots(1,1, figsize=(6,6))
display = pyart.graph.RadarDisplay(radar)
display.plot_ppi('specific_differential_phase', 0, vmin = 0, vmax = 10,
                 ax = ax,  gatefilter = gtfilter)

ax.set_xlim([-50,50])
ax.set_ylim([-50,50])
ax.set_aspect('equal', 'box')

/store/msrad/utils/anaconda3-wolfensb/envs/rainforest_tests/lib/python3.10/site-packages/numpy/core/fromnumeric.py:784: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray.
  a.partition(kth, axis=axis, kind=kind, order=order)

A look at the KDP field shows clusters of very large KDP (> 5 °/km) at the center of the thunderstorm.

Hydrometeor classification

The hydrometeor classification algorithm in Py-ART by Besic et al. (2016) uses ZH, ZDR, RHOHV, KDP and the temperature to classify hydrometeors into one of 8 classes:

• Ice hail, high density Graupel

• Melting hail

• Wet snow

• Vertically oriented ice

• Rain

• Rimed particles

• Light rain

• Crystals

• Aggregates

This algorithm requires centroids of polarimetric variables for the different hydrometeor classes. Below we provide it with centroids specifically suited for the radar of Monte Lema. If left empty, the algorithm will use default centroids at the right frequency band (X, C or S).

centroids = np.array([[13.8231,0.2514,0.0644,0.9861,1380.6],
[3.0239,0.1971,0.,0.9661,1464.1],
[4.9447,0.1142,0.,0.9787,-974.7],
[34.2450,0.5540,0.1459,0.9937,945.3],
[40.9432,1.0110,0.5141,0.9928,-993.5],
[3.5202,-0.3498,0.,0.9746,843.2],
[32.5287,0.9751,0.2640,0.9804,-55.5],
[52.6547,2.7054,2.5101,0.9765,-1114.6],
[46.4998,0.1978,0.6431,0.9845,1010.1]])


hydro = pyart.retrieve.hydroclass_semisupervised(radar, mass_centers = centroids,
                                 refl_field =  'corrected_reflectivity',
                                 zdr_field = 'corrected_differential_reflectivity',
                                 kdp_field = 'specific_differential_phase',
                                 rhv_field = 'uncorrected_cross_correlation_ratio',
                                 temp_field = 'temperature')

radar.add_field('radar_echo_classification', hydro['hydro'])

fig, ax = plt.subplots(1,1, figsize=(6,6))
display = pyart.graph.RadarDisplay(radar)
import matplotlib as mpl

labels = ['NC','AG', 'CR', 'LR', 'RP', 'RN', 'VI', 'WS', 'MH', 'IH/HDG']
ticks = np.arange(len(labels))
boundaries = np.arange(-0.5, len(labels) )
norm = mpl.colors.BoundaryNorm(boundaries, 256)

cax = display.plot_ppi('radar_echo_classification', 0, ax = ax,  gatefilter = gtfilter,
                 norm = norm, ticks = ticks, ticklabs = labels)

ax.set_xlim([-50,50])
ax.set_ylim([-50,50])
ax.set_aspect('equal', 'box')

Note that the plotting commands are slightly more complicated due to the categorical colormap.

A look at the hydrometeor classification reveals the presence of wet hail in the center of the thunderstorm surrounded by rain and by light rain. A few isolated pixels (unfiltered ground clutter) are also classified as hail. There was indeed intense hail at the ground on that day.

QPE

Py-ART provides several QPE algorithms but the most refined relies on the hydrometeor classification and uses different relations between radar variables and precipitation intensities within the different hydrometeor classes.

qpe = pyart.retrieve.est_rain_rate_hydro(radar, refl_field = 'corrected_reflectivity',
                                         hydro_field = 'radar_echo_classification',
                                         a_field = 'specific_attenuation',
                                         thresh=40)

radar.add_field('radar_estimated_rain_rate', qpe)

We will now plot the precipitation intensity on a Cartopy map and add some spatial features (land borders) using RadarMapDisplay

lon_bnds = [8.2, 9.5]
lat_bnds = [45.5, 46.5]

display = pyart.graph.RadarMapDisplay(radar)

fig = plt.figure(figsize=(8,8))
display.plot_ppi_map('radar_estimated_rain_rate', 0, vmin=0, vmax=120.,
          colorbar_label='', title='Precipitation intensity [mm/h]', gatefilter = gtfilter,
          min_lon = lon_bnds[0], max_lon = lon_bnds[1],mask_outside = True,
          min_lat = lat_bnds[0], max_lat = lat_bnds[1],
          lon_lines=np.arange(lon_bnds[0], lon_bnds[1], .2), resolution='10m',
          lat_lines=np.arange(lat_bnds[0], lat_bnds[1], .2),
          lat_0=radar.latitude['data'][0],
          lon_0=radar.longitude['data'][0], embellish=True)

states_provinces = cartopy.feature.NaturalEarthFeature(
                category='cultural',
                name='admin_0_countries',
                scale='10m',
                facecolor='none')
lakes = cartopy.feature.NaturalEarthFeature(
                category='physical',
                name='lakes',
                scale='10m',
                facecolor='blue')
rivers = cartopy.feature.NaturalEarthFeature(
                category='physical',
                name='rivers',
                scale='10m',
                facecolor='blue')
display.ax.add_feature(states_provinces, edgecolor='gray')
display.ax.add_feature(lakes, edgecolor='blue', alpha = 0.25)
display.ax.add_feature(cartopy.feature.RIVERS)

/users/wolfensb/pyrad/src/pyart/pyart/graph/radarmapdisplay.py:287: UserWarning: No projection was defined for the axes. Overridding defined axes and using default axes with projection Lambert Conformal.
  warnings.warn(

<cartopy.mpl.feature_artist.FeatureArtist at 0x2b9150c92f80>

Note that we didn’t estimate precipitation intensity at the ground but only aloft. Within the thunderstorm precipitation intensity is extremely high. This is likely too high because QPE in wet hail is very uncertain. However, even the operational QPE algorithm at MeteoSwiss estimated precipitation intensities at the ground close to 120 mm/h.