NZ Metservice Rain Radar Network


Sebastien Delaux


Author: Sebastien Delaux

This notebook describes how to work with the raw radar files produced by the MetServices rain radar network. It shows how to read them and interact with the key quantities to plot the data. The last part of this notebook shows how multiple radar files (files corresponding to different locations). Can be blended to create a composite 3D or 2D gridded version of the data.

MetService is in the process of deciding under which license the data will be made available. For now only a few sample files have been added to this repository for illustration purposes. Please contact the author to arrange access to a larger sample of files.

Python libraries

We install pyart which is a python library specilised in dealing with raw radar data. We also load matplotlib which will be used to generate plots.

import sys !{sys.executable} -m pip install arm-pyart matplotlib
Collecting arm-pyart Collecting matplotlib Using cached Collecting cycler>=0.10 (from matplotlib) Using cached Collecting pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 (from matplotlib) Using cached Collecting kiwisolver>=1.0.1 (from matplotlib) Using cached Collecting numpy>=1.11 (from matplotlib) Using cached Collecting python-dateutil>=2.1 (from matplotlib) Using cached Collecting six (from cycler>=0.10->matplotlib) Using cached Collecting setuptools (from kiwisolver>=1.0.1->matplotlib) Using cached Installing collected packages: arm-pyart, six, cycler, pyparsing, setuptools, kiwisolver, numpy, python-dateutil, matplotlib Successfully installed arm-pyart-1.11.1 cycler-0.10.0 kiwisolver-1.1.0 matplotlib-3.2.1 numpy-1.18.2 pyparsing-2.4.6 python-dateutil-2.8.1 setuptools-46.1.3 six-1.14.0
import pyart from matplotlib import pyplot as plt import numpy as np
## 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
/usr/local/lib/python3.6/dist-packages/_pytest/mark/ DeprecationWarning: The usage of `cmp` is deprecated and will be removed on or after 2021-06-01. Please use `eq` and `order` instead. @attr.s(cmp=False, hash=False)
Reading the raw radar data

The New Zealand MetService's radar data are generated by a network of 8 doppler radars spread around New Zealand. Each radar does a scan of it's surroundings every 7 to 8 minutes. The radar operate independantly and produce one file for every scan.

Let's read one of the example files and look at the structure of the data. First we use pyart to load the data in the shape of a radar object.

file = './sample_data/BOP170101000556.RAWSURV' radar =

The documentation for radar object just loaded can be found here:

