bird insect partitioning and post-hoc precipitation filter

NAMESPACE

@@ -29,6 +29,9 @@ S3method(check_night,pvol)
 S3method(check_night,vp)
 S3method(check_night,vpi)
 S3method(check_night,vpts)
+S3method(clean_mixture,default)
+S3method(clean_mixture,vp)
+S3method(clean_mixture,vpts)
 S3method(convert_legacy,vp)
 S3method(convert_legacy,vpts)
 S3method(dim,ppi)
@@ -99,6 +102,7 @@ export(calculate_param)
 export(calculate_vp)
 export(check_docker)
 export(check_night)
+export(clean_mixture)
 export(composite_ppi)
 export(convert_legacy)
 export(dbz_to_eta)
@@ -107,6 +111,7 @@ export(download_pvolfiles)
 export(download_vpfiles)
 export(doy)
 export(eta_to_dbz)
+export(filter_precip)
 export(filter_vpts)
 export(get_elevation_angles)
 export(get_iris_raw_task)

R/as.vp.R

@@ -0,0 +1,19 @@
+#' Convert a dataframe into a vp object
+#'
+#' @param data a dataframe created from a VPTS CSV file
+#' @returns a bioRad vp object
+#' @examples
+#' # load vp data as a data.frame:
+#' df <- as.data.frame(example_vp)
+#' # convert the data.frame to a vp object:
+#' as.vp(df)
+#' @export
+as.vp <- function(data) {
+  assertthat::assert_that(inherits(data,"data.frame"))
+
+  vpts <- as.vpts(data)
+
+  assertthat::assert_that(length(vpts$datetime) == 1, msg="multiple timestamps found, data is not a single vertical profile")
+
+  vpts_to_vp(vpts)
+}

R/clean_mixture.R

