Merge pull request #173 from birdflow-science/170-extend

DESCRIPTION

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

NAMESPACE

@@ -29,13 +29,16 @@ export(combine_transitions)
 export(distribution_performance)
 export(drop_transitions)
 export(evaluate_performance)
+export(export_birdflow)
 export(export_rasters)
+export(extend_birdflow)
 export(flatten_raster)
 export(get_coastline)
 export(get_countries)
 export(get_dates)
 export(get_distr)
 export(get_dynamic_mask)
+export(get_mask)
 export(get_metadata)
 export(get_naturalearth)
 export(get_states)

NEWS.md

@@ -1,3 +1,23 @@
+# BirdFlowR 0.1.0.9051
+2023-02-28
+
+* Add `get_mask()`
+* Add `extend_birdflow()`
+* Add internal `read_geom()`
+* Add internal `extend_geom()`
+* Make `export_birdflow()` a public function.
+* Update `import_birdflow()` so that it can re-import a BirdFlow model.
+  Previously the only viable workflow was `preprocess_birdflow()`, fit model
+  in BirdFlowPy, and then`import_birdflow()` (once).  After that we wrote the 
+  model to .rda files.  Now it's possible to do the above but then save it 
+  to a new `hdf5`  with: `export_birdflow()` (to hdf5) and then import again 
+  with `import_birdflow()`.   Note, however, that the file structure of a 
+  re-exported hdf5 will not match that of the output from python, b/c on the 
+  first import R updates a bunch of data structures - for instance renaming the
+  marginals and adding a marginal index, and the rexported model retains those
+  updates.  Thus, the need to modify  `import_birdflow()` to behave 
+  appropriately with both file structures.
+
 # BirdFlowR 0.1.0.9050
 2023-02-22
 

R/export_birdflow.R

@@ -1,12 +1,26 @@
 
-
+#' @rdname export_import_birdflow
+#' @export
+#' @description
+#' `export_birdflow()` saves a BirdFlow object to an
+#' HDF5 (Hierarchical Data Format version 5)
+#' @param bf A BirdFlow model
+#' @param file The file name to write. If not supplied or if the path is to
+#' a directory the file name will default to
+#' `"<species code>_<ebirdst year>_<resolution in (km)>`.
+#' @param format The format to export either `"hdf5"` or `"rds"`.
+#' @param overwrite The default `TRUE` will delete a preexisting file and then
+#' write a new one in it's place.  if `FALSE`, `export_birdflow()` will throw
+#' an error if  `file` already exists.
+#' @return `export_birdflow()` writes a file and invisibly returns `TRUE`
+#' if successful.
 export_birdflow <- function(bf, file = NULL,
                             format = "hdf5", overwrite = TRUE) {
   format <- tolower(format) |>  match.arg(c("hdf5", "rds"))
 
   extension <- switch(format,
                       hdf5 = ".hdf5",
-                      rds = ".Rds")
+                      rds = ".rds")
 
   # If file is a directory or missing add a default file name
   # using [sp. code]_[S&T version year]_[res]km
@@ -19,10 +33,10 @@ export_birdflow <- function(bf, file = NULL,
   }
   file <- normalizePath(file, winslash = "/", mustWork = FALSE)
 
-  if (!grepl(paste0("\\", extension, "$"), file))
+  if (!grepl(paste0("\\", extension, "$"), file, ignore.case = TRUE))
     stop(paste0(
       "If file is not NULL or a directory (ending in a slash) it should end",
-         'in ".hdf5" or ".Rds" consistent with the format argument.'))
+         'in ".hdf5" or ".rds" consistent with the format argument.'))
 
   bf_msg("Exporting BirdFlow model for", species(bf), "to:\n\t", file, "\n")
 
@@ -55,4 +69,6 @@ export_birdflow <- function(bf, file = NULL,
                    write.attributes = FALSE,
                    createnewfile = i == 1)  # TRUE for first object
   }
+
+  return(invisible(TRUE))
 }

R/extend_birdflow.R

