Nothing
## ----setup, echo=FALSE, message=FALSE-----------------------------------------
SHOW_ANSWERS <- FALSE
if (Sys.info()["sysname"] == "Linux") prefix <- "/home/adriaan/" else prefix <- "/Users/amd427/"
if (SHOW_ANSWERS) knitr::opts_knit$set(root.dir = normalizePath(paste(prefix, "git/training/colorado2022/", sep = "")))
# knitr::opts_chunk$set(eval=FALSE)
Sys.setenv(TZ = "UTC")
library(bioRad)
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # make sure you start with a fresh R session
# # load the bioRad package
# library(bioRad)
# # check the package version
# packageVersion("bioRad")
## ----eval=FALSE---------------------------------------------------------------
# # bring up the package general help page:
# ?bioRad
## ----eval=FALSE---------------------------------------------------------------
# # make a new local directory on your machine 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)
# # next, we will create two directories where we will be storing data:
# dir.create("./data_vpts") # here we will store vertical profile time series (vpts data)
# dir.create("./data_pvol") # here we will store polar volumes (pvol data)
# # Finally, we set the local time zone to UTC, so all plotted time axes will be in UTC
# Sys.setenv(TZ = "UTC")
## ----eval=SHOW_ANSWERS, results='hide'----------------------------------------
# # Let's first download the NEXRAD polar volume files for the KHGX radar (Houston)
# # for a 15 minute period in 2017:
# download_pvolfiles(date_min=as.POSIXct("2017-05-04 01:25:00"), date_max=as.POSIXct("2017-05-04 01:40:00"), radar="KHGX", directory="./data_pvol")
# # store the filenames in my_pvolfiles
# my_pvolfiles <- list.files("./data_pvol", recursive = TRUE, full.names = TRUE, pattern="KHGX")
# # print to console our files:
# my_pvolfiles
# # let's load the first of our downloaded files:
# my_pvol <- read_pvolfile(my_pvolfiles[1])
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS-------------------------------------
# # The default summary information of a `pvol` object contains information
# # on the scans (sweeps) and their moments:
# my_pvol
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS-------------------------------------
# # We can also extract the elevation angles from the polar volume as follows:
# get_elevation_angles(my_pvol)
## ----eval=SHOW_ANSWERS, warning=FALSE-----------------------------------------
# # let's extract the scan collected at 1.5 degree elevation from our polar volume:
# my_scan <- get_scan(my_pvol, 0.5)
# # print some information about this scan:
# my_scan
# # let's plot the reflectivity factor parameter of the scan in a range - azimuth coordinate system:
# plot(my_scan, param = "DBZH")
## ----eval=SHOW_ANSWERS, warning=FALSE-----------------------------------------
# # 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 <- project_as_ppi(my_scan)
# # print some information about this ppi:
# my_ppi
# # you can see we projected it on a 500 meter grid
# # (check the manual of the project_as_ppi function to see how you can
# # change the grid size (argument grid_size) and the maximum distance
# # from the radar up to where to plot data (argument range_max))
# #
# # Now we are ready to plot the ppi, for example let's plot reflectivity factor DBZH:
# plot(my_ppi, param = "DBZH")
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS, warning=FALSE----------------------
# # Plot the correlation coefficient (RHOHV):
# plot(my_ppi, param = "RHOHV")
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS, warning=FALSE----------------------
# # Plot the radial velocity (VRADH):
# plot(my_ppi, param = "VRADH")
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS, warning=FALSE----------------------
# # Answer:
# # * Precipitation areas are characterized by high correlation coefficients (typically > 0.95),
# # i.e. the top left corner of the image is precipitation
# # * The blue radial velocity of the precipitation indicates it is moving towards the radar.
# # * The radial velocity field of the biology shows areas south-west of the radar moving
# # towards the radar (blue), and areas north-east of the radar moving away from it (red).
# # The biology is therefore moving towards the north-east
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS, warning=FALSE----------------------
# # Answer:
# # * Precipitation moves with the wind field. Since the biology moves into a different direction
# # than the wind and a very different speed, we can be sure these are birds.
# # The biological scatterers have a high self-propelled speed, which is
# # typical for birds, not for insects.
# # * Note: This is an S-band radar. In C-band radars you will typically see that the texture
# # (spatial variability) of the radial velocity is considerably smoother in areas with
# # precipitation than in areas with biology.
## ----eval=SHOW_ANSWERS, warning=FALSE-----------------------------------------
# # It is often informative to plot radar data on a base layer.
# # First choose a base layer from the list of rosm::osm.types()
# basemap = "osm"
# # then overlay the PPI on the basemap, restricting the color scale from -20 to 40 dBZ:
# map(my_ppi, map = basemap, param = "DBZH", zlim = c(-20, 40))
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # Screen out the reflectivity areas with RHOHV < 0.95
# my_ppi_clean <- calculate_param(my_ppi, DBZH = ifelse(RHOHV > 0.95, NA, DBZH))
# # plot the original and cleaned up reflectivity:
# map(my_ppi, map = basemap, param = "DBZH", zlim = c(-20, 40))
# map(my_ppi_clean, map = basemap, param = "DBZH", zlim = c(-20, 40))
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # apply the MistNet model to the polar volume file and load it as a polar volume (pvol):
# my_pvol <- apply_mistnet(my_pvolfiles[1])
# # mistnet will add additional parameters to the
# # elevation scans at 0.5, 1.5, 2.5, 3.5 and 4.5 degrees
# # let's extract the scan closest to 0.5 degrees:
# my_scan <- get_scan(my_pvol, 0.5)
# # plot some summary info about the scan to the console:
# my_scan
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # let's add depolarization ratio (DR) as a parameter (following Kilambi 2018):
# my_ppi <- calculate_param(my_ppi, DR = 10 * log10((1+ ZDR - 2 * (ZDR^0.5) * RHOHV) /
# (1 + ZDR+ 2 * (ZDR^0.5) * RHOHV)))
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # plot the depolarization ratio, using a viridis color palette:
# map(my_ppi, map = basemap, param = "DR", zlim=c(-25,-5), palette = viridis::viridis(100))
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS-------------------------------------
# # Now let us screen out the reflectivity areas with DR < -12 dB:
# my_ppi_clean <- calculate_param(my_ppi, DBZH = ifelse(DR < -12, NA, DBZH))
# # plot the original and cleaned up reflectivity:
# map(my_ppi, map = basemap, param = "DBZH", zlim = c(-20, 40))
# map(my_ppi_clean, map = basemap, param = "DBZH", zlim = c(-20, 40))
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # as before, project the scan as a ppi:
# my_ppi <- project_as_ppi(my_scan)
# # plot the probability for the WEATHER class
# plot(my_ppi, param = 'WEATHER')
# # plot the final segmentation result:
# # plot the probability for the WEATHER class
# plot(my_ppi, param = 'CELL')
# # let's remove the identified precipitation area (and additional fringe) from the ppi, and plot it:
# my_ppi_clean <- calculate_param(my_ppi, DBZH = ifelse(CELL >= 1, NA, DBZH))
# map(my_ppi_clean, map=basemap, param = 'DBZH')
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS-------------------------------------
# # To remove the additional fringe around the rain segmentation by MistNet
# # we want to screen out CELL values equal to 2 only.
# my_ppi_clean <- calculate_param(my_ppi, VRADH = ifelse(CELL >= 2, NA, VRADH))
# map(my_ppi_clean, map=basemap, param = 'VRADH')
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # Usually we would load processed vertical profiles (vp files) by:
# # my_vplist <- read_vpfiles("./your/directory/with/processed/profiles/goes/here")
# # my_vplist contains after running the command a list of vertical profile (vp) objects
# # To save time, we load these data directly from file
# my_vplist <- readRDS("data_vpts/KBRO20170514.rds")
# # print the length of the vplist object. It should contain 95 profiles
# length(my_vplist)
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # let's extract a profile from the list, in this example the 41st profile:
# my_vp <- my_vplist[[41]]
# # print some info for this profile to the console
# my_vp
# # test whether this profile was collected at day time:
# check_night(my_vp)
# # plot the vertical profile, in terms of reflectivity factor
# plot(my_vp, quantity = "dbz")
# # plot the vertical profile, in terms of (linear) reflectivity
# plot(my_vp, quantity = "eta")
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # let's plot the vertical profile, in terms of bird density
# plot(my_vp, quantity = "dens")
# # print the currently assumed radar cross section (RCS) per bird:
# rcs(my_vp)
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS, warning=FALSE----------------------
# # Answer:
# # All bird densities will be a factor 10 times lower.
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # let's change the RCS to 110 cm^2
# rcs(my_vp) <- 110
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS, warning=FALSE----------------------
# # After changing the RCS above, we simply plot the vertical profile again
# plot(my_vp)
# # Indeed the densities are scaled down by a factor 10
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # convert the list of vertical profiles into a time series:
# my_vpts <- bind_into_vpts(my_vplist)
# # time series objects can be subsetted, just as you may be used to with vectors
# # here we subset the first 50 timesteps:
# my_vpts[1:50]
# # here we extract a single timestep, which gives you back a vertical profile class object:
# my_vpts[100]
# # to plot the full time series:
# plot(my_vpts)
# # 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
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # filter our vpts for night time
# my_vpts_night <- filter_vpts(my_vpts, night=TRUE)
# # plot this smaller time series:
# plot(my_vpts_night)
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS, warning=FALSE----------------------
# # Answer:
# #
# # At 1500 meter 6 UTC the wind barbs have 2 full flags and one half flag.
# # Therefore the ground speed is approximately 2x5 + 2.5 = 12.5 m/s.
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS, warning=FALSE----------------------
# # First extract the profile at 6 UTC:
# vp_6UTC <- filter_vpts(my_vpts_night, nearest = "2017-05-14 06:00")
# plot(vp_6UTC, quantity = "ff")
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS, warning=FALSE----------------------
# # Plot the ground speed (ff):
# plot(vp_6UTC, quantity = "ff")
# # Speeds at 1500 metre is approximately 12 m/s, confirming our earlier reading above of 12.5 m/s.
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # Let's continue with the vpts object created in the previous example.
# # The vertically integrated quantities are calculated as follows:
# my_vpi <- integrate_profile(my_vpts)
# # The my_vpi object you created is a vpi class object, which is an acronym for "vertical profile integrated". It has its own plot method, which by default plots migration traffic rate (MTR):
# plot(my_vpi)
# # you can also plot vertically integrated densities (VID):
# plot(my_vpi, 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(my_vpi, night_shade = FALSE)
# # plot the cumulative number of birds passing the radar, i.e. migration traffic (mt):
# plot(my_vpi, quantity = "mt")
# # execute `?plot.vpi` to open the help page listing all the options.
# ?plot.vpi
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS, warning=FALSE----------------------
# # Answer:
# #
# # VID = (200 birds / km^3) * (1 km) + (100 birds / km^3) * (0.5 km)
# # = 250 birds / km^2
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS, warning=FALSE----------------------
# # Answer:
# #
# # MTR = (200 birds / km^3) * (50 km / hour) * (1 km) + (100 birds / km^3) * (100 km / hour) * (0.5 km)
# # = 10000 birds / km / hour + 5000 birds / km / hour
# # = 15000 birds / km / hour
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS, warning=FALSE----------------------
# # Answer:
# #
# # MT = MTR * (3 hour) = (15000 birds / km / hour) * (3 hour)
# # = 45000 birds / km
# # for a 10 kilometer transect:
# # Number_of_birds = 10 km * 45000 birds / km = 450000 birds
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # instead of vertically integrated density (VID), you can use vertically integrated reflectivity (VIR):
# plot(my_vpi, quantity = "vir")
# # instead of migration traffic rate (MTR), you can use the reflectivity traffic rate (RTR):
# plot(my_vpi, quantity = "rtr")
# # instead of migration traffic (MT), you can use the reflectivity traffic (RT):
# plot(my_vpi, quantity = "rt")
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # load a time series for the KBGM radar in Binghamton, NY
# my_vpts <- readRDS("data_vpts/KBGM20170527-20170602.rds")
# # print the loaded vpts time series for this radar:
# my_vpts
# # plot the bird density over time:
# plot(my_vpts, quantity = "dens")
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS, warning=FALSE----------------------
# # also show meteorological signals:
# plot(my_vpts, quantity = "DBZH")
#
# # Periods with high reflectivities extending to high altitudes indicate precipitation,
# # i.e. second half of the second night, and on and off during the fourth night.
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # define ranges from 0 to 2500000 meter (250 km), in steps of 100 m:
# range <- seq(0, 250000, 100)
#
# # plot the beam height of the 0.5 degree elevation beam:
# plot(range, beam_height(range, 0.5), ylab = "beam height [m]", xlab = "range [m]", type='l')
#
# # let's add the lower and upper altitude of the beam, as determined by the beam width:
# points(range, beam_height(range, 0.5)-beam_width(range)/2, type='l',lty=3)
# points(range, beam_height(range, 0.5)+beam_width(range)/2, type='l',lty=3)
## ----eval=SHOW_ANSWERS, warning=FALSE, results='hide'-------------------------
# # download a polar volume for the KBRO radar in Brownsville, TX
# download_pvolfiles(date_min=as.POSIXct("2017-05-14 05:50:00"), date_max=as.POSIXct("2017-05-14 06:00:00"), radar="KBRO", directory="./data_pvol")
# # Load all the polar volume filenames downloaded so far for the KBRO radar:
# my_pvolfiles <- list.files("./data_pvol", recursive = TRUE, full.names = TRUE, pattern="KBRO")
# # we will process the first one into a vp:
# my_pvolfile <- my_pvolfiles[1]
# # calculate the profile, using MistNet to remove precipitation:
# # we calculate 60 layers of 50 meter width, so up to 30*100=3000 m.
# my_vp <- calculate_vp(my_pvolfile, n_layer=60, h_layer=50, sd_vvp_threshold = 1)
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS, warning=FALSE----------------------
# # plot the profile
# plot(my_vp)
#
# # Answer:
# # * The vertical profile shows altitude bands of about 300 meter width
# # We therefore expect the radar to start having difficulty to resolve this profile
# # when the beam shape becomes broader, so less than 50 km
# # (this is why we typically use 35 km as the maximum range to use in vertical profile estimation.)
# # * the birds fly up to about 2km in this profile.
# # At about 200000 m (200 km) the lower end of the beam no longer overlaps with the migration layer
# # Therefore, if we birds are flying according to the same altitude distribution everywhere,
# # we expect the radar to become fully blind for birds beyond 200 km
## ----eval=SHOW_ANSWERS,warning=FALSE------------------------------------------
# # We will use the piping operator %>% of magrittr package to
# # execute multiple operations in one statement:
# library(magrittr)
# # first, load the polar volume:
# my_pvolfile %>% read_pvolfile() -> my_pvol
# # Next, let's calculate a PPI for the 1.5 degree elevation scan
# # Finally, we calculated the vertically integrated PPI
# my_ppi_integrated <- integrate_to_ppi(pvol=my_pvol,vp=my_vp, res=1000)
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS,warning=FALSE-----------------------
# #
# my_pvol %>%
# get_scan(elev = 1.0) %>%
# project_as_ppi(raster=raster::raster(my_ppi_integrated$data)) ->
# my_ppi
# # let's take the logarithm of the vertically integrated density, so
# # we can compare it more directly to DBZH, which is also a logarithmic quantity:
# my_ppi_integrated <- calculate_param(my_ppi_integrated, logVID=log(VID))
# # plot both ppi's:
# plot(my_ppi)
# plot(my_ppi_integrated, param="logVID", zlim=c(-10,10))
## ----echo=SHOW_ANSWERS, eval=SHOW_ANSWERS-------------------------------------
# # The altitude distribution in this case shows two migration layers.
# # In the 1.5 PPI scan these show up as two concentric rings of density.
# # Since in a conventional PPI a larger distance from the radar also
# # implies a higher altitude, we see the rings showing up at the
# # ranges where the beam intersects each of the migration layers.
# #
# # The vertically integrated PPI is more straightforward to interpret
# # Here we correct for aforementioned beam effects and integrate the
# # density over altitude, giving a more realistic reconstruction of the
# # spatial distribution of migrants.
## ----eval=SHOW_ANSWERS,results='hide',warning=FALSE---------------------------
# # First we download more data, for a total of one additional hour for the same radar:
# download_pvolfiles(date_min=as.POSIXct("2017-05-04 01:40:00"), date_max=as.POSIXct("2017-05-04 02:40:00"), radar="KHGX", directory="./data_pvol")
# # We will process all the polar volume files downloaded so far:
# my_pvolfiles <- list.files("./data_pvol", recursive = TRUE, full.names = TRUE, pattern="KHGX")
# my_pvolfiles
# # create output directory for processed profiles
# outputdir <- "./data_vp"
# dir.create(outputdir)
## ----eval=FALSE, echo=FALSE---------------------------------------------------
# # let's loop over the files and generate profiles
# start=Sys.time()
# for (file_in in my_pvolfiles) {
# # generate the output filename for the input file
# file_out <- paste(outputdir, "/", basename(file_in), "_vp.h5", sep = "")
# # generate the profile, and write it as a hdf5 file (for later reference)
# # we enclose the calculate_vp statement in tryCatch to keep going after errors with specific files
# vp <- tryCatch(calculate_vp(file_in, file_out, mistnet = TRUE, h_layer=50, n_layer=60), error = function(e) NULL)
# }
# end=Sys.time()
# # calculate the processing time that has passed:
# end-start
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # we assume outputdir contains the path to the directory with processed profiles
# my_vpfiles <- list.files(outputdir, full.names = TRUE, pattern="KHGX")
# # print them
# my_vpfiles
# # read them
# my_vplist <- read_vpfiles(my_vpfiles)
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# # make a time series of profiles:
# my_vpts <- bind_into_vpts(my_vplist)
# # plot them between 0 - 3 km altitude:
# plot(my_vpts, ylim = c(0, 3000))
# # because of the rain moving in, our ability to estimate bird profiles slowly drops:
# # let's visualize rain by plotting all aerial reflectivity:
# plot(my_vpts, ylim = c(0, 3000), quantity="DBZH")
#
## ----eval=SHOW_ANSWERS--------------------------------------------------------
# process_file <- function(file_in){
# # construct output filename from input filename
# file_out <- paste(outputdir, "/", basename(file_in), "_vp.h5", sep = "")
# # run calculate_vp()
# vp <- tryCatch(calculate_vp(file_in, file_out, mistnet = FALSE, h_layer=50, n_layer=60), error = function(e) NULL)
# if(is.vp(vp)){
# # return TRUE if we calculated a profile
# return(TRUE)
# } else{
# # return FALSE if we did not
# return(FALSE)
# }
# }
#
# # To process a file, we can now run
# process_file(my_pvolfiles[1])
## ----eval=FALSE---------------------------------------------------------------
# # load the parallel library
# library(parallel)
# # detect how many cores we can use. We will keep 2 cores for other tasks and use the rest for processing.
# number_of_cores = detectCores() - 2
# # Available cores:
# number_of_cores
# # let's loop over the files and generate profiles
# start=Sys.time()
# mclapply(my_pvolfiles, process_file, mc.cores = number_of_cores)
# end=Sys.time()
# # calculate the processing time that has passed:
# end-start
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.