inst/doc/rad_aero_22.R

In bioRad: Biological Analysis and Visualization of Weather Radar Data

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