@@ -0,0 +1,162 @@
+
+
+#' Extend BirdFlow extent
+#'
+#' @param x A single BirdFlow object, or one or more paths to BirdFlow objects
+#' stored as either hdf5 or rds files.
+#' @param y An extent or an object that yields an extent when passed
+#' to [terra::ext()].
+#' @return If `x` is a BirdFlow model object `extend_birdflow()` returns an
+#' extended version of the same model.  If `x` is the path to one or more
+#'  BirdFlow models than those files are modified and a logical vector of the
+#'  same length is returned with TRUE for success.
+#' @export
+#'
+#' @examples
+#' bf <- BirdFlowModels::amewoo
+#'
+#' # Define expanded extent for example
+#' e <-  ext(bf)
+#' buffer <- 3 * res(bf)
+#' e[1] <- e[1] - buffer[1]
+#' e[2] <- e[2] + buffer[1]
+#' e[3] <- e[3] - buffer[2]
+#' e[4] <- e[4] + buffer[2]
+#'
+#' bf2 <- extend_birdflow(bf, e)
+#'
+#'# Compare initial and expandeded extents
+#'data.frame(item = names(as.vector(ext(bf))),
+#'             initial = as.vector(ext(bf)),
+#'             final = as.vector(ext(bf2)))
+#'
+#'\dontrun{
+#'# Plot both versions
+#'library(terra)
+#'plot_distr(get_distr(bf, 1), bf)
+#'plot_distr(get_distr(bf2, 1), bf2)
+#'   }
+#'
+extend_birdflow <- function(x, y) {
+
+  if (inherits(x, "BirdFlow")) {
+    x$geom <- extend_geom(x$geom, y)
+    return(x)
+  }
+
+  if (inherits(x, "character")) { # file path
+
+    # Handle vectors of file names
+    if (length(x) > 1) {
+      success <- rep(FALSE, length(x))
+      for (i in seq_along(x)) {
+        success[i] <- extend_birdflow(x[i], y)
+      }
+      return(success)
+    }
+
+    ## Single files
+
+    stopifnot(file.exists(x))
+
+    # hdf5
+    # only read, modify, and write the geometry component
+    if (grepl("\\.hdf5$", x, ignore.case = TRUE)) {
+      geom <- read_geom(x)
+      geom <- extend_geom(geom, y)
+      rhdf5::h5delete(x, name = "geom")
+      rhdf5::h5write(geom,
+                     file = x,
+                     name = "geom",
+                     native = TRUE,
+                     write.attributes = FALSE,
+                     createnewfile = FALSE)
+      return(TRUE)
+    }
+
+
+    # Rdata file - read - extend - write
+    if (grepl("\\.rds$", x, ignore.case = TRUE)) {
+      bf <- readRDS(x)
+      bf <- extend_birdflow(bf, y = y)
+      saveRDS(bf, file = x)
+      return(TRUE)
+    }
+    stop("if x is not a BirdFlow object it should be the path to an ",
+         ".hdf5 or R .rds file.")
+
+  } # end x is path
+
+  stop("x should be a BirdFlow object or a path to a file containing one.")
+}
+
+
+
+#' Extend geometry component of a BirdFlow object
+#'
+#' This is an internal helper function called twice by [extend_birdflow()]
+#' it adjust the nrow, ncol, ext, and mask elements of the
+#' geom component of a BirdFlow model to expand the extent while preserving
+#' the same number, location, and alignment of the unmasked cells -
+#' thus nothing else in the object needs to change.
+#'
+#' @param geom  The geometry component of a BirdFlow object
+#' @param y An object that returns an extent when passed to [terra::ext()],
+#' this can be an extent, a SpatRaster, or a BirdFlow model.
+#' @return extended geometry (covering larger area)
+#' @keywords internal
+extend_geom <- function(geom, y) {
+
+  eg  <- geom$ext  # extent of geom
+  ey <- ext(y)  # extent of y
+
+  is_aligned <- function(a, res, tolerance = sqrt(.Machine$double.eps)){
+    # Args:
+    # a:  extent coordinates (can be vector)
+    # res: cell size (single value)
+    #
+    # Return: A logical vector of the same length as a, indicating if
+    # each element in a is a multiple of
+    # check to see if a is a multiple of res while allowing for tiny differences
+    a <- as.numeric(as.vector(a)) # double conversion works if a is SpatExtent
+    res <- as.numeric(res)
+    m <- a %% res
+    # Ok if modulo is 0
+    test1 <-
+      sapply(m,  function(x) isTRUE(all.equal(x, 0, tolerance = tolerance)))
+    # Also OK if the modulo is the res which happens if `a` is slightly less
+    # than a multiple of res, but within the tolerance.
+    test2 <-
+      sapply(m, function(x) isTRUE(all.equal(x, res, tolerance = tolerance)))
+    return(test1 | test2)
+  }
+
+  if(!all(is_aligned(ey[c(1,2)], geom$res[1])) ||
+     !all(is_aligned(ey[c(3, 4)], geom$res[2]))) {
+    stop("The new extent (y) does not align with cells in the ",
+         "BirdFlow model (x) ", call. = FALSE)
+  }
+
+  # Check that new extent completely contains old extent
+  stopifnot(all(ey[c(1, 3)] <= eg[c(1, 3)]))
+  stopifnot(all(ey[c(2, 4)] >= eg[c(2, 4)]))
+
+  # Convert mask to SpatRaster
+  x_mask <- geom$mask
+  x_mask  <- terra::rast(x_mask, extent = geom$ext, crs = geom$crs)
+  names(x_mask) <- "mask"
+
+  # Extend the mask
+  r <- terra::extend(x = x_mask, y = ey, fill = FALSE) # raster
+
+  # Update the geom
+  geom$mask <- matrix(as.logical(terra::values(r, mat = FALSE)),
+                   nrow = terra::nrow(r),
+                   ncol = terra::ncol(r),
+                   byrow = TRUE)
+  geom$ncol <- terra::ncol(r)
+  geom$nrow <- terra::nrow(r)
+  geom$ext <- as.numeric(as.vector(ext(y)))
+
+  return(geom)
+}

