P2: Analysis and visualisation of vertical bird profiles

In this practical you will learn to compute, interpret and analyze vertical profiles of birds (VPBs). A VPB consists of the bird density, speed and directions at different altitudes at the location of a single radar. It is typically calculated for all data within a cilinder with a 25 km radius around the radar, i.e. only containing data at relatively short distances, where the radar beam is still sufficiently narrow to resolve altitude information.

1. Preparations

Start again by making a new directory on your local machine that you will use for this practical

# make a new local directory on your machine where to download data for this practical
# replace the string below with the path of that directory.
HOME="your/personal/working/directory/"
# check that the directory exists. If the next statement evaluates to FALSE, something went wrong: the directory does not exist or you didn't specify its path correctly
file.exists(HOME)
# we will make HOME our work directory, the default folder in which R will look
# for files, and where it will output newly generated files.
setwd(HOME)
# Finally, we set the local time zone to UTC, so all plotted time axes will be in UTC
Sys.setenv(TZ="UTC")
# load the bioRad package
library(bioRad)

Your R session is now properly set up

2. Processing polar volume data into vertical profiles: vol2bird algorithm

The following steps take you through the process applying the vol2bird algorithm yourself. You need a working installation of Docker (linux / mac) or Docker for Windows (not Docker Toolbox, this is an older implementation of Docker for Windows operating systems that is not supported).

