Exercise 1

1.1. Preparation

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 provided below 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.13). If you have an older version, update as follows:
library(devtools)
install_github("adokter/bioRad")

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/"
# for example, in Adriaan's case:
HOME="~/Dropbox/teaching/buler2018"
# 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

To work with US NEXRAD data, bioRad currently require 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). If you managed to successfully install Docker, test whether it works

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

If you did not manage to install Docker, you will not be able to load NEXRAD data into R at this time, but you will be able to continue the exercises with European radar data that can be loaded without Docker.

1.2. Downloading and loading polar volume data

US NEXRAD data can be accessed in the Amazon cloud

Let us download this polar volume file for the KPAH radar (Paducah, KY), and store it in your working directory.

If you don’t have a working docker container, you can download this polar volume file, and store it in your working directory. It’s the same file but in a different format (only difference is that it contains only data up to 35 km range)

1.3. The structure of polar volumes

#  let's read in the downloaded volume:
file.in="KPAH20171025_040338_V06"
# or in case you downloaded the file that doesn't require docker, the filename is different:
file.in="volume.h5"
# check that the file is stored in the right location:
file.exists(file.in)
# load the polar volume:
pvol=read.pvol(file.in)
## print some information about the polar volume
pvol
# print information about the polar scans contained in this polar volume:
pvol$scans

1.4. Inspecting radar scans (sweeps)

# let's extract the third scan, which was collected at 0.48 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,cellsize=1000,range.max=100000)
# 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="VRADH",zlim=c(-50,50))
# see plot.ppi for all the plot options for a ppi object:
?plot.ppi

1.5. Overlaying radar scans on maps

# It is often informative to plot radar data on a base layer, such as google earth maps.
# first download the background image:
satelliteImage=basemap(my_ppi,maptype="satellite")
# then overlay the PPI on the satellite image:
map(my_ppi,param="DBZH",map=satelliteImage, zlim=c(-20,15))
# Note that in R, spatial data is often contained in class objects of packag 'sp'
# bioRad also uses these objects in the background, and they can be extracted if you want to.
# The spatial data is stored in the data slot, as in:
my_spatialgrid=my_ppi$data
# you can use the sp package to save the spatial data to all kinds GIS formats, for example ArcGis:
library(sp)
write.asciigrid(my_spatialgrid,"PPI_in_arcgis_format.asc")

Exercise 2

2.1. Processing polar volume data into vertical profiles

We will generate vertical profiles with the automated algorithm vol2bird (https://github.com/adokter/vol2bird), which is included in the bioRad package. 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 process the same file as in section 1.3:
file.in
# check whether the file is still there:
file.exists(file.in)
# run vol2bird
vp=vol2bird(file.in,range.max=35000,sd_vvp=1,dealias=T,dualpol=T)
# vp is now a 'vp', a vertical profile
vp
# alternatively, you may also store the profile on disk as a hdf5 file, which is what we will do next:
# let's first define the name of the output file (we paste the extention ".h5" to the name)
file.out=paste(file.in,".h5",sep="")
# 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 with the name you specified in file.out
# check that we can read this file, and retrieve the vertical profile from it:
vp=readvp(file.out)

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

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

# print the currently assumed radar cross section (RCS) per bird:
rcs(vp)
# plot the vertical profile, in terms of bird density
plot(vp, quantity="dens")

We can change the assumed radar cross section as follows:

# let's change the RCS to 110 cm^2
rcs(vp) = 110
plot(vp, quantity="dens")

Exercise 3

3.1. Loading processed vertical profiles

To save some time, we will continue with some pre-processed vertical profiles for the KPAH radar.

Navigate to https://www.dropbox.com/s/6em8s4yjqutkfle/profiles.zip?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 'profiles' with processed data
# load all the filenames in the new 'profiles' subfolder in your working directory
vp_paths=dir("./profiles", recursive=TRUE,full.names=TRUE)
vp_paths
# read these vertical profiles (hdf5 files) into R (may take a minute to load)
vplist=readvp.list(vp_paths)
# print some information on the vplist object. It should contain 71 profiles 
vplist
# save the object, which allows you to load the data more quickly next time
save(vplist,file="vplist.RData")
# you can restore the vplist object at any time as follows:
load("vplist.RData")

3.2. 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:
# In case your vplist contains profiles of different radars, this function will
# group the profiles and make a separate time-series for each radar.
# (but in our case we have profiles for only one radar)
ts=vpts(vplist)
# print summary information
ts
# plot the time series in terms of reflectivity factors, from 0-2000 metre altitude:
plot(ts, quantity="dbz",ylim=c(0,2000))
# plot the time series in terms of bird density:
plot(ts, quantity="dens",ylim=c(0,2000))
# change the radar cross-section works the same as for single vertical profiles, for example:
rcs(ts)=20
# check the help file for more plotting options
# Because profile timeseries are of class 'vpts', it's associated plotting function
# is plot.vpts:
?plot.vpts

To interpret the wind barbs in the profile time series figure: 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]

Precipitation is known to have a major influence on the timing and intensity of migration, therefore it can be useful to inspect profiles for presence of precipitation. 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.

# plot the time series
plot(ts, quantity="DBZH")
# You can also extract the data from bioRad's class objects into a simple R format:
#
# Extract some data from the time series, e.g. the bird density
fetch(ts, quantity="dens")
# convert all the data in the time series to a standard data.frame:
my_dataframe=as.data.frame(ts)
my_dataframe

3.3 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)
# plot the integrated data to screen:
vintegrated.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")

Now let’s calculate how many birds passed over the radar during the full time series:

# Let's continue with the ts object created in the previous example.
# we can integrate all the traffic rates [unit: birds/km/h] over time, to obtain the number of birds that
# have passed over the radar [unit: birds/km]. Like the traffic rates, these numbers are calculated for 
# a 1 km transect perpendicular to the migratory direction:
mt(ts)