Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

add height quantile option to integrate_profile() #485

Merged
merged 4 commits into from
Dec 16, 2021
Merged

add height quantile option to integrate_profile() #485

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
7417358
add height quantile option for vpts
adokter Dec 14, 2021
01f87fb
add height_quantile option to integrate_profile()
adokter Dec 14, 2021
f0cf26e
revert snaps
adokter Dec 14, 2021
e8a4c66
add example for using new height_quantile argument
adokter Dec 14, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
89 changes: 80 additions & 9 deletions R/integrate_profile.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,9 @@
#' seconds. Traffic rates are set to zero at times \code{t} for which no
#' profiles can be found within the period \code{t-interval_max/2} to
#' \code{t+interval_max/2}. Ignored for single profiles of class \code{vp}.
#' @param height_quantile For default `NA` the calculated height equals
#' the mean flight altitude. Otherwise a number between 0 and 1 specifying a
#' quantile of the height distribution.
#'
#' @return an object of class \code{vpi}, a data frame with vertically
#' integrated profile quantities
Expand Down Expand Up @@ -160,21 +163,32 @@
#' plot(integrate_profile(example_vpts, alt_min = 1000))
#' # plot the (cumulative) migration traffic
#' plot(integrate_profile(example_vpts), quantity = "mt")
#' # calculate median flight altitude (instead of default mean)
#' integrate_profile(example_vp, height_quantile=.5)
#' # calculate the 90% percentile of the flight altitude distribution
#' integrate_profile(example_vpts, height_quantile=.9)
integrate_profile <- function(x, alt_min, alt_max,
alpha = NA, interval_max = Inf) {
alpha = NA, interval_max = Inf,
height_quantile = NA) {
UseMethod("integrate_profile", x)
}

