5 Range bias correction

In a previous section we have explored vertical integration of vertical profiles (vp’s). In this paragraph we will generalize vertical integration to entire radar images (ppi’s). This requires us to account for the changing beam shape of the radar with range.

First, let’s examine the beam shape of the lowest elevation scan of the radar, which is typically around 0.5 degrees.

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

We will start with downloading a polar volume and processing it into a profile:

5.1 Processing a polar volume into a profile

# download a polar volume for the KBRO radar in Brownsville, TX
my_pvol <- get_pvol("KBRO",as.POSIXct("2017-05-14 05:50:00"))
# calculate the profile, correlation-coefficient is used automatically to remove precipitation:
# we calculate 60 layers of 50 meter width, so up to 60*50=3000 m.
my_vp <- calculate_vp(my_pvol, n_layer=60, h_layer=50)
# also clean the pvol object for precipitation using mistnet:
my_pvol <- apply_mistnet(my_pvol)
# set all the DBZH pixels identified by mistnet as precipitation to NA:
my_pvol_clean <- calculate_param(my_pvol, DBZH = ifelse(CELL == 2, NA, DBZH))

Exercise 10: Plot the bird density in the vertical profile you just estimated and compare it with the plots above of the beam height and width of the lowest radar beam. At which approximate range do you expect the radar will no longer be able to resolve the altitude distribution of the migratory birds. And at which range will the radar start to overshoot the migration layer entirely?

To correct for the decreasing ability of the radar to resolve the altitude distributions of birds with range, we will make one important (and likely imperfect!) assumption:

We assume that all birds within the image are distributed vertically according to the same relative vertical profile.

This assumptions simplifies the problem, and allows us to estimate the spatial distribution of the birds, as we will explore in the next paragraph:

5.2 Range bias correction and vertical integration on a map

# We will use this helper function to create a base grid centered around the radar
# distance is the maximum distance from the radar for the grid, res is
# the pixels size, both in meters
define_base_grid <- function(pvol, distance=100000, res=1000, shape_only=FALSE){
  # define the base grid area
  base_grid_area <-
    # define a point at the location of the radar in WGS84 coordinates
    sf::st_sfc(sf::st_point(c(pvol$attributes$where$lon, pvol$attributes$where$lat)), crs = 4326) |>
    # convert to Albers Equal-Area Conic projection, which is good for preserving the relative sizes
    # of areas in the contiguous USA.
    sf::st_transform(crs = 5070) |>
    # create a 100 km buffer radius around the radar location
    sf::st_buffer(dist = distance) |>
    # convert to sf dataframe
    sf::st_as_sf()
  # if shape_only is TRUE, only return the circular shape:
  if(shape_only) return(base_grid_area)
  # otherwise, return a full terra raster:
  # create a raster for the base grid area
  terra::rast(
    ext = terra::ext(base_grid_area),  # extent from the circular area
    resolution = res,          # 1000m resolution
    crs = terra::crs(base_grid_area)   # same CRS as the circular area
  )
}

# define the base grid
base_grid <- define_base_grid(my_pvol_clean)
# Calculated the vertically integrated PPI
my_ppi_integrated <- integrate_to_ppi(pvol=my_pvol_clean,vp=my_vp, raster=base_grid)

Next, we will compare the PPI for the 1.0 degree sweep and the vertically integrated PPI, to highlight the effect of the range correction and vertical integration. We will plot the log-transform of VID for a more direct comparison with the reflectivity factor DBZH, which is also logarithmic in bird density.

# Project the 1 degree sweep on the same basegrid:
my_pvol_clean |>
  get_scan(elev = 1.0) |>
  project_as_ppi(raster=base_grid) ->
  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))

Exercise 11: Explain why we see concentric rings in the 1.0 degree sweep, and not in the vertically integrated PPI.

5.3 Stopover mapping

Birds taking off tend to still be closely associated with ground level. In areas with topography, the common assumption is that birds will still take off simultaneously, resulting in an altitude distribution during the take-off phase that is best estimated relative to ground level.

To calculate a profile referenced to ground level, we first define a base grid and download digital elevation information:

# Download a polar volume for KCCX radar
# This in the Appalachian mountains of Central Pennsylvania
pvol <- getRad::get_pvol("KCCX", as.POSIXct("2013-10-08 23:24:30", tz="UTC"))
# This polar volume is collected 45 mins after local sunset.
# Birds are still in the process of take-off during this time.
pvol$datetime - sunset(pvol$datetime, pvol$geo$lon, pvol$geo$lat)
# define the base grid area:
base_grid_area <- define_base_grid(pvol, shape_only=TRUE)
# download digital elevation grid (z defines resolution), and name the data "HGHT"
base_grid_ground <- elevatr::get_elev_raster(base_grid_area, z=7, clip = "locations")
names(base_grid_ground) = "HGHT"
# Plot the topography as a reference:
plot(rast(base_grid_ground))
# add ground height data as a scan parameter to the polar volume:
pvol_ground <- add_param(pvol, raster=base_grid_ground, param="HGHT")
# verify that ground height information is now added to each elevation scan:
plot(get_scan(pvol_ground,.5), param="HGHT")

Next, we can calculate a vertical profile relative to ground level:

# calculate a profile relative to ground level
vp_ground <- calculate_vp(pvol_ground, height_reference = "ground")
# calculate a range-corrected image, assuming birds maintain a constant
# altitude profile relative to ground level:
my_ppi_ground <- integrate_to_ppi(pvol=pvol_ground,vp=vp_ground, raster=base_grid_ground, height_reference="ground")
# Mask out areas with very strong range correction that are unreliable.
# In this case we require eta_sum_to_VIR<=8:
my_ppi_ground <- calculate_param(my_ppi_ground, VIR=ifelse(eta_sum_to_VIR>8,NA,VIR))
# Plot the result, using viridis color scale:
plot(my_ppi_ground, zlim=c(0,2000)) + viridis::scale_fill_viridis(name = "VIR")

Exercise 12: Create the same stopover map, but instead of assuming birds maintain a constant altitude distribution relative to ground level, assume they maintain a constant altitude distribution relative to sea level. Explain observed differences in the two maps.