P2.1: Basic analysis and visualisation of vertical bird profiles

1. Preparation

Radar data needed for this practical can be downloaded from Dropbox. Don’t download everything right now, but when you are asked to in the exercises (so we spread our Wifi use). https://www.dropbox.com/sh/ajwdo578jfitlsk/AAC6FyCUNGq34K4QdhfQd4wea?dl=0

Throughout the practical we will use the bioRad R package. It’s currently not (yet) in the default CRAN repository, but can be downloaded from Github. See https://github.com/adokter/bioRad for install instructions. Go through steps 1-3 of these install instructions if you haven’t done so already. Execute each of the code examples in Rstudio, and try to complete the exercises.

# make sure you start with a fresh R session
# load the bioRad package
library(bioRad)
# check the package version
packageVersion("bioRad")
# make sure you have the latest version (0.2.1). If you have an older version, download and reinstall it by:
library(devtools)
install_github("adokter/bioRad")

The first half of this practical deals mostly with vertical profile data of birds, i.e. the bird speed and directions at different altitudes. These profiles are extracted from polar volume data by an automated algorithm (vol2bird). In the second half we will have a look at the underlying polar volume data as well.

All the functions of the bioRad package have fairly elaborated manual pages, that you can refer to at any moment:

# bring up the package general help page:
?bioRad

In the help window, you can click on each of the links to pull up help pages of specific functions. Click on the Index link at the bottom of each help page to get an overview of all functions and available manual pages.

Start 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")

Your R session is now properly set up

2. Loading radar data into R

There are two ways to obtain bird profile data: - ENRAM repository: browse to to see what is available. We will start this practical with this approach. - Process from polar volume data into profiles yourself: in this case you need to have access to polar volume data (also known as level 2 data), the lowest level radar data provided by meteorologists. We will leave this as an exercise for at the end of the practical (section 8).

  1. Open your browser and go to the ENRAM repository, and download a zip file with profiles.

We can make this process a little bit easier with the download_vp command:

# Let's download some data for Sweden (country code 'se'), for the Angelholm radar (radar code 'ang'), for the month of January 2017:
download_vp("2017-01-01", "2017-01-31", "se", "ang", localpath = HOME)
  1. Use your favorite file explorer to inspect your HOME directory. It should now list a zip file, which has been unzipped into a directory structure organised as country/radar/year/month/day/hour/

Let’s load some of thse vertical profiles into R:

# first we retrieve the paths of the vertical profile files, say for 14-17 January:
vp_paths=retrieve_vp_paths(HOME, "2017-01-14", "2017-01-17", radar="ang")
# print the list of files:
vp_paths
# read the list of files:
vplist=readvp.list(vp_paths)
# print the vplist object:
vplist

As you might have noticed, the ENRAM repository is still quite empty, but we hope to be ready for upcoming spring.

Instead, we will use vertical profiles generated for the ENRAM validation campaign at Kullaberg in Sweden in 2015. We will downloaded the data from a different location:

# download data
download.file("https://github.com/adokter/ENRAM-training/raw/master/data/seang20150920-30.zip", destfile = "seang20150920-30.zip")
# unzip the data; after unzipping you should have a folder '2015' with all the data
unzip("seang20150920-30.zip")
# load all the filenames
vp_paths=dir(".", recursive=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)]
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")

3. 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 250'st profile:
vp=vplist[250]
# 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 quantity eta is more similar between S-band and C-band (though might be a little different)

# 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?
  2. Same question, but for the reflectivity 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 reflectivity and bird density quantities.

4. 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)
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)
# open the help file for the plotting function.
# Because profile timeseries are of class 'vpts', it's associated plotting function
# is called plot.vpts:
?plot.vpts
  1. 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.
  2. Interpret the wind barbs in the 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 timesteps do:
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.

5. 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.

P2.1: Advanced interpretation of profile data

6. Inspecting precipitation signals

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. At C-band, precipitation typically 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 of the Angelholm radar of the kullaberg campaign into a time series object (see section 2). Make a plot of the migration traffic rate throughout the night (see section 5).

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. Make profile plots for this night (see section 4). Make a plot both for bird density (quantity ‘dens’) and total reflectivity (quantity DBZH, showing birds and precipitation combined). Compare the two plots to visually identify periods and altitude layers with precipitation, and make a qualitative guess how rain may have affected this migration event.

7. Relating profile data to underlying polar volume data

In this section we will inspect a migration event with combined precipitation and biological scattering in more detail, by examining polar volume data.

Before continueing with the next example, first download this volume file from Dropbox. Unzip the file and put the file in your working directory.