R/get_mask.R

@@ -0,0 +1,81 @@
+
+#' Extract mask from BirdFlow model
+#'
+#' `get_mask()` extracts the static mask from a BirdFlow model.  The
+#' static mask is a logical raster indicating which cells are included in the
+#' model (at any timestep).  These are also the cells (in row major order)
+#' that correspond with distribution values, and location indices.
+#'
+#' @param bf A BirdFlow model
+#' @param format One of `'SpatRaster'` for a [terra::SpatRaster] object,
+#'   `'numeric'` for a matrix or array, or`'dataframe'` for raster data
+#'   suitable for plotting with [ggplot2::geom_raster()]
+#' @return
+#' The return type of `get_mask()`depends on the `format` argument:
+#' * `"SpatRaster"` (the default) returns a [terra::SpatRaster] object.
+#' * `"numeric"` returns the mask as a matrix.
+#' * `"dataframe"` will return a data frame suitable for plotting with
+#'   [ggplot2::geom_raster] with columns:
+#'   * `row`, `col` the row and column indices of each cell.
+#'   * `x`, `y` the x and y coordinates of the cell center.
+#'   * `i` the location index (in `bf`) of the cell.
+#'   * `mask`   `TRUE` for cells included in the *model*, `FALSE` for excluded
+#'   cells.
+#' @export
+#' @examples
+#' bf <- BirdFlowModels::amewoo
+#' m <- get_mask(bf)
+#'
+#'\dontrun{
+#' library(terra)
+#' plot(m)
+#'}
+#'
+get_mask <- function(bf, format = "SpatRaster") {
+
+  format <- tolower(format)
+  format <- switch(format,
+                   "spatrast" = "spatraster",
+                   "terra" = "spatraster",
+                   "matrix" = "numeric",
+                   "array" = "numeric",
+                   "raster.data.frame" = "dataframe",
+                   "data.frame" = "dataframe",
+                   format)
+
+  stopifnot("Format must be one of 'SpatRaster', 'numeric', or 'dataframe'" =
+              format %in% c("spatraster", "numeric", "dataframe"))
+
+  if (format == "spatraster") {
+    m <- bf$geom$mask
+    r <- terra::rast(m, extent = bf$geom$ext, crs = bf$geom$crs)
+    names(r) <- "mask"
+    return(r)
+  }
+
+  if (format == "numeric") {
+    return(bf$geom$mask)
+  }
+
+  if (format == "dataframe") {
+    # Data frame for use with ggplot2 geom_raster() it will
+    # have one row per cell in the full raster
+
+    # Get x, y, i, and mask columns for all cells in the full raster
+    #  i will be NA for cells that aren't in the mask
+    rows <- seq_len(nrow(bf))
+    cols <- seq_len(ncol(bf))
+    df <- expand.grid(rows, cols)
+    names(df) <- c("row", "col")
+    df$x <- col_to_x(df$col, bf)
+    df$y <- row_to_y(df$row, bf)
+    df$i <- rc_to_i(df$row, df$col, bf)
+    mv <- match(1:n_active(bf), df$i)
+    n_rast <- ncol(bf) * nrow(bf)
+    df$mask <- FALSE
+    df$mask[mv] <- TRUE
+    return(df)
+  }
+
+  stop("Unrecognized format.") # should not ever get here
+}