High-Resolution Spatial Distribution of Bird Movements Estimated from a Weather Radar Network

Abstract

Weather radars provide detailed information on aerial movements of organisms. However, interpreting fine-scale radar imagery remains challenging because of changes in aerial sampling altitude with distance from the radar. Fine-scale radar imagery has primarily been used to assess mass exodus at sunset to study stopover habitat associations. Here, we present a method that enables a more intuitive integration of information across elevation scans projected in a two-dimensional spatial image of fine-scale radar reflectivity. We applied this method on nights of intense bird migration to demonstrate how the spatial distribution of migrants can be explored at finer spatial scales and across multiple radars during the higher flying en-route phase of migration. The resulting reflectivity maps enable explorative analysis of factors influencing their regional and fine-scale distribution. We illustrate the method’s application by generating time-series of composites of up to 20 radars, achieving a nearly complete spatial coverage of a large part of Northwest Europe. These visualizations are highly useful in interpreting regional-scale migration patterns and provide detailed information on bird movements in the landscape and aerial environment.

1. Introduction

Quantitative information about the spatio-temporal distribution of migratory birds is important for conservation and reducing human-wildlife conflicts [1,2]. Of the variety of methods used to study migration, only a few are suitable to provide insight into spatio-temporal changes in abundance of migrants, especially across large regions and across a broad range of avian taxa; one highly suitable method is the use of networks of weather radars, e.g., [3,4,5]. Existing radar infrastructure maintained by meteorological institutes provides information on bird migration over extended temporal and spatial scales. Weather radars are continuously operated and generally scan the airspace in cycles of five to 15 min. Recent studies using weather radars have provided insight into questions that require a broad-scale perspective on migration and movement, such as the influence of light pollution on broad-front migration [6,7], the seasonal changes in population sizes and routes [3], spatio-temporal variation in migration patterns [5], and the decline of migratory populations [8].

Networks of weather radars cover large parts of the globe, (see Figure 1 in [9]). In Europe, there are already 164 radars that regularly exchange data in the Operational Programme for the Exchange of Weather Radar (OPERA) network [10], while there are 160 in the North American Next-Generation Radar (NEXRAD) network [9], with more radars distributed in many other countries [11]. A quick calculation assuming a detection range of 100 km for biological targets while accounting for overlap shows that the current European network continuously surveys a large part of the European landmass (2/3). In North America (United States and Canada), coverage is one fourth due to a sparser radar network; globally, this is one eighth, (using the reported radars; [12]). Radar networks thus provide ongoing measurements of bird migration across continents. While separating bird echos from others, many studies using multiple weather radars have aggregated the spatially explicit radar images into one single vertical profile of migration speed, direction, and density for each altitude of interest per radar, e.g., [3,5,7,8]. Although profile-based analyses provide great insights on a continental scale, they limit the interpretation on a regional scale and remove any information about spatial heterogeneity within the field of view of a single radar. This interpretation on a regional scale is of great importance to gain insight into responses to habitat and land use, as well as to inform the fine-scale planning of wind energy or aviation flight planning [2].

An important factor complicating the interpretation of radar imagery is the change of observed airspace with distance from the radar. Each elevation scan samples only part of the airspace. As the distance from the radar increases, a given scan measures the airspace at higher altitudes above the Earth’s surface, because of the curvature of the Earth and the elevation angle of the beam. Furthermore, the beam width increases with distance; therefore, measurements at larger distances sample a wider height band. Even the lowest elevation scans miss information at lower altitudes, and broader altitude bands are sampled at larger distances from the radar. As a result, altitudinal variability in the density of migrants can give rise to apparent spatial patterns in individual elevation scans, e.g., high altitude migration layers appear as spatial donuts when projected on the Earth (Figure 1). In many studies of migration using weather radar, the loss of information at lower altitudes with increasing distance from the radar is avoided by analyzing data within a relatively small radius around the radar [13,14]. Even though most commonly used software includes functionality to visualize radar data using “plan position indicator” (PPI) plots, these visualize one elevation scan and do not correct for a difference in observation height or height distribution with increasing distance from the radar. Furthermore, migration at higher altitudes is missed close to the radar in lower elevation scans. These problems are not unique to aeroecology; in radar meteorology, the vertical variability of precipitation can cause apparent horizontal structures in radar imagery. This type of pattern is often seen when the zero degree isotherm is at the height of the radar beam and is known as the bright band, caused by melting snow [15,16].