The radar files contain polar data from two sweeps (two revolution of the radar beam each for a different vertical angle). Each sweep contains data over a number of rays (angle) and gates (segments in the radial direction). Informations about the content of the file can be found using the info function of the radar object.
altitude: data: <ndarray of type: float64 and shape: (1,)> long_name: Altitude standard_name: Altitude units: meters positive: up altitude_agl: None antenna_transition: None azimuth: data: <ndarray of type: float32 and shape: (840,)> units: degrees standard_name: beam_azimuth_angle long_name: azimuth_angle_from_true_north axis: radial_azimuth_coordinate comment: Azimuth of antenna relative to true north elevation: data: <ndarray of type: float32 and shape: (840,)> units: degrees standard_name: beam_elevation_angle long_name: elevation_angle_from_horizontal_plane axis: radial_elevation_coordinate comment: Elevation of antenna relative to the horizontal plane fields: total_power: data: <ndarray of type: float32 and shape: (840, 1000)> units: dBZ standard_name: equivalent_reflectivity_factor long_name: Total power coordinates: elevation azimuth range _FillValue: -9999.0 reflectivity: data: <ndarray of type: float32 and shape: (840, 1000)> units: dBZ standard_name: equivalent_reflectivity_factor long_name: Reflectivity coordinates: elevation azimuth range _FillValue: -9999.0 fixed_angle: data: <ndarray of type: float32 and shape: (2,)> long_name: Target angle for sweep units: degrees standard_name: target_fixed_angle instrument_parameters: unambiguous_range: data: <ndarray of type: float32 and shape: (840,)> units: meters comments: Unambiguous range meta_group: instrument_parameters long_name: Unambiguous range prt_mode: data: <ndarray of type: |S5 and shape: (2,)> comments: Pulsing mode Options are: "fixed", "staggered", "dual". Assumed "fixed" if missing. meta_group: instrument_parameters long_name: Pulsing mode units: unitless prt: data: <ndarray of type: float32 and shape: (840,)> units: seconds comments: Pulse repetition time. For staggered prt, also see prt_ratio. meta_group: instrument_parameters long_name: Pulse repetition time prt_ratio: data: <ndarray of type: float32 and shape: (840,)> units: unitless meta_group: instrument_parameters long_name: Pulse repetition frequency ratio nyquist_velocity: data: <ndarray of type: float32 and shape: (840,)> units: meters_per_second comments: Unambiguous velocity meta_group: instrument_parameters long_name: Nyquist velocity radar_beam_width_h: data: <ndarray of type: float32 and shape: (1,)> units: degrees meta_group: radar_parameters long_name: Antenna beam width H polarization radar_beam_width_v: data: <ndarray of type: float32 and shape: (1,)> units: degrees meta_group: radar_parameters long_name: Antenna beam width V polarization pulse_width: data: <ndarray of type: float32 and shape: (840,)> units: seconds comments: Pulse width meta_group: instrument_parameters long_name: Pulse width latitude: data: <ndarray of type: float64 and shape: (1,)> long_name: Latitude standard_name: Latitude units: degrees_north longitude: data: <ndarray of type: float64 and shape: (1,)> long_name: Longitude standard_name: Longitude units: degrees_east nsweeps: 2 ngates: 1000 nrays: 840 radar_calibration: None range: data: <ndarray of type: float32 and shape: (1000,)> units: meters standard_name: projection_range_coordinate long_name: range_to_measurement_volume axis: radial_range_coordinate spacing_is_constant: true comment: Coordinate variable for range. Range to center of each bin. meters_to_center_of_first_gate: [0.] meters_between_gates: [300.] scan_rate: None scan_type: ppi sweep_end_ray_index: data: <ndarray of type: int32 and shape: (2,)> long_name: Index of last ray in sweep, 0-based units: count sweep_mode: data: <ndarray of type: |S20 and shape: (2,)> units: unitless standard_name: sweep_mode long_name: Sweep mode comment: Options are: "sector", "coplane", "rhi", "vertical_pointing", "idle", "azimuth_surveillance", "elevation_surveillance", "sunscan", "pointing", "manual_ppi", "manual_rhi" sweep_number: data: <ndarray of type: int32 and shape: (2,)> units: count standard_name: sweep_number long_name: Sweep number sweep_start_ray_index: data: <ndarray of type: int32 and shape: (2,)> long_name: Index of first ray in sweep, 0-based units: count target_scan_rate: None time: data: <ndarray of type: float64 and shape: (840,)> units: seconds since 2017-01-01T00:05:56Z standard_name: time long_name: time_in_seconds_since_volume_start calendar: gregorian comment: Coordinate variable for time. Time at the center of each ray, in fractional seconds since the global variable time_coverage_start metadata: Conventions: CF/Radial instrument_parameters version: 1.3 title: institution: references: source: history: comment: instrument_name: b'radbp' original_container: sigmet sigmet_task_name: b'SURVEILLANCE' sigmet_extended_header: false time_ordered: none rays_missing: 0
Plotting the raw data

The raw data from the two sweeps are stored in a single two dimensional array. They can be plotted in the form of two polar plots.

%matplotlib inline ## Create nsweeps polar plots fig, axs = plt.subplots(1, radar.nsweeps, subplot_kw=dict(polar=True)) ## For each plot find start and end indices for ax, istart, iend in zip(axs.flat, radar.sweep_start_ray_index['data'], radar.sweep_end_ray_index['data']): # Convert azimuth to geometric angle geometric_angle = -(radar.azimuth['data'][istart:iend+1]/180.*np.pi+0.5*np.pi) # Create grid from indices R, T = np.meshgrid(radar.range['data']/1000., geometric_angle) # Plot data on grid im = ax.pcolor(T,R,radar.fields['reflectivity']['data'][istart:iend+1,:], vmin=-30, vmax=40) fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) cbar = fig.colorbar(im, cax=cbar_ax) cbar.set_label("Reflectivity (dBZ)", rotation=270)

Those two plots allow to get a good idea of the horizontal spreading of the reflectivity field. However the data recorded by the radar is three-dimensional and slightly spread over time. Indeed, as shows below it takes about 1 minute for the two sweeps to be recorded.

%matplotlib inline plt.plot(radar.time['data']) plt.ylabel("Time since start of scan (s)") plt.xlabel("ID for ray")
Text(0.5, 0, 'ID for ray')
Gridded composite radar fields

