Adding velocity azimuth display plots by bart1 · Pull Request #739 · adokter/bioRad

bart1 · 2025-09-02T20:13:28Z

In this pull request i want to introduce a VAD plotting function. It helps to explain and asses the fit of velocities in a vp. There is considerable flexibility in formatting the plots. Here are some example plots:

devtools::load_all("~/bioRad")
#> ℹ Loading bioRad
#> Welcome to bioRad version 0.11.0.9000
#> 
#> using vol2birdR version 1.1.1 (MistNet not installed)
pvolfile <- system.file("extdata", "volume.h5", package = "bioRad")
example_pvol <- read_pvolfile(pvolfile)
vad(example_pvol,
  range = c(5000, 30000), height = c(200, 400),
  point_geom_args = list(size = 0.5)
)

vp <- calculate_vp(pvolfile)
#> Running vol2birdSetUp
#> Warning: radial velocities will be dealiased...
vad(example_pvol,
  vp = vp, height = c(400, 1200),
  point_geom = ggplot2::geom_bin2d,
  annotate = "sdvvp: {round(sd_vvp,2)} [m/s]",
) + ggplot2::scale_fill_viridis_c()

adokter · 2025-09-04T15:01:10Z

Thanks @bart1, this will be helpful.

…_geom`

bart1 · 2025-09-15T12:28:30Z

I have updated the documentation with suggestions from @CeciliaNilsson709, thanks!

adokter · 2025-09-17T01:25:28Z

+          max_bin_height = hgt + vp$attributes$where$interval,
+          height_bin = bin
+        )),
+        fun = function(x, v, a) cos((x - a) / 180 * pi) * v,


There is an issue with wanting to plot data from multiple elevation scans using the same cosine function as done here. Because the radar measures radial velocity, there should be cos(elevation_angle) factor in this plotting function. E.g. for the extreme case where the radar looks vertically up in a bird bath scan, you will have radial velocities all zero for any speed and direction. Therefore you cannot just plot the same curve to data from different elevation scans.

It seems more natural for this function to not act on a polar volume, but on a scan - in general it would be natural to have a method for a scan object as well.

If you really want a method that plots all data for a volume, you might be able to "correct" the raw radial velocity data for the reduction in radial velocity due to the elevation angle. Then you might be able to plot multiple scans in one graph - but we would have to very carefully document what is actually happening and that the user is no longer looking at pristine raw data, and I see some complications in doing such a correction for velocity data that is folded. For a 20 degree elevation angle scan you are looking at a correction of more than 5%. You probably want some upper/lower limits on elevation angle as well, for example to exclude extremes like the bird bath scan, or other very steep elevation angles that need a lot of correction.

This is the function we fit the radial velocity field to, with theta elevation angle, phi azimuth, and u,v the components of the horizontal speed and w the vertical speed:
VRAD = w sin(theta) + u cos(theta) sin(phi) + v cos(theta) cos(phi)
For birds you can usually assume w is zero (but we do estimate it), so it's just a simple extra cos(theta) factor that needs to be added when plotting a scan. Or alternatively, when aligning multiple scan into one plot, the VRAD values needed to be divided by cos(theta) - but would have to give folding into the Nyquist interval more thought

Good point, I'm thinking a bit more about this. Maybe we can add a vad.scan method to include these corrections in the vp curve? That would keep the flexibility to plot everything for a quick look at a pvol and do the correct thing at a scan level. I agree that correcting the VRAD values would be quite difficult to comprehend for users due to the nyquist folding. What do you think?

adokter · 2025-09-17T01:42:31Z

+#' deviations from the sine function.
+#'
+#' @param x A polar volume from which range gates are extracted.
+#' @param ... Currently not used.


why do we need it if it is not used?

To prevent other arguments from being used unnamed. With this setup you are forced to make all other arguments explicit. I generally like this setup as it avoids error with the calling order of arguments

ok - we don't use this in other functions. Could you add an issue to discuss if we should use this consistently throughout the package, so others can weigh in as well

made an issue: #746

bart1 · 2025-09-18T10:20:53Z

@adokter I now included different ways to correct the cosine error (please check). It also warns for problems with the cosine (e.g. when large errors are likely to occur). It also uses a corrected vp curve for a single scan.

adokter

Not a fully finished review, but these are the most important points, let me know what you think

adokter · 2025-09-23T02:53:48Z

+#' deviations from the sine function.
+#'
+#' @param x A polar volume from which range gates are extracted.
+#' @param ... Currently not used.


ok - we don't use this in other functions. Could you add an issue to discuss if we should use this consistently throughout the package, so others can weigh in as well

adokter · 2025-09-23T02:55:48Z

+#' Use `NULL` if no annotation is desired.
+#' @param annotation_color The color used for the `vp` annotations and line.
+#' @param annotation_size The text size used for the annotation.
+#' @param cosine_correction A character option to select what approach should be taken to correct for the fact


I think the word cosine_correction is misleading - there is no ture correction, the cosine is simply part of the model fit, so I think we shouldn't confuse the reader with there being a correction applied. We should just plot the fit and the datapoints, nothing more.

The main sticking point remains how to plot multiple sweeps into one plot. You either need to adjust the positions of the points to "absorb" the part of the model variation with elevation angle (what you call "range_gates"), so you can keep one fitting line, or - alternatively - plot the points as they are but include multiple fit lines for each sweep (coloring points and fitting lines by sweep number).

To me the more natural user option would be to be able to toggle between those two. The user argument would then by align_scans = TRUE/FALSE, with the FALSE case plotting multiple fit lines (not requiring any adjustments of the points), and the TRUE case adjusting the points and showing only one line.

your probably right the current name is unfortunate. However I would prefer to also keep a version that does not do changes (for example as the corrections might be incorrect for nyquist folded data). So I need to think about the best way to do that with an align_scans argument.

adokter · 2025-09-23T02:58:36Z

+vad.pvol <- function(x, vp = NULL, ...,
+                     range_min = NULL, range_max = NULL,
+                     alt_min = NULL, alt_max = NULL,
+                     range_gate_filter =


I think this argument is overly complicated. Why not give single arguments rho_hv (same name used in calculate_vp() and eta_max to set the two thresholds. I realize this gives you less flexibility, but I think it's better to keep things simple for the user.

I think you then you would still need several extra arguments like reflectivity_quantity and rho_quantity. Maybe a personal preference but I tend to want to keep the number of arguments limited. Using explicit thresholds would also reduce the option to use other filters (e.g. if you have a bird classifier or use differential reflectivity). I think for a function mostly useful for exploration this flexibility is valuable

I see your point, but it's very hard to understand for users what this does and how to modify it - and therefore prone to errors if people do try to change this. Could this be included as a customizable filter function to be included in the package, for which you can change the reflectivity and rho quantities and thresholds? There should be a way to simplify this

dplyr::if_all(utils::head(dplyr::matches(c("^DBZ$", "^DBZH$", "^DBZV$", "^TH$", "^TV$")), 1), \(dbz) dbz_to_eta(dbz, !!x$attributes$how$wavelength) < 36000 ) & dplyr::if_any(dplyr::matches("RHOHV"), \(rhohv) rhohv < .95)

adokter · 2025-09-23T03:00:39Z

+    if (length(elangles) != 1 && cosine_correction == "vp") {
+      warning("The data to plot is based on multiple elevation angles. To correct the `vp` data one elevation angle is needed, it is therefore averaged resulting in a inperfect correction.")
+    }
+    mean_elangle_cos <- cospi(mean(elangles) / 180) # TODO should this be a weighted mean?


I think instead of the mean, points should be corrected according to the elevation scan they belong to, so different cosine corrections for each of the elevation scans of the volume

That is done when you use cosine_correction="range_gates", then the points/range_gates are corrected. Here I adjust the sine curve. What could be improved here I now see is take the mean elevation angle per panel/height bin that we plot. Alternatively you would have to include a sine for each elevation angle that occurs in a height bin but that gets messy. For those small difference I would prefer one sine curve. Alternatively you use cosine_correction="range_gates" to avoid this averaging

Even if it's a little messy, I feel its more accurate than a mean, which we would have to explain carefully. You could try different colors for the points and sine curves of different sweeps, maybe that looks cleaner?

adokter · 2025-09-23T03:13:16Z

+                     annotate = "{round(ff,1)} m/s, {round(dd)}\u00B0",
+                     annotation_size = 4,
+                     annotation_color = "red",
+                     cosine_correction = c("none", "vp", "range_gates")) {


To keep it simple for the user, why don't we use what you now have under "vp" when plotting a single scan (so taking account the cosine in the plotted sine curve), and when you plot a "pvol" what you now have under "range_gates". Alternatively, you could plot multiple sine curves in different colors, and highlight points in the same colors to indicate which scan they belong too.

adokter · 2025-09-23T13:45:49Z

+#'
+#' As for the radial velocity to plot the first of `VRAD`, `VRADH` or `VRADV` is used.
+#'
+#' For these plots no vertical velocity is assumed as aeroecology will predominantly move in the horizontal plane.


I think we should take into account the fitted vertical velocity when plotting the model fit from the vp, it is part of the fit and easy to incorporate - available as quantity w in the profile. The full radial velocity functional form is

VRAD = w sin(theta) + u cos(theta) sin(phi) + v cos(theta) cos(phi)

so you only need to add a term w sin(theta)

Not having looked into the vertical velocity estimates very much I do wonder if this would not add more noise then it resolves. I see the estimates are considered unreliable (?summary.vp) and sometimes are quite high (e.g. -15 m/s in example_vp). On the other hand even these high speeds will not change the sin very much for moderate elevation angles. I will look into this

Right, because the angles are so low high vertical velocities have very little effect on the fit. I just think it would be cleaner to show the actual fit instead of something else. Or have an option to remove vertical velocity in the plot

…outliers

adokter · 2025-09-23T14:30:52Z

+#' @param ... Currently not used.
+#' @param range_min,range_max Numeric. The minimum and maximum range to include, in m. Values are taken from the `vp`
+#' object if provided.
+#' @param alt_min,alt_max Numeric. The minimum and maximum altitude to include, in m. If only `alt_min` value is provided


the function errors with uninformative error message when selecting a layer outside the altitude range where profile data is available, e.g. alt_min=8000, alt_max=8500

Good point, what to do when no range gate data gets selected. We could error however that can be annoying when using the function in a loop the alternative is to return an empty plot

agree - throwing a warning and returning an empty plot sounds good

adokter · 2025-09-23T14:34:52Z

+#' @param range_min,range_max Numeric. The minimum and maximum range to include, in m. Values are taken from the `vp`
+#' object if provided.
+#' @param alt_min,alt_max Numeric. The minimum and maximum altitude to include, in m. Separate panels will be plot for each
+#' altitude layer of the profile `vp`. If only `alt_min` is provided


when providing a single alt_min value, the plot doesn't indicate the altitude bin, as it does when you plot multiple altitude bins. better to be consistent and always show the altitude range as a label on the top.

Related: when I specify a narrow altitude range, like c(700,810), the panels are labeled as 600-800 and 800-1000, while it should be 700-800 and 800-810

Related: when I specify a narrow altitude range, like c(700,810), the panels are labeled as 600-800 and 800-1000, while it should be 700-800 and 800-810

As I'm in a workshop I don't have much time this week. A quick remark here is that the labeling in this case would not reflect the altitude range of the vp sin that is plotted, as is currently the case. I guess either way is imperfect but up to now I have been thinking about the altitude bins as those from the vp and not the height selection of the pvol data

Bart added 3 commits September 2, 2025 13:42

initial draft of vad function

ef4191f

Fix documentation and enable some control for visual elements

530733b

fix cran checks

8101251

adokter added this to the 0.12.0 milestone Sep 4, 2025

adokter added the new feature label Sep 4, 2025

include suggestions to documentation and rename argument to `plotting…

