Merge pull request #176 from birdflow-science/123-bird-flux

123 bird flux

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: BirdFlowR
 Title: Predict and Visualize Bird Movement
-Version: 0.1.0.9053
+Version: 0.1.0.9054
 Authors@R: 
     c(person("Ethan", "Plunkett",  email = "plunkett@umass.edu", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0003-4405-2251")),
@@ -43,6 +43,7 @@ Imports:
     stringr,
     tidyr,
     terra,
+    units,
     utils,
     viridisLite
 VignetteBuilder: knitr

NAMESPACE

@@ -15,6 +15,7 @@ S3method(st_bbox,BirdFlow)
 S3method(st_crs,BirdFlow)
 export(add_dynamic_mask)
 export(animate_distr)
+export(animate_flux)
 export(animate_movement_vectors)
 export(animate_routes)
 export(as_distr)
@@ -23,6 +24,7 @@ export(birdflow_crs)
 export(birdflow_options)
 export(build_collection_index)
 export(build_transitions)
+export(calc_flux)
 export(calc_movement_vectors)
 export(col_to_x)
 export(combine_transitions)
@@ -38,6 +40,7 @@ export(get_countries)
 export(get_dates)
 export(get_distr)
 export(get_dynamic_mask)
+export(get_marginal)
 export(get_mask)
 export(get_metadata)
 export(get_naturalearth)
@@ -74,6 +77,7 @@ export(n_timesteps)
 export(n_transitions)
 export(pad_timestep)
 export(plot_distr)
+export(plot_flux)
 export(plot_movement_vectors)
 export(plot_routes)
 export(preprocess_species)

NEWS.md

@@ -1,3 +1,41 @@
+
+# BirdFlowR 0.1.0.9054
+## Flux I
+Non-directional only, additional changes pending
+
+### Changes
+
+* The first argument of `lookup_date()` is now `x` instead of `timestep` 
+  `timestep` is added as a deprecated third argument so all old usage should 
+  work.  
+* `lookup_dates()` now works with transition names (e.g. "T_01-02") marginal 
+names ("M_02-03"), character timesteps ("T1", "T2"), and, as before, numeric 
+timesteps. Marginal and transition dates are the midpoint between the two 
+timesteps, this falls squarely between two dates and it is left up to 
+`base::mean.Date()` to resolve this, which it does by rounding down.
+
+
+* Add internal function `is_between()` to determines if a set points are on the
+great circles connecting every possible pair of active cells in a BirdFlow 
+model.  It will be used to calculate flux. Currently it is non-directional.
+
+* Add  `get_marginal()`  
+* Add `calc_flux()`, `plot_flux()`, `animate_flux()`  currently they only 
+calculate and work with non-directional flux / net movement.
+
+
+### Pending flux changes
+
+* Test `is_between()` 
+* Test `calc_flux()`
+* Test `get_marginal()`
+* Test `plot_flux()` and `animate_flux()`
+
+* Normalize?  See Dan's comment in the issue.
+* Directional `is_between()`
+* Directional `calc_flux()`
+
+
 # BirdFlowR 0.1.0.9053
 2023-03-15
 
@@ -13,6 +51,7 @@
 * `export_birdflow()` gets new arguments to control output file names. 
   Default values mimic old behavior.
 
+
 # BirdFlowR 0.1.0.9051
 2023-02-28
 

R/animate_flux.R

@@ -0,0 +1,22 @@
+#' Animate flux
+#'
+#' Animate net migration from [calc_flux()].
+#'
+#' @inheritParams plot_flux
+#' @inheritDotParams plot_flux
+#'
+#' @return A [`gganim`][gganimate::gganimate-package] object
+#' @export
+#' @inherit calc_flux examples
+#' @seealso [calc_flux()],[plot_flux()]
+animate_flux <- function(flux, bf, title = species(bf), ...) {
+  p <- plot_flux(flux, bf, ...)
+
+  anim <- p +
+    ggplot2::facet_null() +
+    gganimate::transition_manual(frames = .data$date) +
+    ggplot2::labs(title = title,
+                  subtitle = "{current_frame}")
+
+  return(anim)
+}

R/calc_flux.R

@@ -0,0 +1,136 @@
+#' Estimate net bird movement
+#'
+#' `calc_flux()` estimates the proportion of the species that passes through
+#' a set of points during each transition in a BirdFlow model.
+#'
+#' @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. Should be a data frame with `x` and  `y`
+#' columns containing point coordinates in [crs(bf)][crs()].
+#' @param radius The radius in meters around the points, used to determine if
+#' a path is passing through the point.
+#' @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.
+#' @param format The format to return the results in one of:
+#' \describe{
+#' \item{`"points"`}{Returns a list with `net_movement` a matrix or array of
+#'  movement, and `points` a data frame of either the input `points` or the
+#'  point locations }
+#' \item{`"dataframe"`}{Returns a "long" data frame with columns:
+#' * `x` and `y` coordinates of the points.
+#' * `transition` Transition code.
+#' * `movement` The total abundance moving through the point in the transition.
+#' * `date` The date associated with the transition, will be at the midpoint
+#'          between timesteps.
+#'
+#'  }
+#' \item{`"SpatRaster"`}{Returns a `terra::SpatRaster` with layers for each
+#' transition.}
+#'}
+#'
+#' @return See `format` argument.
+#' @export
+#'
+#' @examples
+#'
+#' \dontrun{
+#' bf <- BirdFlowModels::amewoo
+#' flux <- calc_flux(bf)
+#'
+#' plot_flux(flux, bf)
+#'
+#' animate_flux(flux, bf)
+#' }
+#'
+calc_flux <- function(bf, points = NULL, radius = NULL, n_directions = 1,
+                      format = NULL) {
+
+  if (n_directions != 1)
+    stop("Only one directional flux is supported at the moment.")
+
+  if (is.null(format)) {
+    if (is.null(points)) {
+      format <- "dataframe"
+    } else {
+      format <- "points"
+    }
+  }
+
+  format <- tolower(format)
+  stopifnot(format %in% c("points", "spatraster", "dataframe"))
+
+
+  result <- is_between(bf, points, radius, n_directions)
+  between <- result$between
+  points <- result$points
+
+  timesteps <- lookup_timestep_sequence(bf)
+  transitions <- lookup_transitions(bf)
+  marginals <- gsub("^T", "M", transitions)
+
+  # Empty result matrix, rows are points from between columns are timesteps
+  net_movement <-
+    matrix(NA_real_,
+           nrow = dim(between)[3],
+           ncol = n_transitions(bf),
+           dimnames = c(dimnames(between)[3],
+                        list(transition =  transitions)))
+
+
+  n_pts <- dim(between)[3]
+  for (i in seq_along(marginals)) {
+    from <- timesteps[i]
+    to <- timesteps[i + 1]
+    mar <- get_marginal(bf, marginals[i])
+    fdm <- get_dynamic_mask(bf, from)
+    tdm <- get_dynamic_mask(bf, to)
+
+    # subset betweenness (sb) to conform to this marginals dynamic masks
+    sb <- between[fdm, tdm, ]
+    stopifnot(all(dim(sb)[1:2] == dim(mar)))  # verify
+
+    for (j in seq_len(n_pts)){
+      net_movement[j, i] <- sum(mar[sb[, , j]])
+    }
+  }
+
+
+  if (format == "points") {
+    return(list(net_movement, points))
+  }
+
+  if (format == "spatraster") {
+    raster <- array(data = NA, dim = c(nrow(bf), ncol(bf), ncol(net_movement)))
+    dimnames(raster) <- list(row = NULL, col = NULL,
+                            transition = colnames(net_movement))
+    rc <- cbind(y_to_row(points$y, bf), x_to_col(points$x, bf))
+    colnames(rc) <- c("row", "col")
+    for (i in seq_along(marginals)) {
+      # The line below uses matrix indexing of an array where each
+      # row of the matrix defines a particular cell by treating the
+      # values in that row as an index on the dimensions of the array
+      # eg 1, 10, 40 will index rast[1, 10, 30]
+      rcl <- cbind(rc, loc = i)
+      raster[rcl] <- net_movement[, i]
+    }
+    r <- terra::rast(raster, extent = bf$geom$ext, crs = bf$geom$crs)
+    names(r) <- transitions
+    return(r)
+  }
+
+  if (format == "dataframe") {
+    wide <- cbind(as.data.frame(points)[, c("x", "y")],
+                net_movement)
+    long <- tidyr::pivot_longer(wide, cols = setdiff(names(wide), c("x", "y")),
+                                names_to = "transition", values_to = "movement")
+
+    long$date <- as.character(lookup_date(long$transition, bf))
+
+    return(as.data.frame(long))
+
+  }
+
+  stop(format, "is not a recoginized format.") # shouldn't ever get here
+}

R/get_marginal.R

@@ -0,0 +1,53 @@
+#' Return a marginal matrix from a BirdFlowR model
+#'
+#'   Marginals in BirdFlow models are stored such that the cell \[i, j\]
+#'   represents the probability of the bird being in state i in the prior
+#'   timestep and state j in the next.  Thus the number of rows in the marginal
+#'   equals the number of cells within the dynamic mask for the prior timestep
+#'   and the columns count is equal to the included cells for the following
+#'   timestep.
+#'
+#' @param x A BirdFlow object
+#' @param marginal A marginal code, e.g. "M_01-02"
+#' @param from The first timestep associated with the marginal. Note marginals
+#' are always forward so the second marginal will  be `from + 1` or `1` (when
+#' `from` is the last timestep).
+#' @return A  marginal matrix
+#'
+#' @seealso [lookup_transitions()] will generate a list of the transitions
+#'   needed to predict or route between two points in time.  [get_transition()]
+#'   will return a transition matrix - often calculated on the fly from a
+#'   marginal.
+#'
+#' @export
+get_marginal <- function(x, marginal = NULL, from = NULL) {
+  if (!has_marginals(x))
+    stop("x does not have marginals.")
+
+  if (is.null(marginal)) {
+
+    stopifnot(!is.null(from),
+              is.numeric(from),
+              length(from) == 1,
+              from %in% seq_len(n_timesteps(x)))
+    if (from == n_timesteps(x)) {
+      if(!is_cyclical(x))
+        stop("x isn't cyclical so there's no marginal with from = ", from)
+      to <- 1
+    } else {
+      to <- from + 1
+    }
+
+
+    marginal <- paste0("M_",
+                      pad_timestep(from, x), "-",
+                      pad_timestep(to, x))
+  }
+
+  stopifnot(is.character(marginal), length(marginal) == 1, !is.na(marginal))
+
+  if (!marginal %in% names(x$marginals))
+    stop(marginal, "is not a marginal in x.")
+  return(x$marginals[[marginal]])
+
+}