Exercice Sample Solution
Load required libraries
import pyart
pyart.config.load_config('mch_config.py')
import numpy as np
import cartopy.crs as ccrs
import cartopy
import matplotlib.pyplot as plt
import glob
## You are using the Python ARM Radar Toolkit (Py-ART), an open source
## library for working with weather radar data. Py-ART is partly
## supported by the U.S. Department of Energy as part of the Atmospheric
## Radiation Measurement (ARM) Climate Research Facility, an Office of
## Science user facility.
##
## If you use this software to prepare a publication, please cite:
##
## JJ Helmus and SM Collis, JORS 2016, doi: 10.5334/jors.119
/srv/conda/envs/notebook/lib/python3.9/site-packages/requests/__init__.py:102: RequestsDependencyWarning: urllib3 (1.26.8) or chardet (5.2.0)/charset_normalizer (2.0.10) doesn't match a supported version!
warnings.warn("urllib3 ({}) or chardet ({})/charset_normalizer ({}) doesn't match a supported "
Load all radar files in /data/question_pyart_meteoswiss and merge them into one single radar object
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)
Perform attenuation correction of ZH, using a constant freezing level height of 2700 m.
# Compute attenuation
out = pyart.correct.calculate_attenuation_zphi(radar_merged, fzl = 4200,
phidp_field = 'uncorrected_differential_phase',
temp_ref = 'fixed_fzl')
spec_at, pia, cor_z, spec_diff_at, pida, cor_zdr = out
radar_merged.add_field('corrected_reflectivity', cor_z)
Estimate the QPE with a a polynomial Z-R relation.
qpe = pyart.retrieve.est_rain_rate_zpoly(radar_merged, refl_field = 'corrected_reflectivity')
radar_merged.add_field('radar_estimated_rain_rate', qpe)
Compute a CAPPI of the resulting radar estimate rain rate from 500 to 8000 m above the radar using a vertical resolution of 100 m and a horizontal resolution of 500 m at a x and y distance of up to 100 km to 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 = 100
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_qpe = pyart.map.grid_from_radars(radar_merged, grid_shape=(nz, ny, nx),
grid_limits=((zmin, zmax), (ymin, ymax),
(xmin, xmax)),
fields=['radar_estimated_rain_rate'])
Using numpy, perform a weighted average of all CAPPI levels using the weights. Finally display the results.
weighting = np.exp(-0.5* cappi_qpe.z['data'] / 1000)
qpe_ground = weighting[:,None,None]*cappi_qpe.fields['radar_estimated_rain_rate']['data']
qpe_ground = np.nansum(qpe_ground, axis = 0) / np.sum(weighting)
plt.pcolormesh(cappi_qpe.point_longitude['data'][0],
cappi_qpe.point_latitude['data'][0],
qpe_ground, vmax = 15)
plt.colorbar(label = 'Estimated rain rate at ground [mm/h]')
/tmp/ipykernel_1107/2553980288.py:6: UserWarning: The input coordinates to pcolormesh are interpreted as cell centers, but are not monotonically increasing or decreasing. This may lead to incorrectly calculated cell edges, in which case, please supply explicit cell edges to pcolormesh.
plt.pcolormesh(cappi_qpe.point_longitude['data'][0],
<matplotlib.colorbar.Colorbar at 0x7fec8ff9b520>
Now let’s compare this QPE with the operational QPE of MeteoSwiss at the same timestep:
The agreement near the radar is not too bad, even with such a simple aggregation method. In the south we see some large discrepancies. This is due to the fact that the operational QPE includes many additional steps. In this case, the difference is likely due to the correction for partial beam blocking that is performed by the operational QPE in this mountaineous region south of the radar.