add vertical Integration clean_mixture() quantities

NEWS.md

@@ -1,11 +1,18 @@
 # bioRad 0.12.0.9000
 
-* `apply_mistnet()` now accepts `pvol` objects in addition to pvolfiles.
+## New features
+* `integrate_profile()` now includes quantities added by `clean_mixture()` (#785).
 
-* bugfix `calculate_param()` when using `ifelse` (#770).
+* `apply_mistnet()` now accepts `pvol` objects in addition to pvolfiles (#788).
 
 * Improve axis labels for `map()` and provide option not to plot the radar (#773).
 
+* `clean_mixture()` argument `drop_missing` is now `TRUE` by default (#787).
+
+## Bugfixes
+
+* bugfix `calculate_param()` when using `ifelse` (#770).
+
 * bugfix `integrate_profile()` for height-integrated speed quantities (`u`,`v`,`ff`) in vpts objects (#782).
 
 # bioRad 0.12.0 

R/clean_mixture.R

@@ -21,7 +21,7 @@
 #' @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).
+#' and wind speed are set to NaN when `TRUE` (default), or returned unaltered when `FALSE`.
 #' @param keep_mixture When `TRUE` store original mixture reflectivity and speeds as
 #' renamed quantities with `mixture_` prefix
 #' @returns a named list with cleaned densities and speeds.
@@ -116,7 +116,7 @@ clean_mixture <- function(x, ...){
 
 #' @rdname clean_mixture
 #' @export
-clean_mixture.default <- function(x, slow = 1, fast = 8, drop_slow_component = TRUE, drop_missing = FALSE, keep_mixture = FALSE, u_wind, v_wind, u, v, ...){
+clean_mixture.default <- function(x, slow = 1, fast = 8, drop_slow_component = TRUE, drop_missing = TRUE, keep_mixture = FALSE, u_wind, v_wind, u, v, ...){
   # verify input
   assertthat::assert_that(all(x >= 0, na.rm = TRUE))
   assertthat::assert_that(is.numeric(u))
@@ -183,9 +183,12 @@ clean_mixture.default <- function(x, slow = 1, fast = 8, drop_slow_component = T
   }
 
   if(drop_missing){
-    eta_corr[is.na(f)]=NaN
-    air_u[is.na(f)]=NaN
-    air_v[is.na(f)]=NaN
+    eta_corr[is.na(f)]=NA
+    eta_corr[is.nan(f)]=NaN
+    air_u[is.na(f)]=NA
+    air_u[is.nan(f)]=NaN
+    air_v[is.na(f)]=NA
+    air_v[is.nan(f)]=NaN
   }
 
   # calculate speed and direction
@@ -205,7 +208,7 @@ clean_mixture.default <- function(x, slow = 1, fast = 8, drop_slow_component = T
 
 #' @rdname clean_mixture
 #' @export
-clean_mixture.vpts <- function(x, slow = 1, fast = 8, drop_slow_component = TRUE, drop_missing = FALSE, keep_mixture = FALSE, u_wind="u_wind", v_wind="v_wind", ...){
+clean_mixture.vpts <- function(x, slow = 1, fast = 8, drop_slow_component = TRUE, drop_missing = TRUE, keep_mixture = FALSE, u_wind="u_wind", v_wind="v_wind", ...){
   assertthat::assert_that(inherits(x,"vpts") | inherits(x,"vp"))
   if(inherits(x,"vpts") | inherits(x,"vp")){
     assertthat::assert_that(is.character(u_wind))
@@ -267,9 +270,9 @@ clean_mixture.vpts <- function(x, slow = 1, fast = 8, drop_slow_component = TRUE
   if(sum(presence_test)>0) warning(paste0("Overwriting existing quantities `", paste(quantities[presence_test], collapse="`, `"),"`."))
 
   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$heading=result$heading
   x$data$f=result$f
   if(keep_mixture){
     x$data$mixture_eta=result$mixture_eta
@@ -288,7 +291,7 @@ clean_mixture.vpts <- function(x, slow = 1, fast = 8, drop_slow_component = TRUE
 
 #' @rdname clean_mixture
 #' @export
-clean_mixture.vp <- function(x, ..., slow = 1, fast = 8, drop_slow_component = TRUE, drop_missing = FALSE, keep_mixture = FALSE, u_wind="u_wind", v_wind="v_wind"){
+clean_mixture.vp <- function(x, ..., slow = 1, fast = 8, drop_slow_component = TRUE, drop_missing = TRUE, keep_mixture = FALSE, u_wind="u_wind", v_wind="v_wind"){
   assertthat::assert_that(inherits(x,"vp"))
 
   clean_mixture.vpts(x,slow = slow, fast = fast, drop_slow_component = drop_slow_component,

R/integrate_profile.R

@@ -284,10 +284,10 @@ integrate_profile.vp <- function(x, alt_min = 0, alt_max = Inf, alpha = NA,
     airspeed_u <- get_quantity(x, "u") - get_quantity(x, "u_wind")
     airspeed_v <- get_quantity(x, "v") - get_quantity(x, "v_wind")
     output$airspeed <- stats::weighted.mean(sqrt(airspeed_u^2 + airspeed_v^2), weight_densdh, na.rm = TRUE)
-    output$heading <- stats::weighted.mean((pi / 2 - atan2(airspeed_v, airspeed_u)) * 180 / pi, weight_densdh, na.rm = TRUE)
-    output$heading[which(output$heading<0)]=output$heading[which(output$heading<0)]+360
     output$airspeed_u <- stats::weighted.mean(airspeed_u, weight_densdh, na.rm = TRUE)
     output$airspeed_v <- stats::weighted.mean(airspeed_v, weight_densdh, na.rm = TRUE)
+    output$heading <- (pi / 2 - atan2(output$airspeed_v, output$airspeed_u)) * 180 / pi
+    output$heading[which(output$heading<0)]=output$heading[which(output$heading<0)]+360
     output$ff_wind <- stats::weighted.mean(sqrt(get_quantity(x,"u_wind")^2 + get_quantity(x,"v_wind")^2), weight_densdh, na.rm = TRUE)
     output$u_wind <- stats::weighted.mean(get_quantity(x,"u_wind"), weight_densdh, na.rm = TRUE)
     output$v_wind <- stats::weighted.mean(get_quantity(x,"v_wind"), weight_densdh, na.rm = TRUE)
@@ -387,7 +387,7 @@ integrate_profile.vpts <- function(x, alt_min = 0, alt_max = Inf,
   dh <- pmin(pmin(pmax(x$height+x$attributes$where$interval-alt_min,0),interval),
              pmin(pmax(alt_max-x$height,0),interval)) / 1000
 
-  # helper function that performs a weighted colSum, and that will return
+  # helper function that performs a standard colSum, but returning
   # NaN if all elements in a column are either NA or NaN.
   nan_colSums <- function(x, ...) {
      s <- colSums(x, na.rm = TRUE, ...)
@@ -420,7 +420,7 @@ integrate_profile.vpts <- function(x, alt_min = 0, alt_max = Inf,
   weight_densdh[is.na(weight_densdh)] <- 0
   # Normalize the weight of each vp by its column sum.
   weight_densdh <- sweep(weight_densdh, 2, colSums(weight_densdh), FUN="/")
-  # Find index where no bird are present
+  # Find index where no bird are present - previous line introduces NaN when dividing by zero
   no_bird <- is.na(colSums(weight_densdh))
 
   # create a separate weighting matrix for speed quantities
@@ -487,14 +487,29 @@ integrate_profile.vpts <- function(x, alt_min = 0, alt_max = Inf,
     airspeed_u <- get_quantity(x, "u") - get_quantity(x, "u_wind")
     airspeed_v <- get_quantity(x, "v") - get_quantity(x, "v_wind")
     output$airspeed <- nan_colSums(sqrt(airspeed_u^2 + airspeed_v^2) * weight_ffdh)
-    output$heading <- nan_colSums(((pi / 2 - atan2(airspeed_v, airspeed_u)) * 180 / pi) * weight_ffdh)
-    output$heading[which(output$heading<0)]=output$heading[which(output$heading<0)]+360
     output$airspeed_u <- nan_colSums(airspeed_u * weight_ffdh)
     output$airspeed_v <- nan_colSums(airspeed_v * weight_ffdh)
+    output$heading <- (pi / 2 - atan2(output$airspeed_v, output$airspeed_u)) * 180 / pi
+    output$heading[which(output$heading<0)]=output$heading[which(output$heading<0)]+360
     output$ff_wind <- nan_colSums(sqrt(get_quantity(x,"u_wind")^2 + get_quantity(x,"v_wind")^2) * weight_densdh)
     output$u_wind <- nan_colSums(get_quantity(x,"u_wind") * weight_densdh)
     output$v_wind <- nan_colSums(get_quantity(x,"v_wind") * weight_densdh)
   }
+  if("mixture_eta" %in% names(x$data)){
+    output$mixture_vir <- colSums(get_quantity(x, "mixture_eta") * dh, na.rm = TRUE)
+  }
+  quantities <- c("mixture_u", "mixture_v", "mixture_airspeed", "mixture_heading")
+  present <- quantities[quantities %in% names(x$data)]
+  output[present] <- lapply(present, function(q) {
+    nan_colSums(get_quantity(x, q) * weight_ffdh)
+  })
+  if(all(c("f", "mixture_eta") %in% names(x$data))){
+    eta_slow <- nan_colSums(get_quantity(x, "f") * get_quantity(x, "mixture_eta") * (!is.na(get_quantity(x, "ff"))) * dh)
+    # mixture_vir above includes layers without speeds,
+    # eta_mixture below excludes layers without speeds:
+    eta_mixture <- colSums(get_quantity(x, "mixture_eta") * (!is.na(get_quantity(x, "ff"))) * dh, na.rm = TRUE)
+    output$f <- eta_slow/eta_mixture
+  }
 
   class(output) <- c("vpi", "data.frame")
   rownames(output) <- NULL