@@ -0,0 +1,175 @@
+#' Partition mixtures of animals using assumptions on airspeeds.
+#'
+#' Partition mixtures of animals using assumptions on airspeeds.
+#' @param x a `vp` or `vpts` object.
+#' @param ... `eta`, `u`, `v`, `U`, `V` arguments, taken from object for `vp` or `vpts` class.
+#' @param eta a mixture animal density or linear reflectivity eta.
+#' @param u the mixture's ground speed u component (west to east)
+#' @param v the mixture's ground speed v component (south to north)
+#' @param U the west to east wind component
+#' @param V the south to north wind component
+#' @param fast the fast component's airspeed
+#' @param slow the slow component's airspeed
+#' @param drop_slow_component when TRUE (default) output density, ground speed and
+#' heading for fast component, when FALSE for slow component.
+#' @param drop_missing Values `eta` without an associated ground speed
+#' and wind speed are set to NA when `TRUE`, or returned unaltered when `FALSE` (default).
+#' @param keep_mixture When `TRUE` store original mixture reflectivity and speeds as
+#' renamed quantities with `mixture_` prefix
+#' @return a data.frame with corrected density or reflectivity x, ground speed (u, v)
+#' and heading in clockwise degrees from north.
+#' @export
+#' @name clean_mixture
+#' @examples
+#' # drop the slow component (typically insects)
+#' clean_mixture(100,-13,13,-7,6, fast=8, slow=1)
+#' # drop the fast component (typically birds)
+#' clean_mixture(100,-13,13,-7,6, fast=8, slow=1, drop_slow_component=FALSE)
+#' # keep the original mixture reflectivity and speed components
+#' clean_mixture(100,-13,13,-7,6, fast=8, slow=1, keep_mixture=TRUE)
+#' # keep reflectivity unaltered when one of the speed components is not a number:
+#' clean_mixture(100,-13,13,NaN,6, fast=8, slow=1)["eta"]
+#' # set reflectivity to NaN when one of the speed components is not a number:
+#' clean_mixture(100,-13,13,NaN,6, fast=8, slow=1, drop_missing=TRUE)["eta"]
+NULL
+
+#' @rdname clean_mixture
+#'
+#' @export
+clean_mixture <- function(x, ..., slow = 1, fast = 8, drop_slow_component = TRUE, drop_missing = FALSE, keep_mixture = FALSE){
+  UseMethod("clean_mixture", x)
+}
+
+#' @rdname clean_mixture
+#'
+#' @export
+clean_mixture.default <- function(eta, u, v, U, V, slow = 1, fast = 8, drop_slow_component = TRUE, drop_missing = FALSE, keep_mixture = FALSE){
+  # verify input
+  assertthat::assert_that(all(eta >= 0, na.rm = TRUE))
+  assertthat::assert_that(is.numeric(u))
+  assertthat::assert_that(is.numeric(v))
+  assertthat::assert_that(is.numeric(U))
+  assertthat::assert_that(is.numeric(V))
+  assertthat::assert_that(slow >= 0)
+  assertthat::assert_that(fast > 0)
+  assertthat::assert_that(fast > slow)
+  assertthat::assert_that(assertthat::is.flag(drop_slow_component))
+  assertthat::assert_that(assertthat::is.flag(drop_missing))
+  assertthat::assert_that(assertthat::is.flag(keep_mixture))
+
+  if(!all(is.finite(U))) warning("U component of wind contains non finite values")
+  if(!all(is.finite(V))) warning("V component of wind contains non finite values")
+
+  # define helper quantities:
+  wind_speed <- sqrt(U^2 + V^2)
+  wind_direction <- atan2(V,U)
+  mixture_airspeed <- sqrt((u-U)^2 + (v-V)^2)
+  mixture_heading <- (pi / 2 - atan2(v-V, u-U)) * 180 / pi
+  mixture_heading[which(mixture_heading<0)]=mixture_heading[which(mixture_heading<0)]+360
+
+  # solve mixture equations for slow proportion f
+  p1 <- slow^2 - fast^2
+  p2 <- 2*fast^2 - 2*(slow/wind_speed)*(u*U+v*V - wind_speed^2)
+  p3 <- (u-U)^2 + (v-V)^2 - fast^2
+  f <- (-p2+sqrt(p2^2-4*p1*p3))/(2*p1)
+
+  # catch limiting cases for which f is not defined
+  idx_f_zero <- mixture_airspeed > fast
+  if(length(which(idx_f_zero))>0){
+    warning("Assigning all reflectivity to fast component for mixture airspeeds exceeding airspeed of fast component")
+    f[idx_f_zero]=0
+  }
+  idx_f_one <- mixture_airspeed < slow
+  if(length(which(idx_f_one))>0){
+    warning("Assigning all reflectivity to slow component for mixture airspeeds below airspeed of slow component")
+    f[idx_f_one]=1
+  }
+
+  # initialize corrected eta and speed matrix or vector:
+  eta_corr=eta
+  air_u = u-U
+  air_v = v-V
+
+  # initialize and define indices to correct densities and speeds
+  idx = !is.na(f) & !idx_f_zero & !idx_f_one
+
+  if(drop_slow_component){
+    # fast component airspeed, typically birds:
+    eta_corr[idx]=((1-f)*eta)[idx]
+    air_u[idx]=(((u-U)-(slow/wind_speed)*U*f)/(1-f))[idx]
+    air_v[idx]=(((v-V)-(slow/wind_speed)*V*f)/(1-f))[idx]
+  } else{
+    # slow component airspeed, typically insects:
+    eta_corr[idx]=(f*eta)[idx]
+    air_u[idx]=slow*cos(wind_direction)[idx]
+    air_v[idx]=slow*sin(wind_direction)[idx]
+  }
+
+  if(drop_missing){
+    eta_corr[is.na(f)]=NaN
+    air_u[is.na(f)]=NaN
+    air_v[is.na(f)]=NaN
+  }
+
+  # calculate speed and direction
+  ff <- sqrt((U+air_u)^2 + (V+air_v)^2)
+  dd <- (pi / 2 - atan2(V+air_v, U+air_u)) * 180 / pi
+  dd[which(dd<0)]=dd[which(dd<0)]+360
+
+  output <- list(eta=eta_corr,u=U+air_u,v=V+air_v, airspeed = sqrt(air_u^2+air_v^2), heading=(pi/2-atan2(air_v,air_u))*180/pi,
+                 airspeed_u = air_u, airspeed_v = air_v, f=f)
+
+  if(keep_mixture){
+    output = append(output, list(mixture_eta = eta, mixture_u=u, mixture_v=v, mixture_airspeed=mixture_airspeed, mixture_heading=mixture_heading))
+  }
+
+  output
+}
+
+#' @rdname clean_mixture
+#'
+#' @export
+clean_mixture.vpts <- function(x, ..., slow = 1, fast = 8, drop_slow_component = TRUE, drop_missing = FALSE, keep_mixture = FALSE){
+  assertthat::assert_that(inherits(x,"vpts") | inherits(x,"vp"))
+  assertthat::assert_that(all(c("u_wind","v_wind") %in% names(x$data)), msg="function requires paired wind data, quantities `u_wind` and `v_wind` not found")
+  assertthat::assert_that(assertthat::is.number(x$attributes$how$rcs_bird), msg="radar cross section not defined, please set with `rcs()`")
+
+  # call function
+  result <- clean_mixture.default(x$data$eta, x$data$u, x$data$v, x$data$u_wind, x$data$v_wind,
+                                  slow = slow, fast = fast, drop_slow_component = drop_slow_component,
+                                  drop_missing = drop_missing, keep_mixture = keep_mixture)
+
+  # map results to vp or vpts structure
+  x$data$eta=result$eta
+  x$data$u=result$u
+  x$data$v=result$v
+  x$data$airspeed=result$airspeed
+  x$data$heading=result$heading
+  x$data$airspeed_u=result$airspeed_u
+  x$data$airspeed_v=result$airspeed_v
+  x$data$proportion_slow=result$f
+  if(keep_mixture){
+    x$data$mixture_eta=result$mixture_eta
+    x$data$mixture_u=result$mixture_u
+    x$data$mixture_v=result$mixture_v
+    x$data$mixture_airspeed=result$mixture_airspeed
+    x$data$mixture_heading=result$mixture_heading
+  }
+
+  # trigger rebuild of dens from eta:
+  rcs(x) <- rcs(x)
+
+  # return object
+  x
+}
+
+#' @rdname clean_mixture
+#'
+#' @export
+clean_mixture.vp <- function(x, ..., slow = 1, fast = 8, drop_slow_component = TRUE, drop_missing = FALSE, keep_mixture = FALSE){
+  assertthat::assert_that(inherits(x,"vp"))
+
+  clean_mixture.vpts(x$data$eta, x$data$u, x$data$v, x$data$u_wind, x$data$v_wind,
+                        slow = slow, fast = fast, drop_slow_component = drop_slow_component,
+                        drop_missing = drop_missing, keep_mixture = keep_mixture)
+}