Figure 1 illustrates the challenges arising when trying to interpret a two-dimensional representation of migration density from one elevation scan when migration takes place at a high altitude. On this night, bird migration was concentrated in a layer between 2 and 3 km a.s.l. (see the vertical profile in Figure 1A). Elevation scans of such layering events typically show a donut-shaped reflectivity pattern (Figure 1B,C), with high reflectivities appearing at locations where the radar beam intersects with the high-altitude migration layer. This donut-shaped pattern poorly represents the spatial distribution in the number of migrants aloft.

Earlier work of Buler and Diehl [17] investigated corrections to radar imagery for these distance effects. Their work focused on the distribution of birds at departure after sunset [17,18], to identify areas in the landscape where migratory birds stop to rest and refuel. Flight altitudes are still low at this time, with birds observed primarily in the lowest elevation scans of the radar. The correction method therefore only needs to be applied to the lowest elevation angle scan.

For studies focusing on bird migration en-route, birds fly at multiple altitude layers [19], and a corrected image cannot be obtained from the lowest elevation scans alone. We, therefore, develop an approach that vertically integrates radar reflectivity over all elevation scans to estimate a spatially explicit image of radar reflectivity that accounts for changes in beam shape with distance from the radar. The vertical integration of elevation scans that we present can be applied across radars to get a detailed spatial representation of radar reflectivity from local to regional scales. Our goal is to get more insight into bird movements by providing a two-dimensional representation of radar reflectivity independent of distance. We apply our method to create a spatial distribution of the estimated vertically integrated density (VID) for 20 radars, which we display as a single composite. This composite visualizes the VID over regional and continental scales. We use this composite to explore dynamic visualizations and show its value in interpreting the spatial distribution of the VID and highlight multiple migratory phenomena. We also discuss some of the limitations of this approach and suggest areas for future research.

2. Materials and Methods

Here, we first explain the calculation of the VID map. To demonstrate the value of the method, we discuss observed patterns of several case studies described in more detail below. We emphasize the general application of the algorithm without radar specific optimizations.

2.1. Spatial estimates of Vertically Integrated Density

Our goal is to account for beam shape effects and to estimate a map revealing the spatial distribution of the number of birds aloft. To do so, we will:

• integrate over the vertical dimension by taking into account all radar scans recorded at different elevations
• account for the changing overlap between the radar beams as a function of range

In our analysis, we will treat the (normalized) altitudinal distribution and the number of birds at a radar location separately, because of different expected spatial auto-correlations for these two quantities. Wind and temperature patterns strongly govern the altitudinal distribution of flying birds (e.g., [19,20,21]). Since atmospheric circulation patterns vary relatively gradually geographically [22], we expect birds within the radar domain to select fairly similar flight altitudes [but note the possibility of substantial differences in flight altitude at water-land interfaces, e.g., [23] or Archibald et al. [24] on the descent of nocturnal migrants at dawn over the Great Lakes]. We will, therefore, treat the normalized altitudinal distribution of birds as constant within the radar domain. In contrast, we expect more geographic variation in the total numbers of birds aloft, as a result of variation in the take-off distribution, which reflects variation in land cover and available stopover habitat [17]. Additional spatial variation will be caused by the geometry of landmasses and coastlines that structure the migratory flow.

We estimate the normalized altitudinal distribution of birds using a vertical profiling method to distinguish birds from other sources of reflectivity and estimate density with altitude [13]. This provides us with the average reflectivity by birds $\eta \left(h\right)$ as a function of height h. The vertical profile is estimated only from data collected close to the radar (between 5 and 25 km range), since at longer ranges, the radar beam becomes too broad to resolve altitudinal patterns and may start to overshoot the layer of migration. We apply an often used algorithm; however, our approach is not specific to the vertical profiling method and is generalizable to other profiling methods, e.g., [14,25].

From the profile $\eta \left(h\right)$, we calculate the vertically integrated reflectivity ($\mathrm{vir}$) [26] as:

