Introduction to bioRad

Step 1: loading and visualizing radar scans

The low-level radar data with which bioRad interacts are so-called polar volume data. A polar volume is a collection of full-circle azimuthal scans (also referred to as sweeps) at various elevations of the radar antenna, which together provide a sampling of the atmosphere at all altitudes of interest. bioRad reads polar volumes with the read_pvolfile() function, which returns the polar volume as an object of class pvol. bioRad currently supports HDF5 files (Michelson 2014) that are compliant with the European OPERA Data Information Model (ODIM) (OPERA: Operational Program for Exchange of Weather Radar Information; see Huuskonen et al. 2014), and level-2 data generated by the US Next Generation Weather Radar (NEXRAD) network.

A polar volume (class pvol) contains a list of scans (class scan), each of which consists of a list of scan parameters (class param). A scan parameter is one of the radar’s basic observed quantities, such as reflectivity factor and radial velocity, and for dual-polarization radars additional quantities such as correlation coefficient, differential phase, and differential reflectivity.

Scan parameters can be projected on a georeferenced Cartesian grid in the form of a plan position indicator (PPI) objects (class ppi) using the function project_as_ppi(). These can either be plotted directly using the function plot.ppi() (Fig. 2B and 2C) or overlayed on a customizable basemap using the function map() (Fig. 2D and 2E), which makes use of the ggplot2 (Wickham 2016) and ggmap (Kahle and Wickham 2013) R libraries.

Step 2: processing volumes into vertical profiles

Polar volumes can be processed into vertical profiles of biological targets using the calculate_vp() function, which is a release of the algorithm vol2bird (Dokter et al. 2011), available independently on Github. The function takes in a polar volume file and outputs a vertical profile file and/or a vertical profile (vp) class object. The function has an argument autoconf, which when set to TRUE will select default settings automatically (depending on radar wavelength and availability of polarimetric data).

We describe the most important algorithm parameters and their preferred settings:

1. range_min, range_max: sets the minimum and maximum range (distance from the radar) of data to include. We recommend a minimum range of 5 km, to exclude the closest ranges that typically contain a lot of ground clutter. We recommend a maximum range of 35 km, which for most radars allows coverage up to 3 to 4 km a.s.l., which is the altitude band in which most migration occurs (Bruderer et al. 2018). At longer ranges, the radar beam gets very wide, hampering the radar’s ability to resolve altitudinal distributions.
2. layers, layer_height: sets the number of altitude layers and their thickness, respectively. Altitudes are defined relative to mean sea level, taking into account the antenna height as stored in the original polar volume file. We recommend a thickness of 200 m. Profiles with narrower altitude bin spacings can be extracted (Buler and Diehl 2009), but the finite size of the radar beam precludes resolving altitudinal features smaller than approximately 100–200 m. Profile quantities are estimated based on resolution samples centered within the altitudinal spacing of each layer.
3. dual_pol, rho_hv: the logical dual_pol enables polarimetric filtering of precipitation, which discards contiguous areas with correlation coefficient (ρ_HV) above a threshold rho_hv. We recommend rho_hv = 0.95, since precipitation typically has higher correlation coefficient values (Stepanian et al. 2016) (but note that lower ρ_HV is possible in mixed precipitation, like a combination of snow and rain, cf. Ryzhkov and Zrnic 1998). Single polarization mode is currently only available for C-band radars.
4. dealias, nyquist_min: the logical dealias enables radial velocity dealiasing following the method by Haase and Landelius (2004) when scans are present with a Nyquist velocity smaller than threshold nyquist_min (default 25 m s^–1). The Nyquist velocity is stored in the attributes$how$NI slot of scan class objects. Some radars dealias velocities at acquisition time, e.g. using the dual-PRF technique (Holleman 2005). For such radars we recommend no dealiasing for scans on which this is applied. For data acquired with a single PRF we recommend dealiasing when the Nyquist velocity of a scan is below 25 m s^–1, i.e. if there is a high probability that animal movements will be faster than the Nyquist velocity.
5. sd_vvp_threshold: animal speed and direction are estimated using the Volume Velocity Profiling (VVP) technique (Waldteufel and Corbin 1978, Holleman 2005). VVP also provides the standard deviation of the fit residuals (quantity sd_vvp in a profile). The sd_vvp_threshold parameter sets the threshold for discarding data based on this standard deviation measure. Animal density will be set to zero in altitude layers with a VVP standard deviation sd_vvp < sd_vvp_threshold. We recommend applying this thresholding as a way of removing residual rain contaminations and insects in bird studies using C-band radars, where sd_vvp_threshold = 2 m s^–1 was shown a suitable value (Dokter et al. 2011). We note that sd_vvp may become large in relatively rare cases where the velocity field is highly nonlinear (e.g. strong shear), causing this thresholding criterion to break down. For S-band radars VVP standard deviation thresholding has not been thoroughly evaluated, but radial velocity variability during bird migration may be lower than at C-band in certain cases. We currently recommend a conservative threshold of 1 m s^–1 to retain more biological scatter.
6. rcs: value for the radar cross section (RCS) of an individual. We recommend 11 cm² as a starting point, which was the seasonal average for C-band radars in western Europe during nocturnal passerine migration, according to a calibration experiment (Dokter et al. 2011). Note that radar cross sections depend on target size, body orientation, and radar wavelength (Vaughn 1985).