R/filter_precip.R

@@ -0,0 +1,99 @@
+#' Posthoc precipitation filter
+#'
+#' The posthoc precipitation filter assesses how much of the altitude column has a high
+#' total reflectivity factor (biology + meteorology) consistent with precipitation,
+#' and removes biological signals when there is evidence for the presence of precipitation.
+#' Applied to  vertical profiles ('vp') or time series of vertical profiles ('vpts').
+#'
+#' @param x A `vp` or `vpts` object.
+#' @param dbz The minimum reflectivity factor for precipitation.
+#' @param range The minimum altitude range with total reflectivity factor.
+#' `DBZH` > `dbz` at which the filter is triggered.
+#' @param alt_max Maximum altitude above ground level to consider.
+#' @param drop When `TRUE` the profile is removed from the
+#' @return A `vpts` object or a `vp` object, depending on input `x`.
+#'
+#' @export
+#' @details
+#' During precipitation events usually a high proportion of the altitude
+#' column will show a high total reflectivity `DBZH` (which includes biology +
+#' meteorology), because precipitation usually falls from several kilometers high
+#' to the ground surface. Precipitation events are often obvious in profile plots
+#' of quantity `DBZH` as reflectivity signals extending from ground level to high
+#' altitudes far above the typical altitudes where biology is expected. This filter
+#' identifies and removes these cases.
+#'
+#' The posthoc precipitation filter examines the total reflectivity factor `DBZH`
+#' and calculates the altitude range at which `DBZH` is larger than parameter
+#' `dbz`. Whenever this altitude range is larger than parameter `range` (and `drop`
+#' is `FALSE`), the biology is removed by setting profile quantities `dens` and `eta`
+#' to zero and profile quantity `dbz` to `-Inf`. When parameter `drop`
+#' is `TRUE`, the profile is removed from the time series altogether.
+#'
+#' This posthoc precipitation filter is likely to remove biological scatterers
+#' that co-occur with precipitation, for example biological movements during isolated
+#' convective thunderstorms. It is more aggressive than precipitation filters
+#' applied during vertical profile generation with `calculate_vp()` that attempt to
+#' remove precipitation and retain biology. The posthoc precipitation filter is especially
+#' useful for analyses where you want to minimize the risk of precipitation contamination,
+#' at the cost of filtering out some biology during precipitation events.
+#'
+#' Lowering the default minimum reflectivity (`dbz`) for precipitation
+#' below 7 dBZ is not recommended, as most precipitation has a reflectivity above 7 dBZ.
+#'
+#' Parameter `range` should be chosen carefully, and should be higher than the
+#' typical altitude where biological scatterers still reach a reflectivity factor equal to `dbz`.
+#'
+#' @seealso [eta_to_dbz(), dbz_to_eta()]
+#'
+#' @examples
+#' # rain periods are visible in quantity DBZH as dark vertical lines
+#' # extending above 3 km height:
+#' plot(regularize_vpts(example_vpts), quantity='DBZH')
+#' # Apply posthoc filter to remove profiles during precipitation events:
+#' # (regularization is applied to visualize removed profiles)
+#' my_vpts <- regularize_vpts(filter_precip(example_vpts, drop=TRUE))
+#' # verify that rain events have been removed:
+#' plot(my_vpts, quantity='DBZH')
+#' # the posthoc filter removes some biology during precipitation:
+#' plot(my_vpts, quantity='dens')
+#' # original retains more biology (but at the cost of a higher
+#' # likelihood of occasional precipitation contamination):
+#' plot(regularize_vpts(example_vpts), quantity='dens')
+#' # filter can also be applied to single vp objects:
+#' filter_precip(example_vp)
+filter_precip <- function(x, dbz=7, range=2500, alt_max=3000, drop=FALSE){
+  assertthat::assert_that(is.vp(x) | is.vpts(x))
+  assertthat::assert_that(assertthat::is.number(dbz))
+  assertthat::assert_that(assertthat::is.number(range))
+  assertthat::assert_that(assertthat::is.number(alt_max))
+  assertthat::assert_that(alt_max>0)
+  assertthat::assert_that(assertthat::is.flag(drop))
+  assertthat::assert_that(range<=alt_max)
+  assertthat::assert_that(!(drop & is.vp(x)), msg='parameter `drop` should be `TRUE` for objects of class `vp`')
+  if(dbz<7) warning("dbz value too low for typical precipitation")
+  height_index_max <- ((x$attributes$where$height + alt_max) %/% x$attributes$where$interval)
+  height_index_max <- min(x$attributes$where$levels,height_index_max)
+  if(is.vpts(x)){
+    height_range <- colSums(x$data$DBZH[1:height_index_max,]>dbz,na.rm=T)*x$attributes$where$interval
+  } else{
+    height_range <- sum(x$data$DBZH[1:height_index_max]>dbz,na.rm=T)*x$attributes$where$interval
+  }
+
+  index <- which(height_range > range)
+  if(length(index)==0) return(x)
+  # if remove, drop the profiles
+  if(drop) return(x[-index])
+  # otherwise set the density field to zero, but keep the profile
+  if(is.vpts(x)){
+    x$data$dens[,index] <- 0
+    x$data$eta[,index] <- 0
+    x$data$dbz[,index] <- -Inf
+  } else{
+    x$data$dens[index] <- 0
+    x$data$eta[index] <- 0
+    x$data$dbz[index] <- -Inf
+  }
+
+  x
+}