diff --git a/NEWS.md b/NEWS.md index 904c81a85..f5cade654 100644 --- a/NEWS.md +++ b/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) diff --git a/R/bind_into_vpts.R b/R/bind_into_vpts.R index 229e7f32c..1b2ffa9da 100644 --- a/R/bind_into_vpts.R +++ b/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( diff --git a/man/bind_into_vpts.Rd b/man/bind_into_vpts.Rd index 28898565b..b330d6720 100644 --- a/man/bind_into_vpts.Rd +++ b/man/bind_into_vpts.Rd @@ -33,6 +33,11 @@ Binds vertical profiles (\code{vp}) into a vertical profile time series (\code{vpts}), sorted in time. Can also bind multiple \code{vpts} of a single radar into one \code{vpts}. } +\details{ +\code{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 +} \section{Methods (by class)}{ \itemize{ \item \code{vp}: Bind multiple \code{vp} into a \code{vpts}. diff --git a/tests/testthat/test-bind_into_vpts.R b/tests/testthat/test-bind_into_vpts.R index db59d667a..a9878c0fb 100644 --- a/tests/testthat/test-bind_into_vpts.R +++ b/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) +})