diff --git a/NAMESPACE b/NAMESPACE index b3d6821c..f3cf51b4 100644 --- a/NAMESPACE +++ b/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) @@ -86,6 +89,7 @@ S3method(summary,vpts) export("rcs<-") export("sd_vvp_threshold<-") export(apply_mistnet) +export(as.vp) export(as.vpts) export(attribute_table) export(beam_distance) @@ -99,6 +103,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 +112,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) diff --git a/NEWS.md b/NEWS.md index acf878c5..04be4d85 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,11 +1,20 @@ # bioRad 0.8.1.9000 +## new features -* Correct units specified in plot label for quantity VIR (#674). +* new function `as.vp()` to convert data.frame to vertical profile object (#699) -* Discard profiles with misspecified altitude bins in `as.vpts()` and `read_vpts()` (#684). +* new function `clean_mixture(`) for partitioning bird and insect mixtures (#700) + +* new function `filter_precip()` for posthoc removal of precipitation in vp and vpts objects (#701) * New argument `zoomin` for function `bioRad::map()` to increase basemap resolution (#689). +## bugfixes + +* Corrected units specified in plot label for quantity VIR (#674). + +* Discard profiles with misspecified altitude bins in `as.vpts()` and `read_vpts()` (#684). + * Correct type of gap field in vpts objects for profiles stored in ODIM HDF5 format (.h5) (#635, #691). * Interpret NA values in field DBZH in `integrate_to_ppi()` as pixels that were not irradiated (#658). diff --git a/R/as.vp.R b/R/as.vp.R new file mode 100644 index 00000000..146a9b55 --- /dev/null +++ b/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) +} diff --git a/R/as.vpts.R b/R/as.vpts.R index 9562678e..3530b210 100644 --- a/R/as.vpts.R +++ b/R/as.vpts.R @@ -16,13 +16,14 @@ as.vpts <- function(data) { data <- data %>% dplyr::rename(DBZH = "dbz_all") } - + validate_vpts(data) # sort by datetime and height data <- dplyr::arrange(data, datetime, height) height <- datetime <- source_file <- radar <- NULL + bioRad_internal_levels <- bioRad_internal_interval <- NULL # Throw error if nrows per height are not identical # FIXME: first if statement is a weak check that could fail, could be improved. @@ -30,16 +31,16 @@ as.vpts <- function(data) { if(!remainder_is_zero(dim(data)[1], length(unique(data$height)))){ data %>% dplyr::group_by(radar, datetime) %>% - dplyr::mutate(bioRad_internal_interval = height-lag(height)) %>% + dplyr::mutate(bioRad_internal_interval = height-dplyr::lag(height)) %>% dplyr::add_count(name="bioRad_internal_levels") -> data - interval_median <- median(data$bioRad_internal_interval, na.rm=TRUE) + interval_median <- stats::median(data$bioRad_internal_interval, na.rm=TRUE) interval_unique <- unique(data$bioRad_internal_interval) interval_unique <- interval_unique[!is.na(interval_unique)] if(length(interval_unique)>1){ warning(paste("profiles found with different altitude interval:",paste(sort(interval_unique),collapse=" ")), ", retaining ",interval_median, " only.") data <- dplyr::filter(data, bioRad_internal_interval == interval_median) } - levels_median <- median(data$bioRad_internal_levels) + levels_median <- stats::median(data$bioRad_internal_levels) levels_unique <- unique(data$bioRad_internal_levels) if(length(levels_unique)>1){ warning(paste("profiles found with different number of height layers:",paste(sort(levels_unique),collapse=" ")), ", retaining ",levels_median, " only.") diff --git a/R/clean_mixture.R b/R/clean_mixture.R new file mode 100644 index 00000000..7184595b --- /dev/null +++ b/R/clean_mixture.R @@ -0,0 +1,252 @@ +#' Partition mixtures of animals using assumptions on airspeeds. +#' +#' Partition mixtures of birds and insects using assumptions on their respective airspeeds, +#' following the approach by Shi et al. (2025). +#' @param x a `vp` or `vpts` object, or a mixture animal density or linear reflectivity eta in cm\eqn{^2}/km\eqn{^3}. +#' @param ... `eta`, `u`, `v`, `U`, `V` arguments, taken from object for `vp` or `vpts` class. +#' @param u the mixture's ground speed u component (west to east) in m/s. +#' @param v the mixture's ground speed v component (south to north) in m/s. +#' @param U the west to east wind component in m/s. In the case of `vp` and `vpts` objects +#' the quantity name for the U-component of the wind. +#' @param V the south to north wind component in m/s. In the case of `vp` and `vpts` objects +#' the quantity name for the V-component of the wind. +#' @param fast the fast component's airspeed in m/s, typically the airspeed of insects. +#' @param slow the slow component's airspeed in m/s, typically the airspeed of birds. +#' @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 named list with cleaned densities and speeds. +#' Output differs depending on whether the fast component is retained +#' (`drop_slow_component`=`TRUE`, default) or the slow component (`drop_slow_component`=`FALSE`, default). Output quantities include: +#' * `eta`: cleaned reflectivity in cm^2/km^3. +#' only the fast component (default) or the slow component +#' (when `drop_slow_component` is `TRUE`). +#' * `u`: cleaned ground speed component west to east in m/s. +#' * `v`: cleaned ground speed component south to north in m/s. +#' * `airspeed`: the airspeed of the selected component in m/s. +#' * `airspeed_u`: the u-component (west to east) of the airspeed of the retained component in m/s. +#' * `airspeed_v`: the v-component (south to north) of the airspeed of the retained component in m/s. +#' * `heading`: the heading of the selected component in degrees clockwise from north. +#' * `f`: the reflectivity proportion of the slow component (0-1 range), typically the proportion of insects. +#' +#' For `vp` and `vpts` objects the quantities `eta`,`u`,`v` will be updated, and other +#' quantities listed above will be added. +#' @export +#' @name clean_mixture +#' @examples +#' # convert profile object to data.frame +#' df <- as.data.frame(example_vp, suntime=FALSE) +#' # add wind u and v component wind data +#' # (here a NW wind identical at all altitudes) +#' df$u_wind=3 +#' df$v_wind=-3 +#' # convert back to vp object +#' my_vp <- as.vp(df) +#' # partition the mixture: +#' my_vp_clean <- clean_mixture(my_vp) +#' +#' # drop the slow component (typically insects) +#' clean_mixture(100,u=-13,v=13,U=-7,V=6, fast=8, slow=1) +#' # drop the fast component (typically birds) +#' clean_mixture(100,u=-13,v=13,U=-7,V=6, fast=8, slow=1, drop_slow_component=FALSE) +#' # keep the original mixture reflectivity and speed components +#' clean_mixture(100,u=-13,v=13,U=-7,V=6, fast=8, slow=1, keep_mixture=TRUE) +#' # keep reflectivity unaltered when one of the speed components is not a number: +#' clean_mixture(100,u=-13,v=13,U=NaN,V=6, fast=8, slow=1)["eta"] +#' # set reflectivity to NaN when one of the speed components is not a number: +#' clean_mixture(100,u=-13,v=13,U=NaN,V=6, fast=8, slow=1, drop_missing=TRUE)["eta"] +#' @references +#' * Shi X, Drucker J, Chapman JW, Sanchez Herrera M, Dokter AM +#' Analysis of mixtures of birds and insects in weather radar data. +#' Ornithological Applications. 2025 (in press). +#' * Nussbaumer R, Schmid B, Bauer S, Liechti F. +#' A Gaussian mixture model to separate birds and insects in +#' single-polarization weather radar data. Remote Sensing. 2021 May 19;13(10):1989. +#' @details +#' For a detail description of the methodology see Shi et al. (2025). +#' Most commonly the fast component refers to migrating birds, while +#' the slow component refers to insects. The slow component is always +#' oriented in the direction of the wind by definition. +#' Note that for mixture airspeeds exceeding the airspeed of the fast component, +#' all reflectivity is assigned to the fast component. Similarly, for mixture +#' airspeeds below the airspeed of the slow component, all reflectivity +#' will be assigned to the slow component. In both these limiting cases the +#' +#' +#' ## How to use this function? +#' 1. To apply this function to `vp` or `vpts` data altitudinal wind data +#' needs to be added to the vertical profile data first. This is most easily +#' accomplished by first converting the objects to a `data.frame` with [as.vp()] or [as.vpts()]. +#' Wind data can then be added as a new columns to the data.frame. +#' By default the wind data is expected to be named `u_wind` for the U component +#' and `v_wind` for the V component of the wind. +#' Alternatively, arguments `U` and `V` can be used to specify different names. +#' 2. Realistic assumptions for the expected airspeed for the slow (insect) +#' and fast (bird) components need to be provided, using arguments `slow` and `fast`. +#' See Shi et al. 2025 for recommendations in choosing these values. +#' 3. Use `drop_slow_component` to toggle between retaining the slow or fast component. +#' When `TRUE` the fast (bird) component is retained. When `FALSE` the slow (insect) +#' component is retained. Note that in this case the corrected ground speed direction will be +#' identical to the wind direction, and the magnitude of the ground speed will be equal +#' to the wind speed plus the value of `slow`, due to the underlying assumption +#' of wind following by the slow component. +NULL + +#' @rdname clean_mixture +#' +#' @export +clean_mixture <- function(x, ...){ + UseMethod("clean_mixture", 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, V, u, v, ...){ + # verify input + assertthat::assert_that(all(x >= 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(is.numeric(slow)) + assertthat::assert_that(is.numeric(fast)) + 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=x + 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)*x)[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*x)[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 = x, 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, U="u_wind", V="v_wind", ...){ + assertthat::assert_that(inherits(x,"vpts") | inherits(x,"vp")) + if(inherits(x,"vpts") | inherits(x,"vp")){ + assertthat::assert_that(is.character(U)) + assertthat::assert_that(is.character(V)) + } + assertthat::assert_that(all(c(U,V) %in% names(x$data)), msg=paste0("function requires paired wind data, profile quantities `",U,"` and/or `",V,"` not found.")) + assertthat::assert_that(all(c("u","v") %in% names(x$data)), msg="function requires ground speed data, profile quantities `u` and/or `v` not found.") + assertthat::assert_that("eta" %in% names(x$data), msg="function requires linear reflectivity data, profile quantity `eta` 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, + slow = slow, fast = fast, drop_slow_component = drop_slow_component, + drop_missing = drop_missing, keep_mixture = keep_mixture, + x$data[[U]], x$data[[V]], x$data$u, x$data$v) + + # map results to vp or vpts structure + x$data$eta=result$eta + x$data$u=result$u + x$data$v=result$v + + # check for quantities potentially overwritten + quantities <- c("airspeed","airspeed_u","airspeed_v","heading","f") + if(keep_mixture) quantities <- c(quantities, "mixture_eta","mixture_u","mixture_v","mixture_airspeed","mixture_heading") + presence_test <- quantities %in% names(x$data) + 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$f=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, U="u_wind", V="v_wind"){ + assertthat::assert_that(inherits(x,"vp")) + + clean_mixture.vpts(x,slow = slow, fast = fast, drop_slow_component = drop_slow_component, + drop_missing = drop_missing, keep_mixture = keep_mixture, U=U, V=V) +} diff --git a/R/filter_precip.R b/R/filter_precip.R new file mode 100644 index 00000000..5b0e37ae --- /dev/null +++ b/R/filter_precip.R @@ -0,0 +1,105 @@ +#' 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`. +#' +#' Note that at S-band wavelengths bird migration occurs much more frequently in the reflectivity +#' regime for precipitation than at C-band. Therefore, at C-band lower settings for parameter `dbz` +#' are appropriate than at S-band. +#' +#' @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 +} diff --git a/_pkgdown.yml b/_pkgdown.yml index d6176aba..0f94e415 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -25,7 +25,7 @@ home: reference: - title: "Reading radar data" - desc: "Functions to read polar volume (pvol) data. Requires a locally running Docker daemon for NEXRAD data." + desc: "Functions to read polar volume (pvol) data." contents: - download_pvolfiles - is.pvolfile @@ -67,7 +67,7 @@ reference: - composite_ppi - "`[.ppi`" - title: "Creating vertical profiles of biological targets" - desc: "Functions to process weather radar data (pvol) into vertical profiles (vp) of biological targets. Requires a locally running Docker daemon." + desc: "Functions to process weather radar data (pvol) into vertical profiles (vp) of biological targets." contents: - calculate_vp - nexrad_to_odim # Repeated @@ -84,18 +84,21 @@ reference: - plot.vp - as.data.frame.vp - list_vpts_aloft - - title: "Combining vertical profiles into time series" - desc: "Functions to combine vertical profiles (vp) into time series (vpts) and read, inspect and plot these." + - title: "Manipulating vertical profile data" + desc: "Functions to combine vertical profiles (vp) into time series (vpts) and to post-process, read, inspect and plot these." contents: - c.vp - bind_into_vpts - read_vpts - filter_vpts + - filter_precip + - clean_mixture - summary.vpts # is.vpts - example_vpts - plot.vpts - regularize_vpts - as.data.frame.vp + - as.vp - as.vpts - "`[.vpts`" - get_quantity @@ -116,7 +119,7 @@ reference: - doy_noy - attribute_table - title: "Radar beam geometry" - desc: "Functions relating the radar beam shape to range, distance and height" + desc: "Functions relating the radar beam shape to range, distance and height." contents: - beam_height - beam_width diff --git a/bioRad.Rproj b/bioRad.Rproj index 0dfffde6..87e26f28 100644 --- a/bioRad.Rproj +++ b/bioRad.Rproj @@ -1,5 +1,5 @@ Version: 1.0 -ProjectId: b32cc17d-04d6-4782-95b0-7bdaa77d889b +ProjectId: 78f5cd79-a10d-46e0-8b29-c679c1dee61d RestoreWorkspace: Default SaveWorkspace: Default diff --git a/data/example_vp.rda b/data/example_vp.rda index 983c87d6..9e365471 100644 Binary files a/data/example_vp.rda and b/data/example_vp.rda differ diff --git a/man/as.vp.Rd b/man/as.vp.Rd new file mode 100644 index 00000000..2159e998 --- /dev/null +++ b/man/as.vp.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/as.vp.R +\name{as.vp} +\alias{as.vp} +\title{Convert a dataframe into a vp object} +\usage{ +as.vp(data) +} +\arguments{ +\item{data}{a dataframe created from a VPTS CSV file} +} +\value{ +a bioRad vp object +} +\description{ +Convert a dataframe into a 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) +} diff --git a/man/clean_mixture.Rd b/man/clean_mixture.Rd new file mode 100644 index 00000000..371af80e --- /dev/null +++ b/man/clean_mixture.Rd @@ -0,0 +1,164 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/clean_mixture.R +\name{clean_mixture} +\alias{clean_mixture} +\alias{clean_mixture.default} +\alias{clean_mixture.vpts} +\alias{clean_mixture.vp} +\title{Partition mixtures of animals using assumptions on airspeeds.} +\usage{ +clean_mixture(x, ...) + +\method{clean_mixture}{default}( + x, + slow = 1, + fast = 8, + drop_slow_component = TRUE, + drop_missing = FALSE, + keep_mixture = FALSE, + U, + V, + u, + v, + ... +) + +\method{clean_mixture}{vpts}( + x, + slow = 1, + fast = 8, + drop_slow_component = TRUE, + drop_missing = FALSE, + keep_mixture = FALSE, + U = "u_wind", + V = "v_wind", + ... +) + +\method{clean_mixture}{vp}( + x, + ..., + slow = 1, + fast = 8, + drop_slow_component = TRUE, + drop_missing = FALSE, + keep_mixture = FALSE, + U = "u_wind", + V = "v_wind" +) +} +\arguments{ +\item{x}{a \code{vp} or \code{vpts} object, or a mixture animal density or linear reflectivity eta in cm\eqn{^2}/km\eqn{^3}.} + +\item{...}{\code{eta}, \code{u}, \code{v}, \code{U}, \code{V} arguments, taken from object for \code{vp} or \code{vpts} class.} + +\item{slow}{the slow component's airspeed in m/s, typically the airspeed of birds.} + +\item{fast}{the fast component's airspeed in m/s, typically the airspeed of insects.} + +\item{drop_slow_component}{when TRUE (default) output density, ground speed and +heading for fast component, when FALSE for slow component.} + +\item{drop_missing}{Values \code{eta} without an associated ground speed +and wind speed are set to NA when \code{TRUE}, or returned unaltered when \code{FALSE} (default).} + +\item{keep_mixture}{When \code{TRUE} store original mixture reflectivity and speeds as +renamed quantities with \code{mixture_} prefix} + +\item{U}{the west to east wind component in m/s. In the case of \code{vp} and \code{vpts} objects +the quantity name for the U-component of the wind.} + +\item{V}{the south to north wind component in m/s. In the case of \code{vp} and \code{vpts} objects +the quantity name for the V-component of the wind.} + +\item{u}{the mixture's ground speed u component (west to east) in m/s.} + +\item{v}{the mixture's ground speed v component (south to north) in m/s.} +} +\value{ +a named list with cleaned densities and speeds. +Output differs depending on whether the fast component is retained +(\code{drop_slow_component}=\code{TRUE}, default) or the slow component (\code{drop_slow_component}=\code{FALSE}, default). Output quantities include: +\itemize{ +\item \code{eta}: cleaned reflectivity in cm^2/km^3. +only the fast component (default) or the slow component +(when \code{drop_slow_component} is \code{TRUE}). +\item \code{u}: cleaned ground speed component west to east in m/s. +\item \code{v}: cleaned ground speed component south to north in m/s. +\item \code{airspeed}: the airspeed of the selected component in m/s. +\item \code{airspeed_u}: the u-component (west to east) of the airspeed of the retained component in m/s. +\item \code{airspeed_v}: the v-component (south to north) of the airspeed of the retained component in m/s. +\item \code{heading}: the heading of the selected component in degrees clockwise from north. +\item \code{f}: the reflectivity proportion of the slow component (0-1 range), typically the proportion of insects. +} + +For \code{vp} and \code{vpts} objects the quantities \code{eta},\code{u},\code{v} will be updated, and other +quantities listed above will be added. +} +\description{ +Partition mixtures of birds and insects using assumptions on their respective airspeeds, +following the approach by Shi et al. (2025). +} +\details{ +For a detail description of the methodology see Shi et al. (2025). +Most commonly the fast component refers to migrating birds, while +the slow component refers to insects. The slow component is always +oriented in the direction of the wind by definition. +Note that for mixture airspeeds exceeding the airspeed of the fast component, +all reflectivity is assigned to the fast component. Similarly, for mixture +airspeeds below the airspeed of the slow component, all reflectivity +will be assigned to the slow component. In both these limiting cases the +\subsection{How to use this function?}{ +\enumerate{ +\item To apply this function to \code{vp} or \code{vpts} data altitudinal wind data +needs to be added to the vertical profile data first. This is most easily +accomplished by first converting the objects to a \code{data.frame} with \code{\link[=as.vp]{as.vp()}} or \code{\link[=as.vpts]{as.vpts()}}. +Wind data can then be added as a new columns to the data.frame. +By default the wind data is expected to be named \code{u_wind} for the U component +and \code{v_wind} for the V component of the wind. +Alternatively, arguments \code{U} and \code{V} can be used to specify different names. +\item Realistic assumptions for the expected airspeed for the slow (insect) +and fast (bird) components need to be provided, using arguments \code{slow} and \code{fast}. +See Shi et al. 2025 for recommendations in choosing these values. +\item Use \code{drop_slow_component} to toggle between retaining the slow or fast component. +When \code{TRUE} the fast (bird) component is retained. When \code{FALSE} the slow (insect) +component is retained. Note that in this case the corrected ground speed direction will be +identical to the wind direction, and the magnitude of the ground speed will be equal +to the wind speed plus the value of \code{slow}, due to the underlying assumption +of wind following by the slow component. +} +} +} +\examples{ +# convert profile object to data.frame +df <- as.data.frame(example_vp, suntime=FALSE) +# add wind u and v component wind data +# (here a NW wind identical at all altitudes) +df$u_wind=3 +df$v_wind=-3 +# convert back to vp object +my_vp <- as.vp(df) +# partition the mixture: +my_vp_clean <- clean_mixture(my_vp) + +# drop the slow component (typically insects) +clean_mixture(100,u=-13,v=13,U=-7,V=6, fast=8, slow=1) +# drop the fast component (typically birds) +clean_mixture(100,u=-13,v=13,U=-7,V=6, fast=8, slow=1, drop_slow_component=FALSE) +# keep the original mixture reflectivity and speed components +clean_mixture(100,u=-13,v=13,U=-7,V=6, fast=8, slow=1, keep_mixture=TRUE) +# keep reflectivity unaltered when one of the speed components is not a number: +clean_mixture(100,u=-13,v=13,U=NaN,V=6, fast=8, slow=1)["eta"] +# set reflectivity to NaN when one of the speed components is not a number: +clean_mixture(100,u=-13,v=13,U=NaN,V=6, fast=8, slow=1, drop_missing=TRUE)["eta"] +} +\references{ +\itemize{ +\item Shi X, Drucker J, Chapman JW, Sanchez Herrera M, Dokter AM +Analysis of mixtures of birds and insects in weather radar data. +Ornithological Applications. 2025 (in press). +\item Nussbaumer R, Schmid B, Bauer S, Liechti F. +A Gaussian mixture model to separate birds and insects in +single-polarization weather radar data. Remote Sensing. 2021 May 19;13(10):1989. +} +} diff --git a/man/filter_precip.Rd b/man/filter_precip.Rd new file mode 100644 index 00000000..6d41509b --- /dev/null +++ b/man/filter_precip.Rd @@ -0,0 +1,86 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/filter_precip.R +\name{filter_precip} +\alias{filter_precip} +\title{Posthoc precipitation filter} +\usage{ +filter_precip(x, dbz = 7, range = 2500, alt_max = 3000, drop = FALSE) +} +\arguments{ +\item{x}{A \code{vp} or \code{vpts} object.} + +\item{dbz}{The minimum reflectivity factor for precipitation.} + +\item{range}{The minimum altitude range with total reflectivity factor. +\code{DBZH} > \code{dbz} at which the filter is triggered.} + +\item{alt_max}{Maximum altitude above ground level to consider.} + +\item{drop}{When \code{TRUE} the profile is removed from the} +} +\value{ +A \code{vpts} object or a \code{vp} object, depending on input \code{x}. +} +\description{ +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'). +} +\details{ +During precipitation events usually a high proportion of the altitude +column will show a high total reflectivity \code{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 \code{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 \code{DBZH} +and calculates the altitude range at which \code{DBZH} is larger than parameter +\code{dbz}. Whenever this altitude range is larger than parameter \code{range} (and \code{drop} +is \code{FALSE}), the biology is removed by setting profile quantities \code{dens} and \code{eta} +to zero and profile quantity \code{dbz} to \code{-Inf}. When parameter \code{drop} +is \code{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 \code{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 (\code{dbz}) for precipitation +below 7 dBZ is not recommended, as most precipitation has a reflectivity above 7 dBZ. + +Parameter \code{range} should be chosen carefully, and should be higher than the +typical altitude where biological scatterers still reach a reflectivity factor equal to \code{dbz}. + +Note that at S-band wavelengths bird migration occurs much more frequently in the reflectivity +regime for precipitation than at C-band. Therefore, at C-band lower settings for parameter \code{dbz} +are appropriate than at S-band. +} +\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) +} +\seealso{ +\itemize{ +\item \link{eta_to_dbz} +\item \link{dbz_to_eta} +} +} diff --git a/tests/testthat/test-as.vp.R b/tests/testthat/test-as.vp.R new file mode 100644 index 00000000..1e3297d4 --- /dev/null +++ b/tests/testthat/test-as.vp.R @@ -0,0 +1,6 @@ +test_that("as.vp() returns valid data", { + df <- as.data.frame(example_vp, suntime=FALSE) + expect_s3_class(as.vp(df), "vp") + df <- as.data.frame(example_vpts, suntime=FALSE) + expect_error(as.vp(df),"not a single vertical profile") +}) diff --git a/tests/testthat/test-clean_mixture.R b/tests/testthat/test-clean_mixture.R new file mode 100644 index 00000000..af07206d --- /dev/null +++ b/tests/testthat/test-clean_mixture.R @@ -0,0 +1,37 @@ +# create vp and vpts object and add dummy wind values +vp <- example_vp +vp$data$u_wind <-seq(-1,-8, length.out=dim(vp)[1]) +vp$data$v_wind <-seq(0,-3, length.out=dim(vp)[1]) +clean_mixture(vp, fast=12) + +vpts <- example_vpts +vpts$data$u_wind <- matrix(data=seq(0,10, length.out=dim(vpts)[2]),nrow=dim(vpts)[2],ncol=dim(vpts)[1]) +vpts$data$v_wind <- matrix(data=seq(-2,6, length.out=dim(vpts)[2]),nrow=dim(vpts)[2],ncol=dim(vpts)[1]) +clean_mixture(vpts) + +test_that("clean_mixture() requires correct input data and parameters", { + expect_error(clean_mixture(example_vp), "requires paired wind data") + expect_error(clean_mixture(example_vpts), "requires paired wind data") + expect_error(clean_mixture(vpts, slow=10, fast=1), "fast not greater than slow") +}) + +test_that("clean_mixture() returns vp or vpts object", { + expect_s3_class(clean_mixture(vpts), "vpts") + expect_s3_class(clean_mixture(vp), "vp") +}) + +test_that("clean_mixture() returns warning for infinite wind data", { + vp$data$u_wind[3] <- NA + expect_warning(clean_mixture(vp), "wind contains non finite values") +}) + +test_that("clean_mixture() returns expected values", { + # only birds + expect_equal(clean_mixture(100,u=10,v=0,U=5,V=0, fast=5, slow=1)$f,0) + # only insects + expect_equal(clean_mixture(100,u=6,v=0,U=5,V=0, fast=5, slow=1)$f,1) + # observed airspeed faster than fast component + expect_equal(clean_mixture(100,u=11,v=0,U=5,V=0, fast=5, slow=1)$f,0) + # observed airspeed slower than slow component + expect_equal(clean_mixture(100,u=5.5,v=0,U=5,V=0, fast=5, slow=1)$f,1) +})