$$\mathrm{vir}={\int }_0^{h_{max}}\eta \left(h\right)dh$$

which gives us a measure that is proportional to the total number of birds in a vertical column for the area around the radar for which the profile was calculated. $h_{max}$ is the maximum altitude that will be considered in the analysis, in our case 5 km.

Given a geographic grid with grid cell centers $(x,y)$ in the radar domain, we calculate for each elevation scan i the polar coordinate $(r_i(x,y),\alpha (x,y))$ that defines the point in space vertically above location $(x,y)$ at height $h_i$, with $r_i(x,y)$ range and $\alpha (x,y)$ azimuth. Geographic coordinates $(x,y)$ are converted to polar coordinates $(r_i,\alpha )$ using a spheroid Earth model with an equatorial radius of 6378 km and a polar radius of 6357 km, while accounting for beam refraction in the atmosphere using the 4/3 effective Earth’s radius model [27]. Here, we use a standard model for atmospheric propagation. The work of Buler and Diehl [17] accounted for changes due to non-standard refraction effects using radiosonde data. They noticed only minor improvements in the accuracy of their final model. For this reason and the complications of including additional separate sources of information, we did not include this additional correction. The reflectivity (in linear units (cm²/km³)) measured at this location in scan i we denote as ${\eta }_i(x,y)$ [13,28].

We assumed that the antenna pattern $b(h,r,{\theta }_i)$ had a Gaussian shape, see [29] and could hence be expressed as a function of height h, range r, and beam elevation: ${\theta }_i$ as:

$$b(h,r,{\theta }_i)=\frac{1}{\sqrt{2\pi }}{e}^{-\frac{1}{2}{\left(\frac{h-h_b(r,{\theta }_i)}{w_b\left(r\right)}\right)}^2},$$

with $h_b(r,{\theta }_i)$ the height of the center of the beam (accounting for the radar’s antenna height) and $w_b\left(r\right)$ the beam’s width expressed as a normal standard deviation. Heights are expressed in meters above mean sea level throughout.

$b(h,r,{\theta }_i)$ is plotted in Figure 2 for a volume scan of the Herwijnen radar in the Netherlands, with different elevations ${\theta }_i$ indicated by different colors, and the radiated energy for four different ranges (20 km, 50 km, 75 km, and 100 km). The profiles of radiated energy broaden with range, and the overlap between the lower beams is considerable for heights below 5 km at longer ranges.

For each grid location $(x,y)$, we calculate the reflectivity $\widetilde{{\eta }_i}(x,y)$ that we expect to measure if the altitude distribution of birds aloft at this location would be identical to the vertical profile $\eta \left(h\right)$ measured at close range to the radar:

$$\widetilde{{\eta }_i}(x,y)={\int }_0^{h_{max}}b[h,r_i(x,y),{\theta }_i(x,y)]\eta \left(h\right)dh$$

Finally, we obtain a spatial image of the adjusted vertically integrated reflectivity (${\mathrm{VIR}}_a$), as follows:

$${\mathrm{VIR}}_a(x,y)=\mathrm{vir}\times R(x,y),\mathrm{with} R=\frac{{\sum }_i{\eta }_i(x,y)}{{\sum }_i\widetilde{{\eta }_i}(x,y)}$$

Our spatial image of ${\mathrm{VIR}}_a$ equals the $\mathrm{vir}$ of the profile, times a spatial adjustment factor R. R is the ratio between the total sum of observed reflectivity values ${\eta }_i$ and the total sum of expected reflectivity values ${\tilde{\eta }}_i$ (i.e., expected given the bird’s vertical profile $\eta \left(h\right)$ and the radar’s beam shape $b(h,r,{\theta }_i)$ for each elevation scan at that location). If higher reflectivities are observed at $(x,y)$ than expected based on the average profile $\eta \left(h\right)$, then the ratio $R>1$, and the $\mathrm{vir}$ value is adjusted upward. Likewise, if a lower reflectivity sum is detected at $(x,y)$, then ratio $R<1$, and the $\mathrm{vir}$ value is adjusted downward. The combined effect of this operation is a correction for changing overlap between the radar beams and the layer of birds and an integration over the altitudinal dimension.

