bmtr with efficient continous detection rates

.gitignore

@@ -25,6 +25,11 @@ vignettes/*.html
 vignettes/*.pdf
 vignettes/*_files/
 
+# Rendered output from the data-raw notebooks (sources are tracked, output
+# is regenerated by knitting them)
+data-raw/*.html
+data-raw/*_files/
+
 # OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3
 .httr-oauth
 

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: BirdFlowR
 Title: Predict and Visualize Bird Movement
-Version: 0.1.0.9083
+Version: 0.1.0.9084
 Authors@R: 
     c(person("Ethan", "Plunkett",  email = "plunkett@umass.edu", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0003-4405-2251")),

NAMESPACE

@@ -122,6 +122,7 @@ export(validate_BirdFlow)
 export(validate_BirdFlowIntervals)
 export(validate_BirdFlowRoutes)
 export(validate_Routes)
+export(visualize_distance_weights)
 export(write_intervals)
 export(write_routes)
 export(x_to_col)

NEWS.md

@@ -1,4 +1,68 @@
 
+# BirdFlowR 0.1.0.9084
+2026-07-27
+
+## Continuous BMTR
+
+The focus of this update is calculating Bird Flow Migration Traffic (BMTR)
+with a detection rate that varies continuously, as opposed to the default
+method which has a binary detection rate. 
+
+Previously `calc_bmtr()` did have both methods but the continuous method 
+was prohibitively computationally expensive. This update makes it usable 
+with non-trivial models. 
+
+**Breaking change** This update changes the arguments used to call the 
+continuous version of the function dropping `weighted` and `euclidean`
+arguments.  I believe no users are affected by this change, but 
+I apologize if you were.
+
+## Specific changes
+* `calc_bmtr()`: replaced the `weighted`/`euclidean` arguments with a single
+  `method` argument (`"binary"`, `"continuous"`, or `"continuous-spherical"`)
+  that selects the detection model used to assess whether a transition
+  passes by a point. `"continuous"` uses fast planar (Euclidean) geometry
+  and is the recommended weighted option; `"continuous-spherical"` is the
+  slower great-circle version, kept to assess the impact of switching
+  geometries.
+* `calc_bmtr()`: new `...` argument, forwarded to `calc_dist_weights()` for
+  the two continuous methods, to tune the spread kernel (`kernel`, `gamma`,
+  `kl`, `s1`) used to model uncertainty in a bird's path.
+* `calc_dist_weights()`: renamed the `method` argument to `kernel` to avoid
+  colliding with `calc_bmtr()`'s new `method` argument.
+* Fixed a bug where `calc_bmtr()`'s weighted detection methods ignored a
+  user-supplied `points` argument due to a positional argument-binding bug.
+* Fixed a bug where `calc_bmtr()`'s weighted detection methods ignored a
+  user-supplied `radius` argument, always using `mean(res(bf))` instead.
+* Replaced `weight_between()` and the internal `calc_detection_rate()`
+  dispatcher with `calc_spherical_detection_rate()` and
+  `calc_euclidean_detection_rate()`.
+* `calc_dist_weights()`: `kernel` now supports 3 additional spread kernels
+  (`"m1"`, `"m5"`, `"sq"`) alongside the existing `"m3"` and `"bb"`.
+* New `visualize_distance_weights()`: plots the shape of `calc_dist_weights()`'s
+  spread kernel, as an SD envelope or a weight raster, to help choose
+  `kernel` and its hyperparameters before running `calc_bmtr()`.
+* Fixed a units bug in `calc_dist_weights()`'s `"m1"`, `"m3"`, `"m5"`, and
+  `"sq"` kernels: `dist_along_line`/`line_lengths`/`res_m` were converted to
+  km while `gamma`/`kl` were left in meters, making the kernel's contribution
+  to the spread negligible next to the fixed nugget at any realistic line
+  length. Everything is now kept in meters throughout, matching `gamma`/`kl`.
+  This changes the numeric output of `calc_bmtr()`'s `"continuous"`/
+  `"continuous-spherical"` methods and of `calc_euclidean_detection_rate()`/
+  `calc_spherical_detection_rate()` for all kernels, including `"bb"` (whose
+  `t`/`len` inputs are now meters rather than km, so its `s1` hyperparameter
+  means something different than before).
+* `calc_dist_weights()`/`visualize_distance_weights()`: `gamma`/`kl`/`s1`
+  now default to `3e10`/`9e5`/`200` instead of the old fixed literals
+  (`40000`/`2000`/`10`) tuned for a ~1-2 km demo scale. The three
+  hyperparameters have different units (`gamma` is a variance, `kl` a
+  length, `s1` a sqrt-length). The new values were chosen by eye against
+  `visualize_distance_weights()`'s envelope plots. 
+  See `data-raw/kernel_parameter_tuning.Rmd`. 
+  This changes the default numeric output of `calc_bmtr()`'s `"continuous"`/
+  `"continuous-spherical"` and of
+  `calc_euclidean_detection_rate()`/`calc_spherical_detection_rate()`.
+
 # BirdFlowR 0.1.0.9083
 2026-06-15
 

R/calc_bmtr.R

@@ -18,10 +18,14 @@
 #'
 #' @section Limitations:
 #'
-#' `calc_bmtr()` makes the incorrect simplifying assumption
-#'  that birds follow the shortest (great circle) path
-#' between the center of the the source and destination raster cells.  Caution
-#' should be used when interpreting the results especially around
+#' `calc_bmtr()` makes an incorrect simplifying assumption about the path
+#' birds take: either that they follow the shortest path between the centers
+#' of the source and destination raster cells (`method = "binary"`), or that
+#' they are distributed symmetrically in a Gaussian distribution around that
+#' path (the two continuous methods). That path is a great circle in all
+#' cases except `method = "continuous"`, which uses a straight line in the
+#' model's projected CRS.
+#' Caution should be used when interpreting the results especially around
 #' major geographic features such as coasts, large lakes, mountain ranges, and
 #' ecological system boundaries that might result in non-linear migration paths.
 #'
@@ -41,16 +45,19 @@
 #' standardizing the units based on the entire cell area that that point
 #' represents.
 #'
-#'
 #' @param bf A BirdFlow model
 #' @param points A set of points to calculate movement through. If `points` is
 #' `NULL` they will default to the BirdFlow model cells that are either active
 #' or fall between two active cells. Otherwise a data frame with `x` and  `y`
 #' columns containing point coordinates in [crs(bf)][terra::crs()].
-#' @param radius The radius in meters around the points used to assess whether
-#' a movement line passes by (or through) the point. If a point is farther than
-#' `radius` from a great circle line between two cells centers then it is not
-#' between them.
+#' @param radius The radius in meters around the points used to assess the
+#' detection rate for a movement at the point.
+#' With `method = "binary"`, if a point is within `radius` of the great
+#' circle line between two cell centers then the movement is detected at
+#' that point. For the two continuous detection methods there is a
+#' probability distribution for the location of the bird as it passes by the
+#' point, and the radius defines the band over which that probability is
+#' integrated, giving the detection rate.
 #' @param n_directions The number of directional bins to use for recording
 #' movement direction. Must be either `1` indicating no direction information
 #' or an even number. This is a placeholder, currently only `1` is supported.
@@ -71,32 +78,71 @@
 #' transition.}
 #'}
 #' @inheritParams is_between
-#' @param weighted If `FALSE` use the original and quicker version of bmtr
-#' that sums all the marginal probability for transitions that pass within a
-#' fixed distance of the point.  If `TRUE` assign a weight to the point and
-#' transition combo that then is multiplied by the marginal probability before
-#' summing.  This argument is experimental but the default value is identical
-#' to the old version. The argument name and behavior when set to `TRUE` may
-#' change.
+#' @param method The detection model used to determine how much of a
+#' transition's movement counts towards a point's BMTR:
+#' \describe{
+#' \item{`"binary"`}{(default) Fast and deterministic. A movement line either
+#' does or does not pass within `radius` of the point, per [is_between()].}
+#' \item{`"continuous"`}{Assigns a continuous weight (0 to 1) based on the
+#' probability that a bird's actual path, modeled as spreading away from the
+#' straight line between two cells, passes within `radius` of the point. Uses
+#' planar (Euclidean) geometry in the model's native CRS, and is the
+#' recommended detection model when continuous weighting is desired.}
+#' \item{`"continuous-spherical"`}{The same continuous weighting as
+#' `"continuous"`, but computed with great-circle (spherical) geometry
+#' instead. Much slower, and not recommended for routine use; kept to allow
+#' assessing the impact of switching from spherical to Euclidean geometry.}
+#' }
+#' @param ... For `method = "continuous"` or `"continuous-spherical"`,
+#' additional arguments forwarded to [calc_dist_weights()] to control the
+#' spread kernel used to model uncertainty in a bird's path: `kernel`,
+#' `gamma`, `kl`, and `s1`. See [calc_dist_weights()] for the full list of
+#' supported kernels and their hyperparameters, and
+#' [visualize_distance_weights()] to explore how they shape the spread
+#' before running this (potentially expensive) function. Not applicable,
+#' and an error, for `method = "binary"`.
 #'
 #' @return See `format` argument.
+#' @seealso [visualize_distance_weights()] to explore the spread kernel
+#' hyperparameters accepted via `...`.
 #' @export
 #'
 #' @examples
 #'
 #' \dontrun{
 #' bf <- BirdFlowModels::amewoo
-#' bmtr <- calc_bmtr(bf)
 #'
+#' # Binary detection along shortest path
+#' bmtr <- calc_bmtr(bf)
 #' plot_bmtr(bmtr, bf)
 #'
 #' animate_bmtr(bmtr, bf)
+#'
+#' # Continuous detection with a wider spread kernel
+#' # (gamma defaults to 3e10)
+#' bmtr2 <- calc_bmtr(bf, method = "continuous", gamma = 6e10)
+#' animate_bmtr(bmtr2, bf)
+#'
+#' # Visualize spread for the continuous detection with wider spread
+#' visualize_distance_weights(line_lengths = c(2, 5, 10) * xres(bf),
+#'                            gamma = 6e10, res_m = xres(bf))
+#'
+#' # Visualize spread for the continuous detection with default values
+#' visualize_distance_weights(line_lengths = c(2, 5, 10) * xres(bf),
+#'                            res_m = xres(bf))
+#'
 #' }
 #'
+#'
+#'
+#'
 calc_bmtr <- function(bf, points = NULL, radius = NULL, n_directions = 1,
                       format = NULL, batch_size = 5e5, check_radius = TRUE,
-                      weighted = FALSE) {
+                      method = c("binary", "continuous",
+                                "continuous-spherical"),
+                      ...) {
 
+  method <- match.arg(method)
 
   if (!requireNamespace("SparseArray", quietly = TRUE)) {
     stop("The SparseArray package is required to use calc_bmtr(). ",
@@ -117,15 +163,23 @@ calc_bmtr <- function(bf, points = NULL, radius = NULL, n_directions = 1,
   format <- tolower(format)
   stopifnot(format %in% c("points", "spatraster", "dataframe"))
 
-  # The only difference between is_between() and weight_between() return formats
-  # is in the "between" component.  In the first it's logical (TRUE is between)
-  # in the second it is a weight between 0 and 1 (non zero indicates some
-  # level of betweeness".
-  if (weighted) {
-    result <- weight_between(bf, points, radius, n_directions)
-  } else {
-    result <- is_between(bf, points, radius, n_directions)
-  }
+  # The only difference between is_between() and the continuous detection
+  # functions' return formats is in the "between" component.  In the first
+  # it's logical (TRUE is between) in the second it is a weight between 0
+  # and 1 (non zero indicates some level of betweenness).
+  result <- switch(
+    method,
+    binary = is_between(bf, points = points, radius = radius,
+                        n_directions = n_directions,
+                        batch_size = batch_size, check_radius = check_radius,
+                        ...),
+    continuous = calc_euclidean_detection_rate(
+      bf, points = points, radius = radius, n_directions = n_directions,
+      batch_size = batch_size, check_radius = check_radius, ...),
+    "continuous-spherical" = calc_spherical_detection_rate(
+      bf, points = points, radius = radius, n_directions = n_directions,
+      batch_size = batch_size, check_radius = check_radius, ...)
+  )
 
   between <- result$between
   points <- result$points

R/calc_distance_weights.R

@@ -24,62 +24,92 @@
 #' determine the nugget added to the variance to represent the uncertainty in
 #' the starting and ending location of the transition.
 #'
-#' @param method  The method used for calculating the standard deviation
-#' in the probability distribution.  Currently `"m3"`, Martern 3/2; and
-#' `"bb"`, Brownian bridge are supported.
+#' @param kernel The kernel used for calculating the standard deviation
+#' in the probability distribution:
+#' \describe{
+#' \item{`"m3"`}{(default) Matern 3/2. Driven by `gamma` and `kl`.}
+#' \item{`"m1"`}{Matern 1/2. Driven by `gamma` and `kl`.}
+#' \item{`"m5"`}{Matern 5/2. Driven by `gamma` and `kl`.}
+#' \item{`"sq"`}{Squared-exponential (Gaussian). Driven by `gamma` and `kl`.}
+#' \item{`"bb"`}{Brownian bridge. Driven by `s1`.}
+#' }
+#' See [visualize_distance_weights()] to explore how these kernels and their
+#' hyperparameters shape the spread.
+#' @param gamma Spread magnitude hyperparameter (m^2, a variance) for the
+#' Matern-family and squared-exponential kernels (`"m1"`, `"m3"`, `"m5"`,
+#' `"sq"`). Ignored for `kernel = "bb"`. The default was tuned by eye with
+#' `visualize_distance_weights()`.
+#' @param kl Lengthscale hyperparameter (m) for the Matern-family and
+#' squared-exponential kernels (`"m1"`, `"m3"`, `"m5"`, `"sq"`). Ignored for
+#' `kernel = "bb"`. The default was tuned by eye using
+#' `visualize_distance_weights()`.
+#' @param s1 Spread magnitude hyperparameter (units of sqrt(m), not m) for
+#' the `"bb"` kernel. Ignored for all other kernels. The default was tuned by
+#' eye against `visualize_distance_weights()`.
 #'
 #' @return A vector of weights of the same length as the first three arguments.
+#' @seealso [visualize_distance_weights()]
 #' @keywords internal
 calc_dist_weights <- function(dist_to_line, dist_along_line, line_lengths,
-                              radius_m, res_m, method = "m3") {
+                              radius_m, res_m,
+                              kernel = c("m3", "bb", "m1", "m5", "sq"),
+                              gamma = 3e10, kl = 9e5, s1 = 200) {
 
+  kernel <- match.arg(kernel)
 
+  sd <- calc_dist_weights_sd(dist_along_line, line_lengths, res_m,
+                             kernel = kernel, gamma = gamma, kl = kl, s1 = s1)
 
-  # Spatial information
-  len <- line_lengths / 1000  # length of great circle in km
-  t <- dist_along_line / 1000 # Distance from start of great circle to where
-  # the point projects onto the great circle in km
-  d <- dist_to_line / 1000 # Distance from reference point to the great circle
-
-  res_km <- res_m / 1000
-
-  r <- radius_m / 1000  # Radius of transact in KM.
-
-  # s2 is the standard deviation of the nugget
-  s2 <- res_km / 4 # This initial setting of 1/4 the cell width converted
-  # to KM means the cell boundary (orthagonally) is at 2
-  #  standard deviations from the cell center, the corners
-  # will be at 2.8 SD.
-
-  valid_methods <- c("bb", # Brownian Bridge
-                     "m3") # Martern 3/2
+  # Calculate weight
+  weight <- rep(0, length(sd))
+  in_range <- (dist_to_line - radius_m) < 1.96 * sd
+  low_cum_prob <- stats::pnorm(dist_to_line[in_range] - radius_m,
+                               sd = sd[in_range])
+  high_cum_prob <- stats::pnorm(dist_to_line[in_range] + radius_m,
+                                sd = sd[in_range])
+  weight[in_range] <- high_cum_prob - low_cum_prob
 
-  if (!method %in% valid_methods)
-    stop("Method shold be one of ", paste(valid_methods, collapse = ", "))
+  return(weight)
+}
 
-  if (method == "bb") {
+#' Standard deviation of the spread kernel used by [calc_dist_weights()]
+#'
+#' Computes the standard deviation (in m) of the probability distribution
+#' representing uncertainty in a bird's location at a given point along a
+#' transition line, for one of the kernels supported by
+#' [calc_dist_weights()]. Factored out of [calc_dist_weights()] so
+#' [visualize_distance_weights()] can plot the same spread it uses without
+#' duplicating the kernel-selection logic.
+#'
+#' @param dist_along_line How far along the line the point of interest is (m)
+#' @param line_lengths How long the line is (m)
+#' @param res_m The resolution of the associated bird flow model, used to
+#' determine the nugget added to the variance to represent the uncertainty in
+#' the starting and ending location of the transition.
+#' @inheritParams calc_dist_weights
+#' @return A vector of standard deviations (m), the same length as
+#' `dist_along_line`.
+#' @keywords internal
+calc_dist_weights_sd <- function(dist_along_line, line_lengths, res_m,
+                                 kernel, gamma, kl, s1) {
 
-    # Hyperparameters
-    s1 <-  10 # Tune this visually looking at a few resulting distributions?
+  len <- line_lengths  # length of great circle (m)
+  t <- dist_along_line # Distance from start of great circle to where the
+  # point projects onto the great circle (m)
 
+  # s2 is the standard deviation of the nugget
+  s2 <- res_m / 4 # This initial setting of 1/4 the cell width means the
+  # cell boundary (orthagonally) is at 2 standard deviations from the cell
+  # center, the corners will be at 2.8 SD.
 
+  if (kernel == "bb") {
     # Standard deviation
-    sd <-  sqrt(s1^2 * t * (len - t) / len + s2^2)  # SD in KM for this point.
-
-  }
-
-  if (method == "m3") {
-    sd <- sqrt(calc_martern_variance(t, len, k_m3, gamma = 40000, kl = 2000) +
+    sd <- sqrt(s1^2 * t * (len - t) / len + s2^2)  # SD (m) for this point.
+  } else {
+    k_fun <- get(paste0("k_", kernel))
+    sd <- sqrt(calc_martern_variance(t, len, k_fun, gamma = gamma, kl = kl) +
                  s2^2)
   }
 
-
-  # Calculate weight
-  weight <- rep(0, length(sd))
-  in_range <- (d - r) < 1.96 * sd
-  low_cum_prob <- stats::pnorm(d[in_range] - r, sd = sd[in_range])
-  high_cum_prob <- stats::pnorm(d[in_range] + r, sd = sd[in_range])
-  weight[in_range] <- high_cum_prob - low_cum_prob
-
-  return(weight)
+  sd
 }