We will generate vertical profiles with the automated algorithm vol2bird (https://github.com/adokter/vol2bird), which is included in the bioRad package

# start your local Docker installation
# we first test whether R can communicate with Docker:
checkDocker()

If you get a “Hello from Docker!” welcome message, everything is working and you can start processing

# We will continue with the Swedish data you downloaded during the first practical.
# first, make a list containing all the polar volume files starting with 'seang_pvol_20150927'
files=list.files(pattern="seang_pvol_20150927",full.names=TRUE)
# second, make sure we only use the files that end with 'h5', so hdf5 files with a h5 extension only
files=files[grep("h5$",files)]
# let's process the first file, which has the name 'seang_pvol_20150927T030000Z.h5'
file.in=files[1]
# run vol2bird
vp=vol2bird(file.in)
# vp is now a 'vp' profile object, that you can examine as in the previous exercises
vp
#
# alternatively, you may also store the profile as a hdf5 file, which is what we will do next:
# let's first define the name of the output file (we replace the substring 'pvol' into 'vp')
file.out=gsub("pvol","vp",file.in)
# print the newly generated output file name to which we will write:
# note that this is only a string, you can give it any other name if you want to
file.out
# finally, run vol2bird; this generates an output file as specified in file.out
# we set autoconf to TRUE, to let vol2bird figure out the optimal settings by itself
vol2bird(file.in,file.out,autoconf=TRUE)
# your work directory should now contain a new file 'seang_vp_20150927T030000Z.h5'
# check that we can read this file, and retrieve the vertical profile from it:
vp=readvp(file.out)

3. Loading processed vertical profiles

To save some time, we will continue with some pre-processed vertical profiles generated for the case study you worked on yesterday (Angelholm radar in Sweden).

Navigate to https://www.dropbox.com/sh/8dl30htd9zxdbx3/AAAMOHk4_UlCsXOYb6LDFONOa?dl=0 and download the processed profiles, profiles.zip.

# unzip the profiles.zip file, either by clicking the file or the command below; after unzipping you should have a folder '2015' in your working directory with 4 days of profile data (20,21,26 and 27 September 2015).
unzip("./profiles.zip")
# load all the filenames in the new 2015/09 subfolder in your working directory
vp_paths=dir("./2015/09", recursive=TRUE,full.names=TRUE)
# let's make a subselection of files: only those that have '20150926' or '20150927' in the filename, i.e. profiles from 26 sep 2015 to 27 sep 2015.
vp_paths_sep26=vp_paths[grep("20150926",vp_paths)]
vp_paths_sep27=vp_paths[grep("20150927",vp_paths)]
# paste the two lists together:
vp_paths_sep2627=c(vp_paths_sep26,vp_paths_sep27)
# read these vertical profiles (hdf5 files) into R (may take 1-2 minutes to load)
vplist=readvp.list(vp_paths_sep2627)
# print some information on the vplist object. It should contain 576 profiles 
vplist
# save the object, which allows you to load the data more quickly next time
save(vplist,file="vplist_sep2627.RData")
# you can restore the vplist object at any time as follows:
load("vplist_sep2627.RData")

4. Inspecting single vertical profiles

Now that you have loaded a list of vertical profiles, we can start exploring them. We will start with plotting and inspecting single vertical profiles, i.e. a single profile from the vplist object you have just loaded.

# let's extract a profile from the list, in this example the 25th profile:
vp=vplist[25]
# print some info for this profile to the console
vp
# test whether this profile was collected at day time:
day(vp)
# plot the vertical profile, in terms of reflectivity factor
plot(vp, quantity="dbz")
# plot the vertical profile, in terms of reflectivity
plot(vp, quantity="eta")

These two plots look very different, but they are twice the same data plotted on a different scale.

eta = (radar-wavelength dependent constant) * 10^(dbz/10)

So eta and dbz are closely related, the main difference is that reflectivity factors are logarithmic, and reflectivities linear.

The reflectivity factor dbz is the quantity used by most meteorologist. It has the useful property that at different radar wavelengths (e.g. S-band versus C-band) the same amount of precipitation shows up at equal reflectivity factors. The same holds for insects, as well as any other target that is much smaller than the radar wavelength (S-band = 10 cm, C-band = 5 cm), the so-called Rayleigh-scattering limit

In the case of birds we are outside the Rayleigh limit, because birds are of similar size as the radar wavelength. In this case the reflectivity eta is more similar between S-band and C-band (though might be a little different). eta is also more directly related to the density of birds, since eta can be expressed as (bird density) x (radar cross section per bird). For these two reasons, for weather radar ornithologists reflectivity eta is the more conventional unit (instead of reflectivity factor in dBZ).

# plot the vertical profile, in terms of bird density
plot(vp, quantity="dens")
# print the currently assumed radar cross section (RCS) per bird:
rcs(vp)
  1. If you change your assumption on the bird’s radar cross section in the previous example, and assume it is 10 times as large, what will be the effect on the bird density profile?

The assumed radar cross section can be changed as follows:

# let's change the RCS to 110 cm^2
rcs(vp) = 110
  1. Verify your answers on the previous two questions, by re-plotting the vertical profiles for the bird density quantity.

5. Plotting time series data

We will now examine multiple vertical profiles at once that are ordered into a time series, e.g. the vertical profiles obtained from a single radar over a full day.

# convert the list of vertical profiles into a time series:
ts=vpts(vplist)
# print summary information
ts
# extract the position of the radar from the metadata:
lon=ts$attributes$where$lon
lat=ts$attributes$where$lat
# make a subselection for a certain time range, e.g. the night of 26/27
# september between sunset and 7 UTC. We use the suntime function to calculate sunset
t.start=suntime(lon,lat,"2015-09-26",rise=F) # gives moment of sunset
t.end=as.POSIXct("2015-09-27 07:00", tz="UTC")
# now that we have the start and end dates defined, we can make a subselection:
ts.night=ts[ts$dates>t.start & ts$dates<t.end]
# plot the time series
plot(ts.night)
# check the help file for the plotting function of profile time series
# Because profile timeseries are of class 'vpts', it's associated plotting function
# is called plot.vpts:
?plot.vpts
  1. During practical 1 you made a series of PPI images for this same case (from 3 UTC - 8 UTC). Scroll through yesterday’s images once more with the vertical profile time series next to it. Birds ascend and depart in two different cohorts shortly after each other. Make sure you recognize both in the PPIs and in the profile data.
  2. Plot the time series also for reflectivity (quantity eta) and reflectivity factor (quantity dbz). Use the help documentation to check out the plotting options available.
  3. Interpret the wind barbs in the profile time series figure: what is the approximate speed and direction at 1400 metre at 22 UTC. In the speed barbs, each half flag represents 2.5 m/s, each full flag 5 m/s, each pennant (triangle) 25 m/s.
# time series objects can be subsetted, just as you may be used to with vectors
# here we subset the first 50 timesteps:
ts[1:50]
# here we extract a single timestep, which gives you back a vertical profile class object:
ts[100]
# to extract all the dates of all profiles in the time series:
ts$dates
  1. Extract the vertical profile at 22 UTC from the time series and plot the vertical profile of ground speeds (quantity ff, see manual for the vp class for full list of available quantities). Check whether your answer to the previous question was approximately correct.

6. Vertical integration: surface density & migration traffic rate

Often you will want to sum together all the migrants in the vertical dimension, for example if you want a single index of how many birds are migrating at a certain instant. There are at least two ways in which you can do that

  • by calculating the vertically integrated bird density (VID), which is surface density as opposed to a volume densities you have been plotting in the previous exercises: this number gives you how many migrants are aloft per square kilometer earth’s surface (unit individuals/km\(^{2}\)), and is a vertical integration of the volume densities (unit individuals/km\(^{3}\)).
  • Note that the VID quantity doesn’t depend on the speed of the migrants. A common measure that reflects both the density and speed of the migration is migration traffic rate (MTR). This is flux measure that gives you how many migrants are passing the radar station per unit of time and per unit of distance perpendicular to the migratory direction (unit individuals/km/hour).

We will be using bioRad’s vintegrate function to calculate these quantities

# Let's continue with the ts object created in the previous example.
# The vertically integrated quantities are calculated as follows:
vintegrated.ts = vintegrate(ts)
# The vintegrated.ts object you created is a vivp class object, which is an acronym for Vertically Integrated Vertical Profile. It has its own plot method, which by default plots migration traffic rate (MTR):
plot(vintegrated.ts)
# you can also plot vertically integrated densities (VID):
plot(vintegrated.ts, quantity="vid")
# the gray and white shading indicates day and night, which is calculated 
# from the date and the radar position. You can also turn this off:
plot(vintegrated.ts, nightshade = FALSE)
# execute `?plot.vivp` to open the help page listing all the options.
?plot.vivp

The following questions only require pen and paper. Assume a hypothetical situation in which the volume density of birds from 0-1000 metre above ground is 200 birds per cubic kilometer, and from 1000-1500 metre 100 birds per cubic kilometer. In the lower layer birds fly at 50 km/hour, and in the upper layer at 100 km/hour. Above 1500 metre there are no birds.

  1. What is in this case the bird’s surface density (or as we previously called it, vertically integrated density VID)? Give your answer in units birds/km\(^2\)
  2. What is in this case the migration traffic rate? Give your answer in units birds/km/hour

Both MTR and VID depend on the assumed radar cross section (RCS) per bird. If you are unwilling to make any assumptions on RCS, you can alternively use two closely related quantities that do not depend on RCS:

# print the currently assumed radar cross section:
rcs(vintegrated.ts)
# instead of VID, you can use vertically integrated reflectivity (VIR):
plot(vintegrated.ts, quantity="vir")
# instead of MTR, you can use the reflectivity traffic rate (RTR):
plot(vintegrated.ts, quantity="rtr")

VIR gives you the total cross-sectional area of air-borne targets per square kilometer of ground surface, whereas RTR gives you the total cross-sectional area of targets flying across a one kilometer line perpendicalur to the migratory flow per hour.

  1. Re-plot the vertically integrated quantities discussed so far using a different radar cross section, and verify which quantities are affected by a change in radar cross section.

7. Inspecting precipitation signals in profiles

Precipitation is known to have a major influence on the timing and intensity of migration, therefore it is a useful skill to be able to inspect profiles for presence of precipitation.

Also, although automated bird quantification algorithms become more and more reliable, distinguishing precipitation from birds remains challenging for algoriths in specific cases. It is therefore important to have the skills to inspect suspicious profiles. That may help you to identify potential errors of the automated methods, and prevent your from overinterpreting the data.

An easy way of doing that is plotting the vertical profile of total reflectivity (quantity DBZH), which includes everything: birds, insects and precipitation. Precipitation often has higher reflectivities than birds, and also extends to much higher altitudes.

  1. Load data from 20 Sep 2015 12:00 UTC to 21 Sep 2015 12:00 UTC into a time series object (just as you did earlier for two other days, see last code block section. Make a profile plot of bird for this night (see section 5), as well a plot of the migration traffic rate throughout the night (see section 6).

You should see a strong drop in migration traffic rates around 23 UTC, about halfway through the night, which is relatively early. We will check whether precipitation may have played a role.

  1. Compare your profile plot for bird density (quantity ‘dens’) with a profile plot for total reflectivity (quantity DBZH, showing birds and precipitation combined). Compare the two plots to visually identify periods and altitude layers with precipitation. Make a qualitative guess how rain may have affected this migration event.

8. Processing multiple polar volumes

This section contains an example for processing a directory of polar volumes into profiles

# read the filenames of the polar volumes you want to process
files=list.files("~/tmp/nexrad/",full.names=TRUE)
# print the filenames
files
# create output directory for processed profiles
outputdir="~/tmp/nexrad_out"
dir.create(outputdir)
# let's loop over the files and generate profiles
for(fileIn in files){
  # generate the output filename for the input file
  fileOut=paste(outputdir,"/",basename(fileIn),"_vp.h5",sep="")
  # generate the profile, and write it as a hdf5 file (for later reference)
  # we turn autoconfiguration on, so vol2bird chooses the optimal settings for the file automatically
  vp=vol2bird(fileIn,fileOut,autoconf=T)
}

Having generated the profiles, we can read them into R

# we assume outputdir contains the path to the directory with processed profiles
profileFiles=list.files(outputdir,full.names=TRUE)
# print them
profileFiles
# read them
vplist=readvp.list(profileFiles)

You can now continue with visualizing and post-processing as in section 5 and onwards

# make a time series of profiles:
ts=vpts(vplist)
# plot them between 0 - 3 km altitude:
plot(ts,ylim=c(0,3000))