In order to present biologically meaningful data, we converted integrated reflectivity to bird density by assuming a radar cross-section of 11 cm² [13] for a single bird. This size is representative for passerine birds, which are the dominant nocturnal migrants in this region at this time of year.

An implementation of this approach is readily available within the biorad R package [26,30], through the function integrate_to_ppi.

2.2. Data

We illustrate the method using two case studies of composite radar images. We obtained data directly from the respective national weather services, of the Netherlands, Belgium, and Germany. All radars operated at C-band with a wavelength between 5.25 and 5.35 cm. The polar volume data (named “level II data” in the USA) was converted to the ODIM format where needed [31] and processed into vertical profiles, based on data within the 5–25 km range from the radar [13]. The calculation of vertical profiles used several filtering steps to identify birds, for example filtering out rain by a threshold on the value of the correlation coefficient moment.

There are a few consistent differences in the volume coverage pattern of radars between countries. The radars in the Netherlands have a lower first elevation angle of 0.3°, while both Germany and Belgium have a lowest elevation angle of 0.5°. The Dutch radars furthermore conduct three different scans within a five minute volume at this lowest elevation, two of which have a lower pulse repetition frequency, allowing the radar to look further. Our method accounted for these differences in volume coverage pattern as we calculated the vertically integrated reflectivity based on all available scans.

2.3. Case Studies

The first case study, a dataset of four radars across the Netherlands and Belgium, covered one night (3–4 October 2016) with intense migration across Europe [5]. In these countries, radars produced a full volume scan every five minutes, allowing for detailed interpretation of fine-scale temporal patterns. Sunset was at 17:07 (UTC; all times are UTC), while sunrise was at 05:49. To visualize these data, we made a composite by using estimates of bird densities from the closest radar. We also showed Voronoi polygons around the radar locations and a buffer of 200 km around each radar to indicate which radar was responsible for which area. The polygons ensured we always visualized data from the closest radar. For comparison across images, we used a fixed color scale. We visualized values beyond the color scale as the most extreme value of the color scale. When we could not calculate a vertical profile for biological echoes, for example, because of temporary rain close to a radar, we could not calculate the VID image, and therefore, the data of that radar are not shown.

A second case study is presented for a larger region that includes Germany. We generated a similar visualization to the first case study, for the night of 18 October 2017 at 21:00 that included in total 20 different radars. Exploration of the German data revealed that strong filtering was applied to the processed reflectivity data (quantity DBZH). After this filtering, hardly any biological signal was retained especially further away from the radar. Therefore, we opted to use the total unfiltered reflectivity, quantity: TH; see [31] and removed the part of the signal that was due to static objects as determined by a Doppler notch filter (quantity CCORH; [27]).

3. Results

The effect of vertical integration was evident from the result that the donut shape visible in individual scans (Figure 1B,C) was no longer visible in the VID product (Figure 3). The spatial VID estimate highlighted a region of high migratory density in the southeast of the Netherlands that was only visible in a limited number of scans as for most scans, the radar beam was either too low or too high (e.g., Figure 1C).

In the first case study, many features of migration were visible in the animations of the VID (Figure 4, Supplement S1, [32]). As the height profile varied over the night (Supplement S2), many of these features would either not be visible or more difficult to identify in normal PPI plots of the lowest elevation scan. Just before sunrise, very weak migration was visible (Figure 4A); bird densities were very low (speckled pattern around each radar), with very high reflectivity representing non-biological features (see below). Migration started between 17:40 and 18:00, and VID increased across the whole region surveyed by radars. Immediately thereafter, estimated densities above the sea off the coast of Belgium and the Netherlands also considerably increased as birds flew in the west-southwest direction. Besides these large-scale patterns, more local features were visible, for example, in the northern regions of the Netherlands where the density of birds considerably increased from 22:30 onwards (Figure 4C). This high density of migrants then moved across the IJsselmeer (the large lake in the northwest of the Netherlands) towards the west, seemingly avoiding the longest lake crossings. These waves of high density propagated across the borders between radars, showing that spatio-temporal patterns in migratory density were an actual phenomenon that could be observed in different radars. After sunrise at 5:10, most migration over land stopped. At sea, off the coast of the west of the Netherlands, a dawn ascent was visible starting around 5:30 (Figure 4D), a phenomenon resulting from birds quickly increasing their flight altitude at dawn. We recognized this ascent by the simultaneous increase in bird density across regions above the sea in consecutive time slices. At the same time, there was an increase in daytime movements above land presumably due to migratory take-off of diurnally migrating species.

