Processing of Doppler wind data from a Swiss volumetric scan

In this exercice we will load low-resolution filtered Swiss C-band data and process it to estimate a profile of horizontal wind aloft the radar. The following topics will be tackled.

• Dealiasing of radial velocity

• CAPPI plots and profiles

• PseudoRHI profiles

• Computation of a VAD (velocity azimuth display)

# Imports

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

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

Reading and dealiasing the data

The Swiss operational C-band radar network performs 20 PPIs at 20 different elevations (from 0.2° to 40°) every 5 minutes. All PPIs are stored in separate files.

We will thus read all 20 elevations one after the other for a given timestep and use the Py-ART function join_radar to merge them all into a single radar object which containes 20 sweeps.

To avoid using too much memory we will read the pre-processed MeteoSwiss radar data which has a resolution of 500 m in range.

At the end, we will also dealias the radial velocity field. Indeed, at lower elevations, low PRFs are used which result in low Nyquist velocities between 8-12 m/s. This means that a lot of folding will occur, especially in strong winds. In this example we use the simplest dealiasing method of Py-ART which performs by finding regions of similar velocities and unfolding and merging pairs of regions until all regions are unfolded.

# Read all 20 elevations for one timestep
files_radar = sorted(glob.glob('./data/exercice2_swiss_doppler/MLL221790725*'))
for i,f in enumerate(files_radar):
    radar = pyart.io.read_cfradial(f)
    
    if i == 0:
        radar_merged = radar
    else:
        radar_merged = pyart.util.join_radar(radar_merged, 
                                       radar)
        
corr_vel = pyart.correct.dealias_region_based(radar_merged)
radar_merged.add_field('corrected_velocity', corr_vel)

We will now plot the raw and dealiased velocities at two different elevations to see the effect of the correction.

fig, ax = plt.subplots(2,2, figsize=(10,10), sharex= True, sharey=True)
ax = ax.ravel()
display = pyart.graph.RadarDisplay(radar_merged)
display.plot_ppi('velocity', 2, vmin=-30, vmax=30., ax = ax[0], title='El=1 deg',
                 colorbar_label = 'Mean Doppler velocity (m/s)')
display.plot_ppi('corrected_velocity', 2, vmin=-30, vmax=30., title='El=1 deg',
                      ax = ax[1], colorbar_label = 'corr. Mean Doppler velocity (m/s)')
display.plot_ppi('velocity', 6, vmin=-30, vmax=30., ax = ax[2], title='El=4.5 deg',
                 colorbar_label = 'Mean Doppler velocity (m/s)')
display.plot_ppi('corrected_velocity', 6, vmin=-30, vmax=30., title='El=4.5 deg',
                      ax = ax[3], colorbar_label = 'corr. Mean Doppler velocity (m/s)')
ax[0].set_xlim([-50,50])
ax[0].set_ylim([-50,50])
for a in ax:
    a.set_aspect('equal', 'box')

Indeed the raw velocity shows alternating bands of negative and positive velocities which indicates aliasing. The dealiased velocity looks much less discontinuous. Note however that a major difficulty for these algorithms is the presence of isolated pixels, which tend to get arbitrary values as can be seen in the south of the radar.

CAPPI plots

We will now create a CAPPI (constant altitude PPI) of the reflectivity during this event. The idea is to interpolate the volumetric scan on a 3D Cartesian grid using the function grid_from_radars. Here we will create slices every 500 m from 500 m to 8000 m above the radar.

zmin = 500
zmax = 8000
ymin= xmin = -100000
ymax = xmax = 100000
lat = float(radar.latitude['data'])
lon = float(radar.longitude['data'])
alt = float(radar.altitude['data'])
# number of grid points in cappi
cappi_res_h = 500
cappi_res_v = 500
ny = int((ymax-ymin)/cappi_res_h)+1
nx = int((xmax-xmin)/cappi_res_h)+1
nz = int((zmax-zmin)/cappi_res_v)+1

cappi_zh = pyart.map.grid_from_radars(radar_merged, grid_shape=(nz, ny, nx),
        grid_limits=((zmin, zmax), (ymin, ymax),
                     (xmin, xmax)),
        fields=['reflectivity'])

Now we plot the reflectivity at 4 different altitudes (0.5, 3, 5.5 and 8 km), as well as a profile along as a W-E profile at the radar location throught the thunderstorm.

display = pyart.graph.GridMapDisplay(cappi_zh)
projection = ccrs.PlateCarree()
fig = plt.figure(figsize=(18,14))
ax = plt.subplot(221, projection = projection)
display.plot_grid('reflectivity',0, ax = ax, projection = projection)
ax = plt.subplot(222, projection = projection)
display.plot_grid('reflectivity',5, ax = ax, projection = projection)
ax = plt.subplot(223, projection = projection)
display.plot_grid('reflectivity',10, ax = ax, projection = projection)
ax = plt.subplot(224, projection = projection)
display.plot_grid('reflectivity',15, ax = ax, projection = projection)

ax = fig.add_axes([0.25, -0.20, .5, .25])
display.plot_latitude_slice('reflectivity', lon=lon, lat=lat, ax = ax)

/store/msrad/utils/anaconda3-wolfensb/envs/rainforest_tests/lib/python3.10/site-packages/cartopy/mpl/geoaxes.py:1785: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
  result = super().pcolormesh(*args, **kwargs)