# load the polar volume data from the h5 file you just downloaded
pvol=read.pvol("P2.seang_pvol_20150920200000.h5")
# print some information about the polar volume
pvol
# print information about the polar scans in this polar volume:
pvol$scans
# let's extract the third scan, which was collected at 1.5 degree elevation:
pscan = pvol$scans[[3]]
# print some information about this scan:
pscan
# before we can plot the scan, we need to project it on a Cartesian grid,
# i.e. we need to make a Plan Position Indicator (PPI)
my_ppi = ppi(pscan)
# print some information about this ppi:
my_ppi
# you can see we projected it on a 500 metre grid.
# Check the manual of the ppi function to see how you can change the projection
# Now we are ready to plot the ppi
# plot the reflectivity factor image:
plot(my_ppi, param="DBZH")
  1. This event contains convective precipitation, characterised by localized but intense thunderstorms, as well as biological scattering. Make also a ppi plot of the radial velocity. How does the texture of the radial velocity of the precipitive areas differ from the areas with biological scattering?

  2. Estimate the absolute difference in ground speed of the precipitation and the ground speed of the biological scatterers by visually inspecting the radial velocity plot. Are the biological scatterers birds or insects? Why?

Let’s go back to the ts.night object you created in section 4 (re-execute the code there if you lost the object along the way).

  1. plot the migration traffic rates for the ts.night object (see section 5). What is your biological intepretation of the temporal pattern, specifically in relation to the transition between night and day?
  2. Make a vertical profile plot of the diurnal migration from 27 Sep 2015 4:00 UTC to 7:00 UTC (see section 4).

It should be clear from your last two figures that there is peak of nocturnal migration, followed by a peak of diurnal migration. The diurnal peak shows two subsequent departures of birds that coincide with sunrise. We are going to look deeper into the polar volume data from which the vertical profiles were calculated, to see what these two departures at sunrise might be.

Before continueing with the next example, first download these volume files from Dropbox. Unzip the file and put the newly extracted files in your working directory.

# load the file names of the files you just downloaded
# make sure the h5 files you just downloaded are in your working directory.
files=list.files(paste(HOME, "20150927", sep=""), pattern="seang_pvol_20150927",full.names=TRUE)
# print the file names
files
# we will plot the file at 4:30 UTC, when the first emersion is present, but the second not yet (check whether you came to approximately the same answer in the previous exercise)
pvol=read.pvol(files[4])
# print some information about the polar volume
pvol
# print information about the polar scans in this polar volume:
pvol$scans
# let's extract the 2 degree elevation scan, in this case the fifth scan
pscan = pvol$scans[[5]]
# print some information about this scan:
pscan
# before we can plot the scan, we need to project it on a Cartesian grid,
# i.e. we need to make a Plan Position Indicator (PPI)
my_ppi = ppi(pscan)
# print some information about this ppi:
my_ppi
# you can see we projected it on a 500 metre grid.
# Check the manual of the ppi function to see how you can change the projection
# Now we are ready to plot the ppi
# plot the reflectivity factor (we adjust the plotting scale to -20 to 15 dBZ):
plot(my_ppi, param="DBZH", zlim=c(-20,15))
# even more informative is plotting the data on a google earth map
# first download the background image:
satelliteImage=basemap(my_ppi,maptype="satellite")
# then overlay the PPI on the satellite image:
map(my_ppi,param="DBZH",satelliteImage, zlim=c(-20,15))
  1. Make a series of ppi maps for the 2.5 degree polar scan for all the volume files you downloaded for this migration event, and save them to disk. Cycle through the images and try to figure out if you can distinghuish the two emersions, and whether they have a different spatial origin. You can either do that by hand, or - if you feel confident with R - use below example code to generate and save the images automatically.
# let's store all the polar volumes we want to plot in a vector:
# first, get all the files 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 only
files=files[grep("h5$",files)]
# load the ggplot2 library, because we will use its ggsave function later on to save images
library(ggplot2)
# we will use a for-loop to process the images one by one
# this simple loop prints the contents of variable 'i' iterating over a numeric vector
for (i in c(1,2,4,8)){
  print(paste("in this iteration variable i equals",i))
}
# now let's make a more advanced for-loop that makes and save the images one by one:
for (file in files){
  # 1) read the polar volume 
  vol=read.pvol(file)
  # 2) extract the first scan
  scan=vol$scans[[5]]
  # 3) make a ppi
  ppi = ppi(scan)
  # 4) plot the ppi, and store it in the object 'myPlot'.
  myPlot=map(ppi,param="DBZH",satelliteImage,zlim=c(-20,15))
  # 5) make a string with the filename in which to store the image
  myPlotFilename=paste(file,".jpeg",sep="")
  # 6) print the filename to the console
  print(paste("saving file",myPlotFilename,"..."))
  # 7) save the image
  ggsave(myPlotFilename,myPlot)
  # after ending up here, the code will be executed again from the start for the next image
}

8. 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).

# 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

# first, get all the files 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 only
files=files[grep("h5$",files)]
# let's process the first file
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 profile,
# similar to the profile files stored in the ENRAM repository
# 
# let's make an output file (by replacing the substring 'pvol' into 'vp')
file.out=gsub("pvol","vp",file.in)
# print the output file name to which we will write:
file.out
# finally, run vol2bird; this generates an output file as specified in file.out
vol2bird(file.in,file.out)
# your work directory should now contain a new file 'seang_vp_20150927T055000Z.h5'
# check that we can read this file:
vp=readvp(file.out)
  1. write your own for-loop to process a batch of polar volume files. Hint: use a similar approach as the for-loop of the last example in section 7, in which you processed images.