Non-biological features were also visible in the animations. By combining knowledge from online resources, these could be attributed to several sources. Several rainfall events were visible as patches of very high reflectivity moving with the wind from northeast to southwest and early in the morning turning to northwest (Figure 4A). Off the northwest coast of the Netherlands, shipping lanes were recognizable by a series of moving spots of high reflectivity (Figure 4C). Wind farms, both at sea and land, and oil platforms at sea were recognized as restricted static regions of high reflectivity (Figure 4B). Other more detailed dynamic features were, for example, the smoke of a fire in an industrial building southeast of the Herwijnen radar in the Netherlands between 23:00 and 00:00 (Figure 4C, [33]).

When visualizing over a larger spatial region in the second case study (Figure 5), other phenomena became visible. The reflectivity as caused by birds taking off increased first in the northeast, where the Sun set about 35 min earlier compared to the western part of the visualization. Over the next half hour, the reflectivity increased gradually and spread to the west (Supplement S3, [34]). Once the sun set over the whole region, a fully joint picture of the high bird densities across the region was observed. The resulting image was spatially relatively uniform, despite operational differences between radars, such as differences in elevation angles. Some radars showed an apparent higher bird density on the axis perpendicular to the main migratory direction from northeast to southwest. This was an aspect effect, related to the birds having a higher radar cross-section in side view compared to the frontal or rear view [35]. The radar in the northeast of Germany, close to Berlin, stood out as having a high reflectivity. A detailed inspection of individual PPI plots indicated that this radar was affected by anomalous atmospheric propagation. The fact that a Doppler filtered reflectivity product (DBZH) was not available for this radar for biological echoes further complicated the removal of the resulting ground clutter. Towards the Alps in southern Germany, radars also seemed to have an overestimation of densities at larger distances from the radar.

4. Discussion

The spatial patterns of aerial bird density that we presented reflect the effects of habitat, coastlines, and upstream source areas. These effects are frequently more fine-scale than the distances between radars and would not be identified when aggregating data to a single point of information for each radar. Here, we show that visualizations of vertically integrated range-corrected reflectivity measured across weather radars can be used on a regional scale to gain insights into the distribution of birds and their movement. The high density of the radar network in Europe is advantageous for showing this continuous overview, compared to the sparser network in North America. Our method effectively integrates over the vertical dimension, while accounting for the different scanning angles of each radar station. Many of the features we identify are much harder, if not impossible, to identify if we had not estimated VID images to account for the distance effect or when not using animations of sequences of composites. By removing range effects, it is possible to focus on biological phenomena and avoid the distraction of non-uniform sampling. This work thus makes it possible to generate and test hypotheses on the spatial-temporal distributions of migratory behaviors of all migrating birds. One such phenomenon is the increased density we see in Figure 4C; this can be caused by improved weather conditions northwards of the observed area. Alternatively, birds that depart synchronized at sunset from a stopover region and move in a common direction can have a similar effect (e.g., coastal areas before a sea crossing). Animations over larger spatial regions can help to disentangle these hypotheses. Better understanding of local movement patterns is important for reducing human-wildlife conflicts, for example, when planning wind farms. Detailed investigation of spatial distributions of birds in the air can help to quantify the influence of landscape geometry on local aggregations of birds that should be avoided to reduce impact.

Differences in the radar cross-section of individual birds related to body orientation reduce the accuracy of the algorithm. This problem will be especially relevant when birds have a common orientation, as is the case during migration. The magnitude of this error will not exceed the magnitude of differences in radar cross-section between orientation (e.g., a factor of two, [36]). Robust and automated methods for estimating animal’s body orientation are needed in order to apply spatial corrections for such aspect effects, e.g., based on patterns in dual-polarization moments [37], which would be a valuable future improvement.

At large distances where even the lowest elevation scans start to overshoot the migration layer, both the observed and expected $\eta$ values become increasingly small, which will compromise the accuracy of our methods. At larger distances, weak signals of migratory animals may not reach above the lower detection limit of the radar. In our visualizations of northwest Europe, the average nearest neighbor distance between radars is 125 km; radars are thus frequently only visualized up to half this distance, where such limits on detection are typically not reached.