/store/msrad/utils/anaconda3-wolfensb/envs/rainforest_tests/lib/python3.10/site-packages/cartopy/mpl/geoaxes.py:1785: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
  result = super().pcolormesh(*args, **kwargs)
/store/msrad/utils/anaconda3-wolfensb/envs/rainforest_tests/lib/python3.10/site-packages/cartopy/mpl/geoaxes.py:1785: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
  result = super().pcolormesh(*args, **kwargs)
/store/msrad/utils/anaconda3-wolfensb/envs/rainforest_tests/lib/python3.10/site-packages/cartopy/mpl/geoaxes.py:1785: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
  result = super().pcolormesh(*args, **kwargs)

We will now create a pseudo RHI (altitudinal cross-section through a set of PPIs) of the radial velocity and the reflectivity through the thunderstorm at azimuth 270° (to the west).

pseudorhi = pyart.util.cross_section_ppi(radar_merged, [270])

display = pyart.graph.RadarDisplay(pseudorhi)
fig, ax = plt.subplots(2,1, sharex=True,sharey=True, figsize= (10,10))
display.plot_rhi('corrected_velocity', ax = ax[0], vmin = -30, vmax = 30)
display.plot_rhi('reflectivity', ax = ax[1])
ax[0].set_ylim([0,20])
ax[0].set_xlim([0,100])

(0.0, 100.0)

The convention at MeteoSwiss (which is different from the one used by Py-ART) is that positive velocities are moving away from the radar. In this example we see clearly the downdraft in the center of storm and the updraft in its surroundings.

Velocity azimuth display (VAD)

We will now make a VAD retrieval to estimate the horizontal wind profile above the radar. This technique requires to have measurements in as many azimuths as possible and works better for stratiform rain when the radar coverage is wider. We will load data from a cold front event on the 13th July 2021 near the Albis radar (south of Zürich), that showed widespread precipitation around the radar.

Unfortunately the VAD estimation technique in Py-ART can process only one sweep at a time. So we will average the wind vectors obtained over all sweeps to obtain a more reliable estimate. Note that we skip the first four sweeps which are more prone to ground echoes and have a very low Nyquist velocity.

files_radar = sorted(glob.glob('./data/question_pyart_meteoswiss/MLA211941205*'))
for i,f in enumerate(files_radar):
    radar = pyart.io.read_cfradial(f)
    
    if i == 0:
        radar_merged = radar
    else:
        radar_merged = pyart.util.join_radar(radar_merged, 
                                       radar)
        
corr_vel = pyart.correct.dealias_region_based(radar_merged)
corr_vel['data'] *= -1 
radar_merged.add_field('corrected_velocity_neg', corr_vel)

zlevels = np.arange(100,5000,100)

vad = pyart.retrieve.vad_browning(radar_merged, 'corrected_velocity_neg', z_want = zlevels)

max height 2682.0  meters
min height -54.0  meters
max height 5272.0  meters
min height 2.0  meters
max height 7815.0  meters
min height 4.0  meters

/users/wolfensb/pyrad/src/pyart/pyart/retrieve/vad.py:464: RuntimeWarning: Mean of empty slice
  mean_velocity_per_gate = np.nanmean(velocities, axis=0).reshape(1, -1)
/users/wolfensb/pyrad/src/pyart/pyart/retrieve/vad.py:549: UserWarning: Warning: converting a masked element to nan.
  y_new[i] = np.average(y_in_window, weights=weights)

max height 10403.0  meters
min height 6.0  meters
max height 14259.0  meters
min height 10.0  meters
max height 18533.0  meters
min height 15.0  meters
max height 22800.0  meters
min height 19.0  meters
max height 27060.0  meters
min height 24.0  meters
max height 31310.0  meters
min height 28.0  meters
max height 35549.0  meters
min height 32.0  meters
max height 39777.0  meters
min height 37.0  meters
max height 43991.0  meters
min height 41.0  meters
max height 50263.0  meters
min height 47.0  meters
max height 58621.0  meters
min height 56.0  meters
max height 71002.0  meters
min height 69.0  meters
max height 87118.0  meters
min height 85.0  meters
max height 106740.0  meters
min height 105.0  meters
max height 125476.0  meters
min height 124.0  meters
max height 143293.0  meters
min height 143.0  meters
max height 160000.0  meters
min height 160.0  meters

orientation = np.rad2deg(np.arctan2(vad.u_wind, vad.v_wind))%360
speed = np.sqrt(vad.u_wind**2 + vad.v_wind**2)

Note that because the convention at MeteoSwiss is different than the one in Py-ART we have to flip the sign of the radial velocity field.

Finally we do a plot of the vertical profiles or horizontal wind speed and direction

fig,ax = plt.subplots(1,2, sharey=True)
ax[0].plot(speed*2, zlevels + radar_merged.altitude['data'])
ax[1].plot(orientation, zlevels + radar_merged.altitude['data'])
ax[0].set_xlabel('Wind speed [m/s]')
ax[1].set_xlabel('Wind direction [deg]')
ax[0].set_ylabel('Altitude [m]')
fig.suptitle('Wind profile obtained from VAD')

Text(0.5, 0.98, 'Wind profile obtained from VAD')

Now let’s compare this wind profile with the one recorded by the nearest radiosounding operated in Payerne (around 130 km west from the radar):

Though there are some discrepancies the match is not bad, given the distance and the very different ways of measuring wind!