Though the radars are not fully in sync, the frequency at which they perform their scans is about 7 minutes and in many cases it is OK to assume that all the nearest measurments in time of all New Zealand radars can be collated together to create a snapshot representation of the reflectivity field all over New Zealand at a given time. In that case we have to depart from a polar few of the reflectivity field and blend the data from the different radar over a rectilinear grid. The present section gives an example of how two radar files, one coming from the Auckland radar (AKL...) and one coming from the Bay of Plenty radar (BOP...) can be blended together onto a single grid. In the process some minor filtering is applied to the data. An azimithal equidistant projection with longitude/latitude origin (175.5, -37.0) output coordinate system is used as well as 41 levels in the vertical. More information on the regridding function can be found here:

## Path to the 2 files to collate radar_files = ['./sample_data/AKL170101000558.RAWSURV', './sample_data/BOP170101000556.RAWSURV'] ## Loading the data radars = [ for file in radar_files] ## Set filters for the radar object filters = [] for radar in radars: gatefilter = pyart.correct.despeckle.despeckle_field(radar, 'reflectivity', threshold=-100, size=20) gatefilter.exclude_transition() gatefilter.exclude_above('reflectivity', 80) filters.append(gatefilter) ## Do regridding grid =, gatefilters=filters, grid_shape= (41, 2250, 2000), grid_limits= ([0.0, 20000.0], [-460000.0, 460000.0], [-400000.0, 400000.0]), fields=['reflectivity'], max_refl= 80., copy_field_data= True, grid_origin= (-37.0, 175.5), roi_func= 'dist_beam', min_radius= 500.0, h_factor= 1.0, nb= 1.0, bsp= 1.0 )
/home/sebastien/.local/lib/python3.6/site-packages/pyart/map/ DeprecationWarning: Barnes weighting function is deprecated. Please use Barnes 2 to be consistent with Pauley and Wu 1990. " Pauley and Wu 1990.", DeprecationWarning)

Now we can plot five cross-sections of the gridded data at 5 fixes vertical levels (every 10 vertical levels of the grid). And get an idea of the 3D structure of the data.

%matplotlib inline fig, axs = plt.subplots(1, 5) for ax, data in zip(axs.flat, grid.fields['reflectivity']['data'][::10,:,:]): ax.pcolor(data, vmin=0, vmax=35)

We can also plot cross-sections along the north-south and east-west directiontp get an idea of the vertical structure of the data

%matplotlib inline fig, axs = plt.subplots(1, 2) ax = axs.flat[0] ax.pcolor(grid.x['data']/1000., grid.z['data'][:30]/1000., grid.fields['reflectivity']['data'][:30,1125,:]) ax.set(xlabel="Distance from radar", ylabel="Height agl (km)") ax = axs.flat[1] plot = ax.pcolor(grid.y['data']/1000., grid.z['data'][:30]/1000., grid.fields['reflectivity']['data'][:30,:,1000]) ax.set(xlabel="Distance from radar")
[Text(0.5, 0, 'Distance from radar')]

3D data can be quite heavy to work with and something one only consider the maximum reflectivity over the atmospheric column.

%matplotlib inline plt.pcolor(grid.x['data'], grid.y['data'], np.amax(grid.fields['reflectivity']['data'][:,:,:], axis=0), vmin=0, vmax=35) plt.xlabel('Easting') plt.ylabel('Northing') cbar = plt.colorbar() cbar.set_label("Radius of influence", rotation=270)

One can always keep some information on the vertical structure by storing the height at which the maximum in reflectivity was found

max_reflectivity_height = grid.z['data'][np.argmax(grid.fields['reflectivity']['data'][:,:,:], axis=0)]
%matplotlib inline plt.pcolor(grid.x['data'], grid.y['data'], max_reflectivity_height) plt.xlabel('Easting') plt.ylabel('Northing') cbar = plt.colorbar() cbar.set_label("Height of maximum reflectivity (m)", rotation=270)

An additional field that is produced during the regridding process is the radius of influence. This quantify the distance between the data point in the grid and the instrument that recorded it. This is an important quantity as the resolution and hence quality of the data drops with the distance from the instrument. On the plot below one can clearly see that two sources were used to created the grid, their location and how the radius of influence increases with the distance from the radars.

%matplotlib inline plt.pcolor(grid.x['data'], grid.y['data'], grid.fields['ROI']['data'][0,:,:]) plt.xlabel('Easting') plt.ylabel('Northing') cbar = plt.colorbar() cbar.set_label("Radius of influence", rotation=270)