bioRad: biological analysis and visualization of weather radar data

from both the European radar network (OPERA) and the radar network of the United States (NEXRAD). bioRad aims to make weather radar studies in ecology easier and more reproducible, allowing for better inter-comparability of studies.


Introduction
Weather surveillance radars continuously survey the airspace of many countries around the globe to detect precipitation and severe weather. This meteorological infrastructure also has a great and still underappreciated potential for quantifying biological phenomena in the airspace (Chilson et al. 2012a, Shamoun-Baranes et al. 2014, Bauer et al. 2017. Weather radar can measure aerial movements of various biological taxa, including birds (Gauthreaux Jr and Belser 1998), bats (Stepanian and Wainwright 2018) and insects (Rennie 2014). Because of their year-round operation and organization in networks with continental-scale coverage, radar networks can provide standardized monitoring data at unprecedented temporal and spatial scales.
Following a proliferation of advances in information technologies, data infrastructure, and open data policies, access to low-level weather radar data has greatly improved over the last decade (Huuskonen et al. 2014, Ansari et al. 2018. These low-level data consist of scans (sweeps) in polar coordinates of each of the observed quantities by the radar, collected at multiple beam elevations (in the European OPERA network called single-site polar volumes, in the US NEXRAD network called level II data). Large advances have also been made in the development of methods to extract biologically relevant information from low-level radar data (Dokter et al. 2011, Stepanian andHorton 2015). This technological and methodological push, combined with an increasing need to understand and predict how animals are using the airspace, have led to a steep increase in the use of weather radar in ecology over the last decade.
Analysis of weather radar data for biological purposes has remained challenging nonetheless, requiring a variety of computer and programming skills as well as a basic understanding of how radars sample the atmosphere. Here we aim to improve the accessibility to tools and methods for biological analysis of weather radar data through the bioRad package for R (R Core Team), arguably the most widely used highlevel open source software language in biology and ecology. This paper describes a roadmap for analyzing weather radar data using bioRad. It also provides an overview of the various measures found in the literature for quantifying animal movement using weather radar and gives some good practices for reporting these measures.

Basic weather radar measures of animal movement
The movements and amount of animals in the airspace are often summarized in terms of vertical profiles. Vertical profiles can be generated by bioRad from low-level radar data and provide for each altitude above mean sea level (ASL) quantities like ground speed (ff), ground speed direction (dd), reflectivity (η), and animal density (dens). These profile quantities can be combined into multiple measures summarizing the number and passage of animals aloft. In the literature a large variety of measures can be found to report the amount of biological targets detected in the airspace by radar, like reflectivity factor z (Buler and Diehl 2009), reflectivity η (Dokter et al. 2011, Chilson et al. 2012b, vertically integrated reflectivity VIR (Gasteren et al. 2008, Shamoun-Baranes et al. 2011, McLaren et al. 2018, vertically integrated density VID (Buler and Diehl 2009, Dokter et al. 2011, Horton et al. 2014, migration traffic rate MTR (Nilsson et al. 2019), or migration traffic MT (Dokter et al. 2018). This section and Fig. 1 give an overview of these measures and their interrelation. Throughout this paper we provide recommendations for when to use which measure, and how to report them in a standardized way using bioRad.
While in radar meteorology reflectivity factors z (or Z in dBZ) are the conventional unit (for its useful property of being independent of radar wavelength in the case of small scatterers like precipitation, Doviak and Zrnić 1993), for larger animals like birds a more useful unit is reflectivity η (Dokter et al. 2011, Chilson et al. 2012b, which is more directly proportional to aerial animal density (see caption Fig. 1

for conversions).
A first choice is whether to use measures that are closely related to the reflectivity measurements of the radar (Fig. 1, left box), or measurements that are explicit in the numbers of individuals aloft (Fig. 1, right box). The advantage of reflectivity-explicit measures (Fig. 1, left box) is that they do not rely on assumptions of how to convert reflectivity to aerial animal densities, which may be information that is not available or has high uncertainty. The disadvantage is that these measures are less readily interpretable from a biological point of view. Individual-explicit measures (Fig. 1, right box) require knowledge of or explicit assumptions about the typical radar cross section (RCS) of individuals aloft (Vaughn 1985, Dokter et al. 2011, Mirkovic et al. 2016, Drake et al. 2017. The RCS of an object is the apparent area from which the object back-scatters radar waves emitted by the radar. It depends on the object's refractive index, shape and radar wavelength (Vaughn 1985). RCS also varies with aspect angle (body orientation relative to the radar beam), but since profile data is usually averaged over all azimuths, we can suffice with a single average RCS value for a given animal or animal type. When reporting numbers of individuals, it is important to always report accompanying RCS values. For C-band radars in western Europe a seasonal average RCS of 11 cm 2 has been determined in a calibration experiment (Dokter et al. 2011), which we recommend as a good starting point for nocturnal migration of passerines. This value may be refined using more detailed knowledge about which species are migrating, e.g. from information on phenology or other independent measurements (Horton et al. 2018). A second choice is the level of data aggregation, with studies often presenting multiple levels of data aggregation. The most basic profile data is specific for a certain altitude and time ( Fig. 1, top row). Data can be summarized further firstly by accumulating over (a range of ) altitudes (Fig. 1, middle row), and secondly by accumulating data in time (Fig. 1, bottom row).

854
A third choice is whether to use measures that are dependent on the ground speed of the animals (Fig. 1, right column within boxes) or measures that are speed independent ( Fig. 1, left column within boxes). Especially in the context of animal migration, the number of animals passing through an area depend both on the density of animals aloft and their speed. All else being equal, higher speeds represent higher migration intensity since more animals fly through a given area per unit of time. Intensity measures that are products of ground speed and density are therefore common in the literature, most notably the migration traffic rate (MTR) (Lowery 1951, Bruderer 1971, Schmaljohann et al. 2008, for which we introduce here a reflectivity-based equivalent for weather radar (RTR, reflectivity traffic rate) (Fig. 1). Traffic rate measures have the important additional advantage of suppressing stationary (non-migratory) signal components in weather radar data: reflectivity signal components with zero velocity will bias velocity estimates down by the same amount as their contribution to the total reflectivity, hence, measures that are based on the product of speed and reflectivity, like MTR and RTR, are effectively insensitive to these zero-velocity signal components (Dokter et al. 2018).
The migration traffic rate (MTR) for an altitude band is effectively the number of individuals crossing a transect per unit of transect length (usually 1 km) and per unit of time (usually 1 h). In most studies the transect is taken perpendicular to the ground speed direction of movement. Defined as such, MTR is always a positive quantity, defined as: with t time, h 1 the lower altitude and h 2 the upper altitude of interest, and dens(t,h) and ff(t,h) the animal density and speed at altitude h and time t, respectively. Because the transect is perpendicular to the direction of movement, it rotates along with shifting ground speed directions of the animals. The transect direction can also be fixed to a single angle, in which case with dd(t,h) the ground speed direction and α the transect direction (Supplementary material Appendix 1 Fig. A1). The angle α starts at 0 for a west-to-east transect (which has a northward perpendicular direction) and are defined clockwise from north. Note that this equation evaluates to the previous equation when α = dd, as required. In this definition, MTR α is a classical flow rate, giving the numbers of individuals moving into a direction of interest per unit time  855 and per unit transect length. Individuals moving into the direction α contribute positively to MTR α , while targets moving in the opposite direction contribute negatively. MTR α can thus be positive or negative, depending on the direction of movement (cf. Fig. 2J and 2K). For a transect α = 0 in the northern hemisphere, spring migration is typically positive and autumn migration negative. By integrating the migration traffic rates over a time period (from time t 1 to t 2 ), we obtain the migration traffic: the number of individuals that passed the one km transect during the time period: The definitions of RTR and RT are identical to those of MTR and MT above, except density (dens) should be replaced by reflectivity (η). Instead of the numbers of individuals, these measures give the cumulative cross-sectional area crossing the transect per unit time (RTR), or in a period of time (RT). We recommend using the traffic measures dependent on transect angle α when estimating the actual passage across a geographic transect line of interest. Examples are the estimation of influx or efflux from a geographic region (Dokter et al. 2018a), or when comparing weather radar data to other sensors surveying along a stationary geographic line or plane, like a fixed vertically rotating ship radar (Fijn et al. 2015). The measures independent of transect angle are most appropriate when quantifying traffic irrespective of the direction of movement, e.g. when comparing the amount of migration across large areas over which the general direction of movement varies (Nilsson et al. 2019).

General package structure and functionality
The functionality of bioRad is summarized in Fig. 2. Essentially, the package allows users to: 1) Load, inspect and visualize low-level radar data (polar volume data, also called level-II data in the US) of C-band or S-band weather radars, formatted in either the European OPERA (ODIM hdf5) or US NEXRAD data standard.
2) Extract biological information (speed, direction and density) at different altitudes.
3) Visualize, aggregate, and summarize this biological information over specific altitudes and time periods.
In bioRad, class objects are used for storing low-level data and data products, shown as blue/green boxes in Fig. 2. R has multiple class object systems, and bioRad uses the S3 object system (Chambers 2016). Most of these class objects have an associated plot method for making quick visualizations. The right-hand side of Fig. 2 shows examples of the output of these plot methods, for two migration events of similar intensity, one in Europe and one in the US. bioRad is able to extract vertical profiles of speed, direction, and density at different flight altitudes from low-level radar data, while offering standardized tools for post-processing and further analysis. Spatial variation in the horizontal plane is averaged out in profiles, and data is usually processed up to 25-35 km from the radar. Vertical profiles are generated in bioRad with the vol2bird algorithm (available at < https://github.com/ adokter/vol2bird >), originally developed for single and dualpolarization C-band radars (Dokter et al. 2011).
For this publication the underlying C-code for the algorithm has been refactored for compatibility with European and US radar formats, and for improved structure and readability of the code base. Additional support has been added for dual-polarization S-band radars, like the US WSR-88D/ NEXRAD radars, as well for dealiasing radial velocities. The package does not yet support automated removal of precipitation signals for single-polarization S-band radar. For these radars the generated profiles should be manually screened for precipitation contamination (cf. step 4 analysis workflow).

Analysis workflow
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), cf. Fig. 1. 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 (Fig. 2B, C) or overlayed on a customizable basemap using the function map (Fig. 2D, 1E), which makes use of the ggplot2 (Wickham 2016) and ggmap (Kahle and Wickham 2013) R libraries.

Step 2: processing volumes into vertical profiles
Volumes can be processed into vertical profiles using the calculate_vp function, which is a release of the algorithm vol2bird (Dokter et al. 2011), available independently on Github (< https://github.com/adokter/vol2bird >). 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 (Supplementary material Appendix 1).
3) dual_pol, rho_hv: the logical dual_pol enables polarimetric filtering of precipitation, which discards contiguous areas with correlation coefficient (r 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 andCorbin 1978, Holleman 2005). VVP also provides the standard deviation of the fit residuals (see Supplementary material Appendix 1, 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 2 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, as shown in Fig. 2F and 2G for density. These profile plots and Fig. 2D, E 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 2 ) (Fig. 2H, I). 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.

858
The plot 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 (Fig. 1). Figure 2J and 2K show plots of migration traffic rate, both MTR (variable transect angle, Eq. 1) and MTR 0 and MTR 90 (fixed transect angle, Eq. 2). We note as before 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 0 , while the southward autumn migration (European case, Fig. 2J) is negative. For the US case, migration is directed mostly northward, therefore MTR 0 is much larger than MTR 90 , while in the European case, migration is mostly westward, therefore (in absolute value) MTR 0 is smaller than MTR 90 .
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 in Eq. 3 (cf. output columns mt and rt of integrate_profile). For example, for the European case we find MT = 55 × 10 3 , MT 0 = -28 × 10 3 and MT 90 = -45 × 10 3 for the time night-time period. This means that -assuming a radar cross section (RCS) per individual of 11 cm 2 -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 0 ), and a net 45 thousand birds per km flew westward per km over a north-to-south transect (MT 90 ). For these specific definitions, MT ≤ √(MT 0 2 + MT 90 2 ), 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) for further analysis outside of bioRad.

Conclusions, recommendations and outlook
bioRad provides a set of functions to extract biological information from weather radar data, to present the information in graphical form, and to aggregate it in useful summary statistics. bioRad streamlines the reporting of analysis results according to existing conventions in the literature.
For larger-bodied animals like birds, we recommend the following measures when reporting data in aggregated form: 1) To quantify the instantaneous intensity of migration, or other large-scale directed movements for which a rate of passage is of interest: MTR (if RCS unknown: RTR).
2) To quantify the number of migrants passing in a certain time period: MT (if RCS unknown: RT).
3) To quantify the instantaneous number of animals aloft: VID (if RCS unknown: VIR). This measure is especially useful for cases that lack a large-scale directed movement, for example at the moment of a synchronized exodus of flight (Shamoun-Baranes et al. 2011, Buler andDawson 2014), or near a roost (Stepanian and Wainwright 2018).
Traffic measures (MTR, RTR, MT, RT) can be conditional on the choice of a transect line across which animals are counted (angle α, cf. Eq. 2, 3). We therefore recommend reporting whether a fixed transect was used or not, and if so, its direction. These vertically-integrated quantities are also conditional on an altitude band, therefore it should be clear whether they refer to the full altitude column, or only part of it.
The bioRad R software package aims to facilitate radar aeroecology research and make weather radar a more accessible tool in biological research to a broad range of researchers. We encourage extending weather radar as a tool beyond its current main application in quantifying songbird migration, e.g. towards larger flocking and soaring birds, bats and insects. Quantification of the movements of these species groups will require further calibration experiments (Dokter et al. 2011, Nilsson et al. 2019) and theoretical simulation work (Mirkovic et al. 2016) to identify their radar signatures and validate quantifications, in which dual-polarization information will likely be invaluable (Stepanian andHorton 2015, Stepanian et al. 2016). We also encourage the use of this package in biological radar studies in other countries with extensive weather radar networks, such as Australia, Canada, China, India, Japan and Russia, whose data might be readily explored once converted to a standardized (OPERA or NEXRAD) format. Expanding the use of weather radar networks for biological studies around the world will require continuing improvements in access to data and standardization of data formats, and in raising awareness of the value of collecting, distributing and archiving clear-air biological data with radar operators. There is a broad potential for using weather data in fundamental ecological research, in 859 applications for mitigating wildlife-human conflicts, and in conservation (Bauer et al. 2017). We expect its use will therefore only increase in the near future, and we hope these software tools will facilitate the further adoption of weather radar in the toolkit of biologists, conservationists and policy makers alike.
To cite bioRad or acknowledge its use, cite this Software note as follows, substituting the version of the application that you used for 'version 0': Dokter