We assume the vertical distribution of birds to be constant within the field of view of a single radar; however, there might be cases where this assumption is not valid. Vertical distributions of birds are likely driven by weather patterns [20], and weather conditions are fairly smooth over the range of a single radar, except perhaps in mountainous terrain or along coastlines [38]. Hence, our assumption is expected to be valid when meteorological conditions, which strongly determine the vertical profile, are correlated over the domain of a radar. Such conditions will be more common in relatively homogenous topography and in the absence of direct interaction with precipitation, which might force birds abruptly down. Vertical distribution of birds can also differ between land masses and large water bodies, e.g., [23,39]. Vertical profiles are derived close to the radar and thus mostly over land and will therefore be more adequate for migration over land. It is possible to account for these heterogeneous altitude distributions between land and sea by estimating a separate vertical profile above the sea if the radar is sufficiently close to water. Furthermore, birds might follow the contours of the landscape and fly at a constant altitude above ground level, while in other cases, birds may maintain a constant altitude above mean sea level regardless of the underlying topography, e.g., [13,19,40]. Since the elevation of the terrain can change considerably over the area of observation, this might create inaccuracies in the VID estimate. Here, we have mostly shown applications in a relatively flat environment where topography effects are expected to be limited. However, in southern Germany, radars are placed within mountainous regions or at least cover mountains in their respective domains. In these radars, the VID image seems to capture the spatial distribution of birds less accurately. Here, the vertical profile shows high bird densities close to the ground, while further, towards the Alps, birds might have been forced up by the topography resulting in an overestimation by the VID. Furthermore, radars might suffer from (partial) beam blockage in mountainous areas, e.g., [41], resulting in underestimates of bird densities. Additionally, the light emitted from cities at night can influence the altitudinal distribution of birds moving during the night. A recent study investigating this effect found a relatively small effect meaning that these changes will not have a strong effect on our VID estimates [42].

In recent years, smaller, more mobile dedicated bird radars have become more common [43]. These provide insight into the local vertical distribution of migratory birds, up to a limited height. Future research using these radars will provide more insight into the spatial variability of the vertical distribution, by providing more insight into the effects of elevation or habitat on the vertical distribution. This work paves the way toward developing methods that make more accurate assumptions on the local vertical profile compared to the nearest neighbor interpolation (taking the vertical profile of the closest radar) that is currently used; for example, by taking into account the full auto-correlation structure between profiles [44]. Combining weather radar with local bird radars is a powerful combination in which weather radar can provide the context for bird radars. Alternatively, weather radars can identify interesting phenomena that warrant a more detailed investigation using bird radars.

By inspecting the areas at the border between two radars, the accumulation of all possible inaccuracies can be investigated. Our animations show that much of the temporal variation is consistent between radars and that events of high bird densities are visible across radars. In the animation of the Netherlands and Belgium, the case of 3–4 October 2016; see Figure 4 and Supplement S1; [32], we see twice that a high-density area arises in the northern part of the Netherlands (around 18:45 and 23:30) and can be observed later on the radar in the center of the Netherlands. Similarly, the presence of redundant information in the area of overlap between two radars offers opportunities to improve accuracy locally. There are however consistent differences between radars. The radars in Belgium, for example, seem to have consistently shorter detection ranges compared to Dutch radars, which is likely related to the different acquisition setting of the respective radar processors.

To conclude, our method for estimating spatial VID images will enable the exploration and qualitative analysis of migration density within the field of view of a single radar. The method also facilitates the creation of composite maps of vertically integrated density similar to meteorological composites [10]. These provide a powerful tool for experts to gain insight into a range of migratory phenomena. Future improvements are approaches that can account for regional variations in altitudinal distributions [44,45], as well as tools to exclude non-biological signals to ease the interpretation and work toward fully quantitative analysis. Lin et al. [46] already took important initial steps by segmenting out rain from the projected images using machine learning; alternatively, dual-polarization metrics can be used, e.g., [47]. We see an important future for investigating regional variation in migratory patterns using weather radar.