2c9ac32

…_geom`

bart1 marked this pull request as ready for review September 15, 2025 12:28

bart1 requested a review from adokter September 15, 2025 12:28

adokter requested changes Sep 17, 2025

View reviewed changes

Bart added 10 commits September 17, 2025 08:41

remove commented code

62206b1

add comments to orient in function

bd06c4b

error in checking

52d1dba

rename to annotation_color

de356a6

towards min max alt and range

403bc2a

VRAD flexible and use eta to filter

30b606a

dbz and rhohv flexible amd document

cd3e8d4

add basic testing

7d17ef2

style test

d65eb78

include different cosine correction methods

9ee1a72

bart1 requested a review from adokter September 18, 2025 10:21

Bart and others added 5 commits September 18, 2025 12:39

typo and correction tests

f2453a8

fix nexrad bucket in tests

8f30c8d

remove httr2::req_url_query(location = "") line

287cea2

only retrieve head

0de84cf

reviewing manual entry for vad()

3250b8d

adokter requested changes Sep 23, 2025

View reviewed changes

adokter reviewed Sep 23, 2025

View reviewed changes

change to coloring by DBZH, more common use case to identify weather …

e754315

…outliers

adokter reviewed Sep 23, 2025

View reviewed changes

Comment thread R/vad.R

adokter reviewed Sep 23, 2025

View reviewed changes

clarifying description of alt_min argument

692055c

adokter reviewed Sep 23, 2025

View reviewed changes

bart1 mentioned this pull request Oct 7, 2025

Should we use ... to force making some arguments explicit #746

Open

Bart added 5 commits October 7, 2025 13:24

add velocity_quanity argument

068c721

improve documentation

6b65b66

check velocity quantity

1e4a4b6

fix check error

3d0b052

typo

bcaa1cb

Conversation

bart1 commented Sep 2, 2025

Uh oh!

adokter commented Sep 4, 2025

Uh oh!

bart1 commented Sep 15, 2025

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

adokter Sep 17, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

adokter Sep 17, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

bart1 commented Sep 18, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

adokter left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

adokter Sep 23, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

adokter Sep 23, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

bart1 Oct 7, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

bart1 Oct 7, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

adokter Sep 17, 2025 •

edited

Loading

adokter Sep 17, 2025 •

edited

Loading

bart1 commented Sep 18, 2025 •

edited

Loading

adokter Sep 23, 2025 •

edited

Loading

adokter Sep 23, 2025 •

edited

Loading

bart1 Oct 7, 2025 •

edited

Loading

bart1 Oct 7, 2025 •

edited

Loading

adokter Sep 23, 2025 •

edited

Loading