#' @describeIn integrate_profile Vertically integrate a vertical profile.
#'
#' @export
integrate_profile.vp <- function(x, alt_min = 0, alt_max = Inf, alpha = NA,
interval_max = Inf) {
interval_max = Inf, height_quantile = NA) {
stopifnot(inherits(x, "vp"))
stopifnot(is.numeric(alt_min) | alt_min=="antenna")
stopifnot(is.numeric(alt_max))
stopifnot(is.na(alpha) || is.numeric(alpha))

assert_that(is.scalar(height_quantile))
if(!is.na(height_quantile)){
assert_that(is.number(height_quantile))
assert_that(height_quantile>0 && height_quantile<1)
}

if (alt_min=="antenna"){
alt_min = x$attributes$where$height
}
Expand Down Expand Up @@ -217,7 +231,28 @@ integrate_profile.vp <- function(x, alt_min = 0, alt_max = Inf, alpha = NA,
weight_densdh[is.na(weight_densdh)] <- 0
weight_densdh <- weight_densdh / sum(weight_densdh)

height <- weighted.mean(get_quantity(x, "height") + interval / 2, weight_densdh, na.rm = TRUE)
if(is.na(height_quantile)){
# default (no height_quantile specified) is calculating the mean altitude
height <- weighted.mean(get_quantity(x, "height") + interval / 2, weight_densdh, na.rm = TRUE)
}
else{
# calculate a quantile of the flight altitude distribution
# 1) integrate over altitude
denscum=cumsum(weight_densdh)
denscum[is.na(denscum)]=0
# 2) find lowerbound index:
height_index_lower=findInterval(height_quantile, denscum)
# 3) find the two height bins closest to the quantile of interest
height_lower=x$data$height[height_index_lower] + interval / 2
height_upper=x$data$height[min(height_index_lower+1,length(denscum))] + interval / 2
height_quantile_lower <- denscum[height_index_lower]
height_quantile_upper <- denscum[min(height_index_lower+1,length(denscum))]
# 4) do a linear interpolation to estimate the altitude at the quantile of interest
delta_linear_interpolation <- (height_quantile-height_quantile_lower)*(height_upper-height_lower)/(height_quantile_upper-height_quantile_lower)
if(is.na(delta_linear_interpolation)) delta_linear_interpolation=0
# 5) store the quantile flight altitude as height
height <- height_lower+delta_linear_interpolation
}

u <- weighted.mean(get_quantity(x, "u"), weight_densdh, na.rm = TRUE)
v <- weighted.mean(get_quantity(x, "v"), weight_densdh, na.rm = TRUE)
Expand All @@ -240,16 +275,18 @@ integrate_profile.vp <- function(x, alt_min = 0, alt_max = Inf, alpha = NA,
output$heading <- weighted.mean((pi / 2 - atan2(airspeed_v, airspeed_u)) * 180 / pi, weight_densdh, na.rm = TRUE)
output$airspeed_u <- weighted.mean(airspeed_u, weight_densdh, na.rm = TRUE)
output$airspeed_v <- weighted.mean(airspeed_u, weight_densdh, na.rm = TRUE)
output$ff_wind <- weighted.mean(sqrt(u_wind^2 + v_wind^2), weight_densdh, na.rm = TRUE)
output$u_wind <- weighted.mean(u_wind, weight_densdh, na.rm = TRUE)
output$v_wind <- weighted.mean(v_wind, weight_densdh, na.rm = TRUE)
output$ff_wind <- weighted.mean(sqrt(get_quantity(x,"u_wind")^2 + get_quantity(x,"v_wind")^2), weight_densdh, na.rm = TRUE)
output$u_wind <- weighted.mean(get_quantity(x,"u_wind"), weight_densdh, na.rm = TRUE)
output$v_wind <- weighted.mean(get_quantity(x,"v_wind"), weight_densdh, na.rm = TRUE)
}

class(output) <- c("vpi", "data.frame")
rownames(output) <- NULL
attributes(output)$alt_min <- alt_min
attributes(output)$alt_max <- alt_max
attributes(output)$alpha <- alpha
attributes(output)$interval_max <- interval_max
attributes(output)$height_quantile <- height_quantile
attributes(output)$rcs <- rcs(x)
attributes(output)$lat <- x$attributes$where$lat
attributes(output)$lon <- x$attributes$where$lon
Expand All @@ -261,7 +298,8 @@ integrate_profile.vp <- function(x, alt_min = 0, alt_max = Inf, alpha = NA,
#'
#' @export
integrate_profile.list <- function(x, alt_min = 0, alt_max = Inf,
alpha = NA, interval_max = Inf) {
alpha = NA, interval_max = Inf,
height_quantile = NA) {
vptest <- sapply(x, function(y) is(y, "vp"))
if (FALSE %in% vptest) {
stop("requires list of vp objects as input")
Expand All @@ -279,6 +317,8 @@ integrate_profile.list <- function(x, alt_min = 0, alt_max = Inf,
attributes(output)$alt_min <- alt_min
attributes(output)$alt_max <- alt_max
attributes(output)$alpha <- alpha
attributes(output)$interval_max <- interval_max
attributes(output)$height_quantile <- height_quantile
attributes(output)$rcs <- rcs(x)
# TODO set lat/lon attributes
return(output)
Expand All @@ -289,12 +329,19 @@ integrate_profile.list <- function(x, alt_min = 0, alt_max = Inf,
#'
#' @export
integrate_profile.vpts <- function(x, alt_min = 0, alt_max = Inf,
alpha = NA, interval_max = Inf) {
alpha = NA, interval_max = Inf,
height_quantile = NA) {
stopifnot(inherits(x, "vpts"))
stopifnot(is.numeric(alt_min) | alt_min=="antenna")
stopifnot(is.numeric(alt_max))
stopifnot(is.na(alpha) || is.numeric(alpha))

assert_that(is.scalar(height_quantile))
if(!is.na(height_quantile)){
assert_that(is.number(height_quantile))
assert_that(height_quantile>0 && height_quantile<1)
}

# Integrate from antenna height
if (alt_min=="antenna"){
alt_min = x$attributes$where$height
Expand Down Expand Up @@ -341,7 +388,29 @@ integrate_profile.vpts <- function(x, alt_min = 0, alt_max = Inf,
# Find index where no bird are present
no_bird <- is.na(colSums(weight_densdh))

height <- colSums( (get_quantity(x, "height") + interval / 2) * weight_densdh, na.rm = T)
if(is.na(height_quantile)){
# default (no height_quantile specified) is calculating the mean altitude
height <- colSums( (get_quantity(x, "height") + interval / 2) * weight_densdh, na.rm = T)
}
else{
# calculate a quantile of the flight altitude distribution
# 1) integrate over altitude
denscum=apply(weight_densdh, 2, cumsum)
denscum[is.na(denscum)]=0
# 2) find lowerbound index:
height_index_lower=apply(denscum,2,findInterval,x=height_quantile)
# 3) find the two height bins closest to the quantile of interest
height_lower=x$height[height_index_lower] + interval / 2
height_upper=x$height[pmin(height_index_lower+1,nrow(denscum))] + interval / 2
height_quantile_lower <- denscum[seq(0,nrow(denscum)*(ncol(denscum)-1),nrow(denscum))+height_index_lower]
height_quantile_upper <- denscum[seq(0,nrow(denscum)*(ncol(denscum)-1),nrow(denscum))+pmin(height_index_lower+1,nrow(denscum))]
# 4) do a linear interpolation to estimate the altitude at the quantile of interest
delta_linear_interpolation <- (height_quantile-height_quantile_lower)*(height_upper-height_lower)/(height_quantile_upper-height_quantile_lower)
delta_linear_interpolation[is.na(delta_linear_interpolation)]=0
# 5) store the quantile flight altitude as height
height <- height_lower+delta_linear_interpolation
}

height[no_bird] <- NA
u <- colSums( get_quantity(x, "u") * weight_densdh, na.rm = T)
u[no_bird] <- NA
Expand Down Expand Up @@ -384,6 +453,8 @@ integrate_profile.vpts <- function(x, alt_min = 0, alt_max = Inf,
attributes(output)$alt_min <- alt_min
attributes(output)$alt_max <- alt_max
attributes(output)$alpha <- alpha
attributes(output)$interval_max <- interval_max
attributes(output)$height_quantile <- height_quantile
attributes(output)$rcs <- rcs(x)
attributes(output)$lat <- x$attributes$where$lat
attributes(output)$lon <- x$attributes$where$lon
Expand Down
2 changes: 1 addition & 1 deletion man/download_vpfiles.Rd
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

26 changes: 22 additions & 4 deletions man/integrate_profile.Rd
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 1 addition & 1 deletion man/write_pvolfile.Rd
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.