The sd_vvp_threshold and rcs parameters can be changed using the sd_vvp_threshold() and rcs() functions (in step 3 and up) without having to reprocess the vertical profile (step 2).

Step 3: visualizing and interpreting individual profiles

The various quantities in a vertical profile (e.g. dens: animal density, ff: ground speed, dd: ground speed direction, eta: reflectivity) can be visualized with plot.vp(), as shown in Fig. 2F and 2G for density. These profile plots and Fig. 2D and 2E are for the same moment in time. Note that both profiles show layering of birds: a density concentration at high altitude (here at approx. 1.5 km) (cf. Dokter et al. 2013). These layers show up as concentric rings in Fig. 2D and 2E. These rings appear because at an increasing distance from the radar, measurements are made at higher altitudes, because of the positive beam elevation and the curvature of the earth.

Also note that the peak densities of the two cases are similar, on the order of 100 individuals km^–3 (assuming RCS = 11 cm²) (Fig. 2H and 2I). The reflectivity factors (in dBZ scale, not to be confused with reflectivity η (Dokter et al. 2011, Chilson et al. 2012b) are however much higher for the US case than the European case. This is related to the difference in radar wavelength (Dokter et al. 2011), with NEXRAD radars in the US being S-band and European radars being mostly C-band.

Step 4: analyzing and visualizing vertical profiles as time series

After processing volume data into profiles, the profile data of consecutive volume scans of a radar can be organized into a time series of vertical profiles. The function bind_into_vpts() binds vertical profile objects (class vp) into time series object (class vpts), for which the default plot is shown in figure 2H and 2I. The dotted line indicates the time slice of Fig. 2B–G.

The plot.vpts() method overlays one of the reflectivity-based quantities (e.g. dens, eta or dbz) with a barb indicating the animals ground speed and direction. This follows meteorological conventions for graphically displaying wind speed and direction (with north being up). The number of barb flags indicate the speed (ff) while its tip points into the direction where animals are moving (dd).

Another useful profile quantity to inspect as time series is DBZH. This is the reflectivity factor for all scatterers, including meteorological targets like precipitation. Time periods with rain are often clearly visible as high DBZH values over the full altitude column. We recommend making plots of DBZH as a way of screening for precipitation contaminations and quality control, which is often a useful way to check remarkable altitude patterns in the biological data (e.g. the layering of birds at 1.5 km can also be seen in Fig. 2I) or short spikes with high values that might be due to rain contamination.

bioRad provides multiple functions to further aggregate and summarize time series data. We can integrate over the altitude dimension using integrate_profile(), which outputs a specially classed data.frame (class vpi) containing altitudinally integrated or averaged quantities. Figure 2J and 2K show plots of migration traffic rate, both MTR (variable transect angle, Eq. 1) and MTR₀ and MTR₉₀ (fixed transect angle, Eq. 2). We note that MTR is always positive, but MTR_α definitions can become negative depending on the migratory direction in relation to α. For example, the northward spring migration (US case, Fig. 2K) result in a positive MTR₀, while the southward autumn migration (European case, Fig. 2J) is negative. For the US case, migration is directed mostly northward, therefore MTR₀ is much larger than MTR₉₀, while in the European case, migration is mostly westward, therefore (in absolute value) MTR₀ is smaller than MTR₉₀.

Vertically-integrated time series can be further accumulated in time into measures summarizing migration traffic having passed the radar station during a time period, like MT (migration traffic) or RT (reflectivity traffic) (cf. output columns mt and rt of integrate_profile()). For example, for the European case we find MT = 55 × 10³, MT₀ = –28 × 10³ and MT₉₀ = –45 × 10³ for the time night-time period. This means that – assuming a radar cross section (RCS) per individual of 11 cm² – 55 thousand birds per 1 km transect flew over the radar station in this night (irrespective of direction). Decomposing the migration traffic into two perpendicularly oriented components, we find a net 28 thousand birds flew southward per km over a west-to-east transect (MT₀), and a net 45 thousand birds per km flew westward per km over a north-to-south transect (MT₉₀). For these specific definitions, MT ≤ √(MT₀² + MT₉₀²), with the left- and right-hand side being equal when migration directions dd all point into a sector of at most 180 degrees wide, as is usually the case for periods confined to a single spring or fall.

Both the vp and vpts class objects can be exported to standard R data frames (using as.data.frame.vp() and as.data.frame.vpts()) for further analysis outside of bioRad.