bind profiles with different number of altitude layers (#343)

NEWS.md

@@ -2,6 +2,8 @@
 
 * Warn when more then one scan could be returned (#414) and provide alternative
 
+* `bind_into_vpts()` now works for vp's and vpts's with different heights (#343)
+
 * bugfix ylim argument in `composite_ppi()` (#389)
 
 * `scan_to_spatial()` now creates points for cell centers (#430)

R/bind_into_vpts.R

@@ -13,6 +13,10 @@
 #'
 #' @export
 #'
+#' @details `bind_into_vpts()` currently requires profiles to have aligning altitude
+#' layers that are of equal width. Profiles are allowed to differ in the number
+#' of altitude layers, i.e. the maximum altitude
+#'
 #' @examples
 #' # load example time series of vertical profiles:
 #' data(example_vpts)
@@ -86,9 +90,12 @@ bind_into_vpts.vpts <- function(..., attributes_from = 1) {
   if (length(radars) > 1) {
     stop("Vertical profiles are not from a single radar")
   }
-  if (length(unique(lapply(vptss, "[[", "height"))) > 1) {
-    stop("Vertical profiles have non-aligning altitude layers")
+  height_list <- lapply(vptss, "[[", "height")
+  if (length(unique(height_list)) > 1) {
+    target_heights <- combined_heights(height_list)
+    vptss <- lapply(vptss, add_heights_vpts, target = target_heights)
   }
+
   if (length(unique(lapply(vptss, function(x) names(x$"data")))) > 1) {
     stop("Vertical profiles have different quantities")
   }
@@ -172,7 +179,41 @@ vplist_to_vpts <- function(x, radar = NA) {
     ))
   }
 }
-
+combined_heights <- function(x) {
+  assert_that(is.list(x))
+  unique_height_diff <- unique(unlist(lapply(lapply(x, diff), unique)))
+  if (length(unique_height_diff) != 1) {
+    stop("Not all data has the same size of altitude bins")
+  }
+  height_alignment <- unique(unlist(x) %% unique_height_diff)
+  if (length(height_alignment) != 1) {
+    stop("Not all data has the same alignment of altitude bins")
+  }
+  return(sort(unique(unlist(x))))
+}
+add_heights_vpts <- function(x, target) {
+  if (identical(x$height, target)) {
+    return(x)
+  }
+  old <- target %in% x$height
+  x$height <- target
+  x$attributes$where$levels <- length(target)
+  m <- matrix(nrow = length(target), ncol = length(x$datetime))
+  x$data <- lapply(x$data, function(x, m, s) {
+    m[s, ] <- x
+    return(m)
+  }, s = old, m = m)
+  return(x)
+}
+add_heights_vp <- function(x, target) {
+  if (identical(x$data$height, target)) {
+    return(x)
+  }
+  x$data <- data.frame(rbindlist(list(x$data, data.frame(height = target[!(target %in% x$data$height)])), fill = TRUE))
+  x$data <- x$data[order(x$data$height), ]
+  x$attributes$where$levels <- length(target)
+  return(x)
+}
 vp_to_vpts_helper <- function(vps) {
   datetime <- .POSIXct(do.call("c", lapply(vps, "[[", "datetime")), tz = "UTC")
   daterange <- .POSIXct(c(min(datetime), max(datetime)), tz = "UTC")
@@ -181,11 +222,10 @@ vp_to_vpts_helper <- function(vps) {
   datetime <- .POSIXct(do.call("c", lapply(vps, "[[", "datetime")), tz = "UTC")
   difftimes <- difftime(datetime[-1], datetime[-length(datetime)], units = "secs")
   profile.quantities <- names(vps[[1]]$data)
-  if (length(unique(lapply(vps, function(x) x$data$height))) > 1) {
-    stop(paste(
-      "Vertical profiles of radar", vps[[1]]$radar,
-      "have non-aligning altitude layers."
-    ))
+  height_list <- lapply(vps, function(x) x$data$height)
+  if (length(unique(height_list)) > 1) {
+    target_heights <- combined_heights(height_list)
+    vps <- lapply(vps, add_heights_vp, target = target_heights)
   }
   if (length(unique(lapply(vps, function(x) names(x$"data")))) > 1) {
     stop(paste(

tests/testthat/test-bind_into_vpts.R

@@ -14,3 +14,32 @@ test_that("vpts based on vpts or multiple vpts", {
   expect_s3_class(bind_into_vpts(example_vpts[1:10]), "vpts")
   expect_s3_class(bind_into_vpts(example_vpts[1:10], example_vpts[11:20]), "vpts")
 })
+
+test_that("vp with different heights", {
+  data(example_vp)
+  lower <- example_vp
+  lower$data <- lower$data[1:20, ]
+  lower$attributes$where$levels = 20
+  higher <- example_vp
+  higher$data <- higher$data[3:25, ]
+  lower$attributes$where$levels = 23
+  expect_s3_class(bind_into_vpts(example_vp, lower), "vpts")
+  expect_s3_class(bind_into_vpts(list(example_vp, higher)), "vpts")
+  expect_s3_class(bind_into_vpts(example_vp, higher), "vpts")
+  vpts <- bind_into_vpts(lower, higher)
+  expect_s3_class(vpts, "vpts")
+  expect_equal(vpts$attributes$where$levels, example_vpts$attributes$where$levels)
+  higher$data$height <- higher$data$height + 1
+  expect_error(bind_into_vpts(example_vp, higher))
+})
+
+test_that("vpts with different heights", {
+  data(example_vpts)
+  lower <- example_vpts
+  lower$height <- lower$height[1:20]
+  lower$attributes$where$levels = 20
+  lower$data <- lapply(lower$data, head, 20)
+  vpts <- bind_into_vpts(lower, example_vpts)
+  expect_s3_class(vpts,'vpts')
+  expect_equal(vpts$attributes$where$levels, example_vpts$attributes$where$levels)
+})