diff --git a/DESCRIPTION b/DESCRIPTION
index 490ad68b..7d2d1c8b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: BirdFlowR
 Title: Predict and Visualize Bird Movement
-Version: 0.1.0.9031
+Version: 0.1.0.9032
 Authors@R: 
     c(person("Ethan", "Plunkett",  email = "plunkett@umass.edu", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0003-4405-2251")),
diff --git a/NAMESPACE b/NAMESPACE
index 96082864..414d75a8 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -14,6 +14,7 @@ S3method(st_as_sf,BirdFlowRoutes)
 S3method(st_bbox,BirdFlow)
 S3method(st_crs,BirdFlow)
 export(add_dynamic_mask)
+export(animate_distr)
 export(animate_movement_vectors)
 export(animate_routes)
 export(as_distr)
@@ -65,6 +66,7 @@ export(n_distr)
 export(n_parameters)
 export(n_timesteps)
 export(n_transitions)
+export(plot_distr)
 export(plot_movement_vectors)
 export(plot_routes)
 export(preprocess_species)
@@ -99,6 +101,8 @@ importClassesFrom(Matrix,dgCMatrix)
 importClassesFrom(Matrix,dgRMatrix)
 importClassesFrom(Matrix,sparseMatrix)
 importFrom(Matrix,Matrix)
+importFrom(grDevices,gray)
+importFrom(grDevices,grey)
 importFrom(grDevices,rgb)
 importFrom(methods,as)
 importFrom(methods,is)
diff --git a/NEWS.md b/NEWS.md
index 4271a3fd..78e03483 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,15 @@
+# BirdFlowR 0.1.0.9032
+2023-09-20
+
+* BREAKING change to column names in data frame returned by `rasterize_distr(format = "dataframe")` 
+  * Old `"time"` is now `"label"` and is an ordered factor it is derived from
+  the column names in the distribution which are often but not always dates.
+  * Old `density` is now value. 
+  
+* New functions to visualize distributions:
+  * `plot_distr()` 
+  * `animate_distr()` 
+
 # BirdFlowR 0.1.0.9031
 2023-09-08
 
diff --git a/R/animate_distr.R b/R/animate_distr.R
new file mode 100644
index 00000000..5cd452f0
--- /dev/null
+++ b/R/animate_distr.R
@@ -0,0 +1,79 @@
+#' Animate distributions
+#'
+#' Animate distributions as produced by `get_distr()`,
+#' `as_distr()`, or `predict()`.  The distributions will be displayed
+#' in the column order in `distr` and column labels will be used as
+#' plot subtitle.
+#'
+#' @param distr A set of distributions it should be a matrix with `n_active(bf)`
+#' rowsand a column for each distribution. The animation will proceed in the
+#' in column order and column names will be used a subtitles in the plot.
+#' @param bf The BirdFlow object that `distr` is associated with.
+#' @param title The title to use for the animation. The default is
+#' the common name of the species.
+#' @inheritDotParams plot_distr -distr -bf -title
+#'
+#' @return A gganimate object that can be displayed with `print()` or
+#'  or `gganimate::animate()`.  See example for how to export to a file.
+#' @export
+#'
+#' @examples
+#'
+#' # Animate distributions from BirdFlow object - derived from
+#' # eBird Status and Trends:
+#'
+#' bf <- BirdFlowModels::amewoo
+#' ts <- lookup_timestep_sequence(bf, season = "prebreeding")
+#' distr <- get_distr(bf, ts)
+#' anim <- animate_distr(distr, bf,  show_dynamic_mask = TRUE)
+#'
+#' anim
+#'
+#' ### Project a distribution
+#'
+#' # Make starting distribution
+#' # Since we define the point in WGS84 (not crs(bf)) we also have to provide
+#' # the crs.
+#' point <- data.frame(x = -90, y = 35)
+#' d1 <- as_distr(point, bf, crs = "EPSG:4326" )
+#'
+#' # Project - density will spread over type resulting in a vastly different
+#' # range of values
+#' density_spread <- predict(bf, d1, season = "prebreeding")
+#'
+#' # Have the color gradient rescaled to the range of data in each
+#' # individual frame  - density scaling is dynamic.
+#' spread_anim <- animate_distr(density_spread, bf,   dynamic_scale= TRUE)
+#'
+#' # Or put in values to use for the limits of the color scale - values outside
+#' # of the limits will be truncated
+#' spread_anim <- animate_distr(density_spread, bf,  limit = c(0, 0.05))
+#'
+#'
+#' \dontrun{
+#'   # example render fo file
+#'   gif <- gganimate::animate(spread_anim,
+#'                             device = "ragg_png", # ragg_png is fast and pretty
+#'                             width = 7, height = 6,
+#'                             res = 150, units = "in")
+#'   # Display
+#'   print(gif)
+#'
+#'   # mv
+#'   gif_file <- tempfile("animation", fileext = ".gif")
+#'   gganimate::save_animation(gif, gif_file)
+#'   file.remove(gif_file) # cleanup
+#'   }
+animate_distr <- function(distr, bf, title = species(bf), ...){
+
+   p <- plot_distr(distr, bf, ...)
+
+  # Drop faceting and add animation
+  a <- p +
+    ggplot2::facet_null() +
+    gganimate::transition_manual(frames = .data$label) +
+    ggplot2::labs(title = title,
+                  subtitle = "{current_frame}")
+
+  return(a)
+}
diff --git a/R/is_location_valid.R b/R/is_location_valid.R
index 1614a254..8f56c485 100644
--- a/R/is_location_valid.R
+++ b/R/is_location_valid.R
@@ -10,7 +10,7 @@
 #' one should be used: either `i` a vector of location indices
 #' (see [xy_to_i()]), or `x` and `y`.
 #'
-#' For both functions time shold be input either as `timestep` or `date`.
+#' For both functions time should be input either as `timestep` or `date`.
 #'
 #' The number of timesteps or dates should either be 1 or match to the number
 #' of locations or distributions; if singular it will be applied to all.
diff --git a/R/plot_distr.R b/R/plot_distr.R
new file mode 100644
index 00000000..3ab8fbaf
--- /dev/null
+++ b/R/plot_distr.R
@@ -0,0 +1,278 @@
+
+#' Plot distributions
+#'
+#' Return a [ggplot2::ggplot] object with plots of one or more
+#' distributions.
+#'
+#' @details
+#' ## Masks
+#'
+#' BirdFlow objects have both a mask and a dynamic
+#' mask. For any given timestep the model only covers cells that aren't excluded
+#' by either of these masks.
+#'
+#' * **mask** The static mask, AKA "mask", indicates which cells
+#' in the raster extent are included in the model at any timestep these are
+#' also called active cells (e.g. [n_active()]). The number of elements in a
+#' distribution always match the number of unmasked cells in the (static) mask.
+#' Cells are masked if they always have zero abundance in the Bird S&T data
+#' for every week of the year.
+#'
+#' * **dynamic mask** The dynamic mask indicates which of the cells not excluded
+#' by the static mask are modeled within each timestep. It,like distributions,
+#' only has values for active cells. Thus the dimensions of objects returned
+#' by [get_dynamic_mask(bf)](get_dynamic_mask()) and get
+#' [get_distr(bf)](get_distr()) will be identical.
+#' The purposoe of the dynamic mask is to improve efficiency of fitting,
+#' storing, and using BirdFlow models by eliminating unlikely locations
+#' in the model.  The dynamic mask includes cells that have a zero in the
+#' eBird S&T abundance for the associated week.
+#'
+#' To display both masks set `show_dynamic_mask` to `TRUE` and leave
+#' `show_mask` at it's default (`TRUE`).  Showing the dynamic mask relies on
+#' matching the columns names of the distribution to timesteps in the model so
+#' requires that the column names match the column names in `get_distr()`.
+#'
+#' ## Subtitles
+#' Plot subtitles will be derived from the column names in `distr`. It may be
+#' useful to define single distributions as a one column matrix so that the
+#' label information is there. For example if `d` is a distribution in a
+#' matrix with multiple columns use `drop=FALSE` while subsetting,
+#' `plot_distr(d[, 1, drop = FALSE], bf)`; and to add a label to a vector
+#' distribution use `d <- matrix(d, ncol = 1);  colnames(d) <- "new label"`
+#'
+#' @param distr a vector or matrix with distribution values corresponding to
+#'   active cells in `bf`.
+#' @param bf A BirdFlow object
+#' @param subset Defines an optional subset of `distr` that should be plotted.
+#'   Use either column numbers, column names, or a logical vector.
+#' @param show_mask If `TRUE` (the default) the static mask that indicates which
+#'   cells are active in the model at any timestep will be shown.
+#' @param show_dynamic_mask Defaults to `FALSE`. Set to `TRUE` to visualize the
+#'   dynamic mask. This is achieved by overwriting cells that are dynamically
+#'   masked with NA. For `show_dynamic_mask = TRUE` to work the column names in
+#'   `distr` should all be in `colnames(get_distr(bf))`.  This is
+#'   true for distributions returned by [`predict()`](predict.BirdFlowR)
+#'   and [get_distr()].
+#' @param limits The range of density values to map `gradient_colors` to. The
+#'   default is the range of the values in `distr` after applying `subset`. If
+#'   you want to standardize the range across multiple models of a single
+#'   species you might want to set to `c(0, max)` where `max` is the maximum
+#'   observed value across all models. Alternatively if the range is highly
+#'   variable among the columns in ``distr`` as when density spreads out from a
+#'   single point in the results of [`predict(bf)`](predict.BirdFlowR) you may
+#'   want to set this smaller than the full range in which case the values will
+#'   be trunctated to the limits (see examples).
+#' @param dynamic_scale Set to `TRUE` to have the range of the data in each
+#'   distribution mapped to the full color gradient. This allows visualizing the
+#'   full range of values within each timestep optimally at the cost of
+#'   consistency in the color scale among the facets - or animation frames
+#'   if using [animate_distr()].
+#' @param coast_linewidth The line width to use when plotting the coast. Default
+#'   is `0.25`.  If `NULL` the coast will not be plotted.
+#' @param coast_color The color to use for plotting the coastline. If `NULL` the
+#'   coast will not be plotted.
+#' @param gradient_colors A color gradient that will be used to plot the density
+#'   values. Leave or set to `NULL` for the default of
+#'   `ebirdst::abundance_palette(10, season = "weekly")`.
+#' @param active_cell_color The background color for active cells in the
+#'   landscape. Only used if `show_mask` is `TRUE`. These cells will only be
+#'   visible if there are `NA` values in `distr` or if `show_dynamic_mask`
+#'   is `TRUE`.
+#' @param inactive_cell_color The color to use for inactive cells in the landscape.
+#'   These are cells that are always masked. Only relevant when
+#'   `show_mask = TRUE`.
+#' @param title The title for the plot. It defaults to the common name of the
+#'   species (`species(bf)`).
+#' @param value_label The label used for the values in the distribution.
+#'   Defaults to "Density"
+#'
+#'
+#' @return [ggplot2::ggplot()] object.  Use `print()` to render it.
+#' @export
+#' @importFrom grDevices gray grey
+#' @examples
+#' bf <- BirdFlowModels::amewoo
+#' p <- plot_distr(get_distr(bf, c(1,11, 21)), bf, show_dynamic_mask = TRUE)
+#' @seealso
+#' * [animate_distr()] for animating distributions.
+#' * [plot_routes()] and  [animate_routes()] for visualizing routes.
+#' * [as_distr()], [get_distr()], [`predict()`](predict.BirdFlow)
+#' for functions that produce distributions.
+plot_distr <- function(distr,
+                       bf,
+                       subset = NULL,
+                       show_mask = TRUE,
+                       show_dynamic_mask = FALSE,
+                       limits = NULL,
+                       dynamic_scale = FALSE,
+                       coast_linewidth =  .25,
+                       coast_color = gray(0.5),
+                       gradient_colors = NULL,
+                       active_cell_color = rgb(1, 1, 1, .3),
+                       inactive_cell_color = rgb(0, 0, 0, .2),
+                       title = species(bf),
+                       value_label = "Density"){
+
+  if(dynamic_scale){
+    distr <- apply(distr, 2, function(x) x / max(x, na.rm = TRUE) )
+  }
+
+
+  if(!missing(subset)){
+    if(is.null(dim(distr)) || length(dim(distr)) != 2 ){
+      stop("subset should only be used if distr contains multiple distributions")
+    }
+
+    # Apply subset
+    distr <- distr[, subset, drop = FALSE]
+
+    # If a single distr drop dimensions and add attributes
+    if(ncol(distr) == 1){
+      name <- colnames(distr)
+      distr <- as.vector(distr)
+      attr(distr, "time")<- name
+    }
+  }
+
+
+  if(is.null(limits)){
+    limits <- range(distr, na.rm = TRUE)
+  } else {
+    stopifnot(is.numeric(limits), length(limits) == 2, all(!is.na(limits)),
+              limits[1] < limits[2])
+    # Truncate to limits
+    distr[distr < limits[1]] <- limits[1]
+    distr[distr > limits[2]] <- limits[2]
+  }
+
+  if(show_dynamic_mask){
+    bf <- add_dynamic_mask(bf)
+    dm <- get_dynamic_mask(bf)
+
+    multiple <- is.matrix(distr) || is.array(distr)
+    if(multiple){
+      if(!all(colnames(distr) %in% colnames(dm)))
+        stop(" Cannot apply dynamic mask. ",
+             "Not all column labels in distr match time labels from bf.")
+      dm <- dm[, match(colnames(distr), colnames(dm))]
+      stopifnot(all(colnames(distr) == colnames(dm)))
+      distr[!dm] <- NA
+    } else {
+      label <- attr(distr, "time")
+      if(is.null(label) || !label %in% colnames(dm))
+        stop(" Cannot apply dynamic mask. ",
+             "distr label is missing or not a time label from bf")
+
+      dm <- dm[, colnames(dm) == label]
+      distr[!dm] <- NA
+    }
+  }
+
+
+  r <- rasterize_distr(distr, bf, format = "dataframe" )
+  multiple <- length(unique(r$order)) > 1
+
+  names(r)[names(r) == "value"] <- value_label
+
+  # Add original order (column index) to the data and define a labeller tp
+  # convert that to the column names.
+  # Note: if you use the column names directly to create the facets then
+  # they end up in alphabetical order, not the order in the the distr
+  # object.
+  if(multiple){
+    order <- data.frame(order = 1:ncol(distr), label = colnames(distr))
+    order_to_label <- function(x) {
+      order$label[match(as.character(x), as.character(order$order))]}
+    order_labeller  <-  ggplot2::as_labeller(order_to_label)
+  }
+
+  coast <- get_coastline(bf)
+
+  if(is.null(gradient_colors))
+    gradient_colors <- ebirdst::abundance_palette(10, season = "weekly")
+
+  p <-
+    ggplot2::ggplot(r, ggplot2::aes(x = .data$x,
+                                    y = .data$y,
+                                    fill = .data[[value_label]])) +
+    ggplot2::geom_raster()
+
+
+  if(dynamic_scale){
+    p <- p +  ggplot2::scale_fill_gradientn(colors = gradient_colors,
+                                            na.value = active_cell_color,
+                                            limits = limits,
+                                            breaks = c(0, 1),
+                                            labels = c("Min.", "Max."))
+  } else {
+    p <- p + ggplot2::scale_fill_gradientn(colors = gradient_colors,
+                                           na.value = active_cell_color,
+                                           limits = limits)
+  }
+
+  if(!is.null(coast_color) && !is.null(coast_linewidth)){
+
+    p  <- p +
+      ggplot2::geom_sf(data = coast,
+                       inherit.aes = FALSE,
+                       linewidth = coast_linewidth,
+                       color = coast_color)
+
+  }
+
+  # coord_sf is required to adjust coordinates while using geom_sf
+  # Here we are preventing expanding the extent of the plot.
+  # Setting the CRS is only necessary when the coastline isn't plotted because
+  # of NULL coast_color or coast_linewidth
+  p <- p +
+    ggplot2::coord_sf(expand = FALSE, crs = crs(bf)) +
+    ggplot2::theme(axis.title = ggplot2::element_blank(),
+                   strip.background = ggplot2::element_blank())
+
+  if(multiple) {
+    p <- p +
+      ggplot2::facet_wrap(facets = ggplot2::vars(order),
+                          labeller = order_labeller ) +
+      ggplot2::ggtitle( label = title )
+
+  } else {
+    subtitle <- r$label[1]
+    p <- p +
+      ggplot2::ggtitle(label = title, subtitle = subtitle)
+
+  }
+
+  if(show_mask){
+    # This adds
+
+    # Note you can't have more than one geom_raster layer in a ggplot, but
+    # you can add a annotation_raster.  It's a matrix with the colors in the
+    # cells - often it's used to overlay images on the plot.
+
+    # I also played with annotate(geom = "raster") and got it to work with
+    # single plots but not facets.
+
+
+    # Make a color matrix for the mask
+    mask <- bf$geom$mask
+    col_mask <- matrix("", nrow(mask), ncol(mask))
+    col_mask[mask] <- active_cell_color
+    col_mask[!mask] <- inactive_cell_color
+
+    # Add it to the plot
+    p <- p + ggplot2::annotation_raster(col_mask, xmin = xmin(bf),
+                               xmax = xmax(bf),
+                               ymin = ymin(bf),
+                               ymax = ymax(bf))
+
+
+    # Move the new annotation layer to the first layer so it draws under others
+    n_layers <- length(p$layers)
+    p$layers <- p$layers[c(n_layers, seq_len(n_layers - 1 ))]
+  }
+
+  return(p)
+
+}
+
diff --git a/R/plot_routes.R b/R/plot_routes.R
index af8138c3..fb2d59b4 100644
--- a/R/plot_routes.R
+++ b/R/plot_routes.R
@@ -178,8 +178,8 @@ plot_routes <- function(routes, bf, facet = FALSE, max_stay_len = NULL,
 
   # Make raster showing which cells are active in the model
   rast <- rasterize_distr(rep(TRUE, n_active(bf)), bf, format = "dataframe")
-  rast$density[is.na(rast$density)] <- FALSE
-  rast$value <- rast$density
+  rast$value[is.na(rast$value)] <- FALSE
+  rast$value <- rast$value
 
   # Coastline for this model
   coast <- get_coastline(bf)
diff --git a/R/predict.R b/R/predict.R
index 3f9f478c..f54cc002 100644
--- a/R/predict.R
+++ b/R/predict.R
@@ -3,7 +3,7 @@
 #' `predict()` projects bird distributions into the future or past. Given an
 #'  initial distribution and time period specified via `...`, `predict()`
 #'  generates probability distributions for each timestep.
-#'
+#' @aliases predict
 #' @param object A BirdFlow model object
 #' @param distr  a starting distribution
 #' @inheritDotParams lookup_timestep_sequence -x
diff --git a/R/rasterize_distr.R b/R/rasterize_distr.R
index c256cc67..6c6bcefb 100644
--- a/R/rasterize_distr.R
+++ b/R/rasterize_distr.R
@@ -1,10 +1,10 @@
 #' @name rasterize
-#' @title Functions to convert a BirdFlow object or distribution into a
-#' SpatRaster
-#' @description `rasterize_distr()` creates a [SpatRaster][terra::SpatRaster]
-#' similar to those created by [terra::rast()] from one or more distributions
-#' in their compact vector form. `rast()` converts a BirdFlow object directly
+#' @title Convert a BirdFlow distribution into a raster
+#' @description `rast()` converts a BirdFlow object directly
 #' to a [SpatRaster][terra::SpatRaster].
+#' `rasterize_distr()` converts a [distribution](as_distr) into a
+#' [SpatRaster][terra::SpatRaster], numeric matrix or array, or a raster data
+#'  frame.
 #'
 #' @param distr a distribution in its vector form or a
 #' matrix in which each column represents a different distribution.
@@ -14,16 +14,36 @@
 #'   `'numeric'` for a matrix or array, or`'dataframe'` for raster data
 #'   suitable for plotting with [ggplot2::geom_raster()]
 #' @inheritParams get_distr
-#' @return A [terra::SpatRaster] object.
+#' @return
+#' For `rasterize_distr()` the return type depends on the `format` argument:
+#' * `"SpatRaster"` (the default) returns a [terra::SpatRaster] object.
+#' * `"numeric"` returns a matrix (one distribution) or array (multiple
+#' distributions).  In either case the first two dimensions will be y (rows),
+#' and x (columns).
+#' * `"dataframe"` will return the raster information in a data frame with a
+#' row for every value (long format) with columns:
+#'   * `x`, `y` the x and y coordinates of the cell center.
+#'   * `i` the location index (in `bf`) of the cell.
+#'   * `label` The label associated with the distribution, taken from column
+#'   names of `distr` - typically it indicates time. It is an ordered factor
+#'    with the level order matching the order of the distribution in `distr`.
+#'    The ordered factor is helpful when animating.
+#'   * `value`  The cell value, typically density
+#'   * `order`  The column index of the distribution in `distr` or if only one
+#'   distribution `1`.
+#'   The object is suitable for plotting with [[ggplot2::geom_raster]].
+#'
+#'
 #' @importMethodsFrom terra rast
 #' @export
-rasterize_distr <- function(distr, bf, format = "SpatRast") {
+rasterize_distr <- function(distr, bf, format = "SpatRaster") {
 
   format <- tolower(format)
 
   # allow alternate output format names
   format <- switch(format,
-                   "terra" = "spatrast",
+                   "spatrast" = "spatraster",
+                   "terra" = "spatraster",
                    "matrix" = "numeric",
                    "array" = "numeric",
                    "raster.data.frame" = "dataframe",
@@ -31,8 +51,8 @@ rasterize_distr <- function(distr, bf, format = "SpatRast") {
                    format
                    )
 
-  stopifnot("Format must be one of 'SpatRast', 'numeric', or 'dataframe'" =
-              format %in% c("spatrast", "numeric", "dataframe"))
+  stopifnot("Format must be one of 'SpatRaster', 'numeric', or 'dataframe'" =
+              format %in% c("spatraster", "numeric", "dataframe"))
 
   if (format == "dataframe") {
     # Data frame for use with ggplot2 geom_raster() it will
@@ -56,24 +76,32 @@ rasterize_distr <- function(distr, bf, format = "SpatRast") {
     if (distr_is_vector) {
       distr2 <- rep(NA, n_rast)
       distr2[mv] <- distr
-      distr2 <- data.frame(density = distr2)
+      distr2 <- data.frame(value = distr2)
       # Add time column if appropriate
       time <- attr(distr, "time")
-      if (!is.null(time))
-        distr2 <- cbind(data.frame(time = time), distr2)
+      if (is.null(time))
+        time <- ""
+      distr2 <- cbind(data.frame(label = time), distr2)
       df <- cbind(all[, c(-1, -2)], distr2)
-    } else {
+      df$label <- ordered(df$label)
+
+    } else {  # multiple distributions
       distr2 <- matrix(NA, nrow = n_rast, ncol = ncol(distr))
       distr2[mv, ] <- distr
       colnames(distr2) <- colnames(distr)
       distr2 <- as.data.frame(distr2)
       df <- cbind(all[, c(-1, -2)], distr2)
       df <-   tidyr::pivot_longer(df,  cols = setdiff(names(df), names(all)),
-                                  names_to = "time",
-                                  values_to = "density")
-    }
+                                  names_to = "label",
+                                  values_to = "value")
+      # Label column has original distribution column names as a an ordered
+      # factor with the level in the
+      df$label <- ordered(df$label, levels = colnames(distr))
 
+    }
 
+    # 1:n_distr  in original column order
+    df$order <- as.numeric(df$label)
 
     return(df)
   }
@@ -84,7 +112,7 @@ rasterize_distr <- function(distr, bf, format = "SpatRast") {
     return(m)
   }
 
-  if (format == "spatrast") {
+  if (format == "spatraster") {
     r <- terra::rast(m, extent = bf$geom$ext, crs = bf$geom$crs)
     n_dim <- length(dim(m)) # no of dimensions of output
     if (n_dim == 3) # two or more distributions
@@ -109,5 +137,6 @@ rast.BirdFlow <- function(x, which = "all") {
 }
 # nolint end
 #' @rdname rasterize
+#' @return `rast()` returns a [terra::SpatRaster]
 #' @export
 setMethod(rast, "BirdFlow", rast.BirdFlow)
diff --git a/man/animate_distr.Rd b/man/animate_distr.Rd
new file mode 100644
index 00000000..bf287a40
--- /dev/null
+++ b/man/animate_distr.Rd
@@ -0,0 +1,121 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/animate_distr.R
+\name{animate_distr}
+\alias{animate_distr}
+\title{Animate distributions}
+\usage{
+animate_distr(distr, bf, title = species(bf), ...)
+}
+\arguments{
+\item{distr}{A set of distributions it should be a matrix with \code{n_active(bf)}
+rowsand a column for each distribution. The animation will proceed in the
+in column order and column names will be used a subtitles in the plot.}
+
+\item{bf}{The BirdFlow object that \code{distr} is associated with.}
+
+\item{title}{The title to use for the animation. The default is
+the common name of the species.}
+
+\item{...}{
+  Arguments passed on to \code{\link[=plot_distr]{plot_distr}}
+  \describe{
+    \item{\code{subset}}{Defines an optional subset of \code{distr} that should be plotted.
+Use either column numbers, column names, or a logical vector.}
+    \item{\code{show_mask}}{If \code{TRUE} (the default) the static mask that indicates which
+cells are active in the model at any timestep will be shown.}
+    \item{\code{show_dynamic_mask}}{Defaults to \code{FALSE}. Set to \code{TRUE} to visualize the
+dynamic mask. This is achieved by overwriting cells that are dynamically
+masked with NA. For \code{show_dynamic_mask = TRUE} to work the column names in
+\code{distr} should all be in \code{colnames(get_distr(bf))}.  This is
+true for distributions returned by \href{predict.BirdFlowR}{\code{predict()}}
+and \code{\link[=get_distr]{get_distr()}}.}
+    \item{\code{limits}}{The range of density values to map \code{gradient_colors} to. The
+default is the range of the values in \code{distr} after applying \code{subset}. If
+you want to standardize the range across multiple models of a single
+species you might want to set to \code{c(0, max)} where \code{max} is the maximum
+observed value across all models. Alternatively if the range is highly
+variable among the columns in \code{distr} as when density spreads out from a
+single point in the results of \href{predict.BirdFlowR}{\code{predict(bf)}} you may
+want to set this smaller than the full range in which case the values will
+be trunctated to the limits (see examples).}
+    \item{\code{dynamic_scale}}{Set to \code{TRUE} to have the range of the data in each
+distribution mapped to the full color gradient. This allows visualizing the
+full range of values within each timestep optimally at the cost of
+consistency in the color scale among the facets - or animation frames
+if using \code{\link[=animate_distr]{animate_distr()}}.}
+    \item{\code{coast_linewidth}}{The line width to use when plotting the coast. Default
+is \code{0.25}.  If \code{NULL} the coast will not be plotted.}
+    \item{\code{coast_color}}{The color to use for plotting the coastline. If \code{NULL} the
+coast will not be plotted.}
+    \item{\code{gradient_colors}}{A color gradient that will be used to plot the density
+values. Leave or set to \code{NULL} for the default of
+\code{ebirdst::abundance_palette(10, season = "weekly")}.}
+    \item{\code{active_cell_color}}{The background color for active cells in the
+landscape. Only used if \code{show_mask} is \code{TRUE}. These cells will only be
+visible if there are \code{NA} values in \code{distr} or if \code{show_dynamic_mask}
+is \code{TRUE}.}
+    \item{\code{inactive_cell_color}}{The color to use for inactive cells in the landscape.
+These are cells that are always masked. Only relevant when
+\code{show_mask = TRUE}.}
+    \item{\code{value_label}}{The label used for the values in the distribution.
+Defaults to "Density"}
+  }}
+}
+\value{
+A gganimate object that can be displayed with \code{print()} or
+or \code{gganimate::animate()}.  See example for how to export to a file.
+}
+\description{
+Animate distributions as produced by \code{get_distr()},
+\code{as_distr()}, or \code{predict()}.  The distributions will be displayed
+in the column order in \code{distr} and column labels will be used as
+plot subtitle.
+}
+\examples{
+
+# Animate distributions from BirdFlow object - derived from
+# eBird Status and Trends:
+
+bf <- BirdFlowModels::amewoo
+ts <- lookup_timestep_sequence(bf, season = "prebreeding")
+distr <- get_distr(bf, ts)
+anim <- animate_distr(distr, bf,  show_dynamic_mask = TRUE)
+
+anim
+
+### Project a distribution
+
+# Make starting distribution
+# Since we define the point in WGS84 (not crs(bf)) we also have to provide
+# the crs.
+point <- data.frame(x = -90, y = 35)
+d1 <- as_distr(point, bf, crs = "EPSG:4326" )
+
+# Project - density will spread over type resulting in a vastly different
+# range of values
+density_spread <- predict(bf, d1, season = "prebreeding")
+
+# Have the color gradient rescaled to the range of data in each
+# individual frame  - density scaling is dynamic.
+spread_anim <- animate_distr(density_spread, bf,   dynamic_scale= TRUE)
+
+# Or put in values to use for the limits of the color scale - values outside
+# of the limits will be truncated
+spread_anim <- animate_distr(density_spread, bf,  limit = c(0, 0.05))
+
+
+\dontrun{
+  # example render fo file
+  gif <- gganimate::animate(spread_anim,
+                            device = "ragg_png", # ragg_png is fast and pretty
+                            width = 7, height = 6,
+                            res = 150, units = "in")
+  # Display
+  print(gif)
+
+  # mv
+  gif_file <- tempfile("animation", fileext = ".gif")
+  gganimate::save_animation(gif, gif_file)
+  file.remove(gif_file) # cleanup
+  }
+}
diff --git a/man/is_location_valid.Rd b/man/is_location_valid.Rd
index fc96bfc0..f5512f39 100644
--- a/man/is_location_valid.Rd
+++ b/man/is_location_valid.Rd
@@ -44,7 +44,7 @@ For \code{is_location_valid()} there are 2 ways of inputting locations. Only
 one should be used: either \code{i} a vector of location indices
 (see \code{\link[=xy_to_i]{xy_to_i()}}), or \code{x} and \code{y}.
 
-For both functions time shold be input either as \code{timestep} or \code{date}.
+For both functions time should be input either as \code{timestep} or \code{date}.
 
 The number of timesteps or dates should either be 1 or match to the number
 of locations or distributions; if singular it will be applied to all.
diff --git a/man/plot_distr.Rd b/man/plot_distr.Rd
new file mode 100644
index 00000000..275f6e26
--- /dev/null
+++ b/man/plot_distr.Rd
@@ -0,0 +1,142 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_distr.R
+\name{plot_distr}
+\alias{plot_distr}
+\title{Plot distributions}
+\usage{
+plot_distr(
+  distr,
+  bf,
+  subset = NULL,
+  show_mask = TRUE,
+  show_dynamic_mask = FALSE,
+  limits = NULL,
+  dynamic_scale = FALSE,
+  coast_linewidth = 0.25,
+  coast_color = gray(0.5),
+  gradient_colors = NULL,
+  active_cell_color = rgb(1, 1, 1, 0.3),
+  inactive_cell_color = rgb(0, 0, 0, 0.2),
+  title = species(bf),
+  value_label = "Density"
+)
+}
+\arguments{
+\item{distr}{a vector or matrix with distribution values corresponding to
+active cells in \code{bf}.}
+
+\item{bf}{A BirdFlow object}
+
+\item{subset}{Defines an optional subset of \code{distr} that should be plotted.
+Use either column numbers, column names, or a logical vector.}
+
+\item{show_mask}{If \code{TRUE} (the default) the static mask that indicates which
+cells are active in the model at any timestep will be shown.}
+
+\item{show_dynamic_mask}{Defaults to \code{FALSE}. Set to \code{TRUE} to visualize the
+dynamic mask. This is achieved by overwriting cells that are dynamically
+masked with NA. For \code{show_dynamic_mask = TRUE} to work the column names in
+\code{distr} should all be in \code{colnames(get_distr(bf))}.  This is
+true for distributions returned by \href{predict.BirdFlowR}{\code{predict()}}
+and \code{\link[=get_distr]{get_distr()}}.}
+
+\item{limits}{The range of density values to map \code{gradient_colors} to. The
+default is the range of the values in \code{distr} after applying \code{subset}. If
+you want to standardize the range across multiple models of a single
+species you might want to set to \code{c(0, max)} where \code{max} is the maximum
+observed value across all models. Alternatively if the range is highly
+variable among the columns in \code{distr} as when density spreads out from a
+single point in the results of \href{predict.BirdFlowR}{\code{predict(bf)}} you may
+want to set this smaller than the full range in which case the values will
+be trunctated to the limits (see examples).}
+
+\item{dynamic_scale}{Set to \code{TRUE} to have the range of the data in each
+distribution mapped to the full color gradient. This allows visualizing the
+full range of values within each timestep optimally at the cost of
+consistency in the color scale among the facets - or animation frames
+if using \code{\link[=animate_distr]{animate_distr()}}.}
+
+\item{coast_linewidth}{The line width to use when plotting the coast. Default
+is \code{0.25}.  If \code{NULL} the coast will not be plotted.}
+
+\item{coast_color}{The color to use for plotting the coastline. If \code{NULL} the
+coast will not be plotted.}
+
+\item{gradient_colors}{A color gradient that will be used to plot the density
+values. Leave or set to \code{NULL} for the default of
+\code{ebirdst::abundance_palette(10, season = "weekly")}.}
+
+\item{active_cell_color}{The background color for active cells in the
+landscape. Only used if \code{show_mask} is \code{TRUE}. These cells will only be
+visible if there are \code{NA} values in \code{distr} or if \code{show_dynamic_mask}
+is \code{TRUE}.}
+
+\item{inactive_cell_color}{The color to use for inactive cells in the landscape.
+These are cells that are always masked. Only relevant when
+\code{show_mask = TRUE}.}
+
+\item{title}{The title for the plot. It defaults to the common name of the
+species (\code{species(bf)}).}
+
+\item{value_label}{The label used for the values in the distribution.
+Defaults to "Density"}
+}
+\value{
+\code{\link[ggplot2:ggplot]{ggplot2::ggplot()}} object.  Use \code{print()} to render it.
+}
+\description{
+Return a \link[ggplot2:ggplot]{ggplot2::ggplot} object with plots of one or more
+distributions.
+}
+\details{
+\subsection{Masks}{
+
+BirdFlow objects have both a mask and a dynamic
+mask. For any given timestep the model only covers cells that aren't excluded
+by either of these masks.
+\itemize{
+\item \strong{mask} The static mask, AKA "mask", indicates which cells
+in the raster extent are included in the model at any timestep these are
+also called active cells (e.g. \code{\link[=n_active]{n_active()}}). The number of elements in a
+distribution always match the number of unmasked cells in the (static) mask.
+Cells are masked if they always have zero abundance in the Bird S&T data
+for every week of the year.
+\item \strong{dynamic mask} The dynamic mask indicates which of the cells not excluded
+by the static mask are modeled within each timestep. It,like distributions,
+only has values for active cells. Thus the dimensions of objects returned
+by \href{get_dynamic_mask()}{get_dynamic_mask(bf)} and get
+\href{get_distr()}{get_distr(bf)} will be identical.
+The purposoe of the dynamic mask is to improve efficiency of fitting,
+storing, and using BirdFlow models by eliminating unlikely locations
+in the model.  The dynamic mask includes cells that have a zero in the
+eBird S&T abundance for the associated week.
+}
+
+To display both masks set \code{show_dynamic_mask} to \code{TRUE} and leave
+\code{show_mask} at it's default (\code{TRUE}).  Showing the dynamic mask relies on
+matching the columns names of the distribution to timesteps in the model so
+requires that the column names match the column names in \code{get_distr()}.
+}
+
+\subsection{Subtitles}{
+
+Plot subtitles will be derived from the column names in \code{distr}. It may be
+useful to define single distributions as a one column matrix so that the
+label information is there. For example if \code{d} is a distribution in a
+matrix with multiple columns use \code{drop=FALSE} while subsetting,
+\code{plot_distr(d[, 1, drop = FALSE], bf)}; and to add a label to a vector
+distribution use \verb{d <- matrix(d, ncol = 1);  colnames(d) <- "new label"}
+}
+}
+\examples{
+bf <- BirdFlowModels::amewoo
+p <- plot_distr(get_distr(bf, c(1,11, 21)), bf, show_dynamic_mask = TRUE)
+}
+\seealso{
+\itemize{
+\item \code{\link[=animate_distr]{animate_distr()}} for animating distributions.
+\item \code{\link[=plot_routes]{plot_routes()}} and  \code{\link[=animate_routes]{animate_routes()}} for visualizing routes.
+\item \code{\link[=as_distr]{as_distr()}}, \code{\link[=get_distr]{get_distr()}}, \href{predict.BirdFlow}{\code{predict()}}
+for functions that produce distributions.
+}
+}
diff --git a/man/predict.BirdFlow.Rd b/man/predict.BirdFlow.Rd
index 828a539c..11dccc2f 100644
--- a/man/predict.BirdFlow.Rd
+++ b/man/predict.BirdFlow.Rd
@@ -2,6 +2,7 @@
 % Please edit documentation in R/predict.R
 \name{predict.BirdFlow}
 \alias{predict.BirdFlow}
+\alias{predict}
 \title{Predict bird distributions}
 \usage{
 \method{predict}{BirdFlow}(object, distr, ...)
diff --git a/man/rasterize.Rd b/man/rasterize.Rd
index c99a5153..33898491 100644
--- a/man/rasterize.Rd
+++ b/man/rasterize.Rd
@@ -4,10 +4,9 @@
 \alias{rasterize}
 \alias{rasterize_distr}
 \alias{rast,BirdFlow-method}
-\title{Functions to convert a BirdFlow object or distribution into a
-SpatRaster}
+\title{Convert a BirdFlow distribution into a raster}
 \usage{
-rasterize_distr(distr, bf, format = "SpatRast")
+rasterize_distr(distr, bf, format = "SpatRaster")
 
 \S4method{rast}{BirdFlow}(x, which = "all")
 }
@@ -29,11 +28,34 @@ e.g. \code{"2019-02-25"}; \code{\link[base:Dates]{Date}} objects;    or \code{"a
 return distributions for all timesteps.}
 }
 \value{
-A \link[terra:SpatRaster-class]{terra::SpatRaster} object.
+For \code{rasterize_distr()} the return type depends on the \code{format} argument:
+\itemize{
+\item \code{"SpatRaster"} (the default) returns a \link[terra:SpatRaster-class]{terra::SpatRaster} object.
+\item \code{"numeric"} returns a matrix (one distribution) or array (multiple
+distributions).  In either case the first two dimensions will be y (rows),
+and x (columns).
+\item \code{"dataframe"} will return the raster information in a data frame with a
+row for every value (long format) with columns:
+\itemize{
+\item \code{x}, \code{y} the x and y coordinates of the cell center.
+\item \code{i} the location index (in \code{bf}) of the cell.
+\item \code{label} The label associated with the distribution, taken from column
+names of \code{distr} - typically it indicates time. It is an ordered factor
+with the level order matching the order of the distribution in \code{distr}.
+The ordered factor is helpful when animating.
+\item \code{value}  The cell value, typically density
+\item \code{order}  The column index of the distribution in \code{distr} or if only one
+distribution \code{1}.
+The object is suitable for plotting with [\link[ggplot2:geom_tile]{ggplot2::geom_raster}].
+}
+}
+
+\code{rast()} returns a \link[terra:SpatRaster-class]{terra::SpatRaster}
 }
 \description{
-\code{rasterize_distr()} creates a \link[terra:SpatRaster-class]{SpatRaster}
-similar to those created by \code{\link[terra:rast]{terra::rast()}} from one or more distributions
-in their compact vector form. \code{rast()} converts a BirdFlow object directly
+\code{rast()} converts a BirdFlow object directly
 to a \link[terra:SpatRaster-class]{SpatRaster}.
+\code{rasterize_distr()} converts a \href{as_distr}{distribution} into a
+\link[terra:SpatRaster-class]{SpatRaster}, numeric matrix or array, or a raster data
+frame.
 }
diff --git a/tests/testthat/_snaps/rasterize_distr.md b/tests/testthat/_snaps/rasterize_distr.md
index 540036cc..1aa5fa9b 100644
--- a/tests/testthat/_snaps/rasterize_distr.md
+++ b/tests/testthat/_snaps/rasterize_distr.md
@@ -3,30 +3,30 @@
     Code
       hdf
     Output
-                   x       y    i      time      density
-      110 -1364968.8 4862048  936 January 4 5.290266e-07
-      155 -1284973.0 4862048  937 January 4 9.646747e-06
-      394  -884994.1 3742106 1442 January 4 3.851673e-07
+                   x       y    i     label        value order
+      110 -1364968.8 4862048  936 January 4 5.290266e-07     1
+      155 -1284973.0 4862048  937 January 4 9.646747e-06     1
+      394  -884994.1 3742106 1442 January 4 3.851673e-07     1
 
 ---
 
     Code
       hdf
     Output
-      # A tibble: 3 x 5
-                x        y     i time           density
-            <dbl>    <dbl> <int> <chr>            <dbl>
-      1 -1364969. 4862048.   936 January 4  0.000000529
-      2 -1364969. 4862048.   936 January 11 0.000000537
-      3 -1364969. 4862048.   936 January 18 0.000000666
+      # A tibble: 3 x 6
+                x        y     i label            value order
+            <dbl>    <dbl> <int> <ord>            <dbl> <dbl>
+      1 -1364969. 4862048.   936 January 4  0.000000529     1
+      2 -1364969. 4862048.   936 January 11 0.000000537     2
+      3 -1364969. 4862048.   936 January 18 0.000000666     3
 
 # rasterize_distr() to dataframe works
 
     Code
       d
     Output
-                   x       y    i      time      density
-      110 -1364968.8 4862048  936 January 4 5.290266e-07
-      155 -1284973.0 4862048  937 January 4 9.646747e-06
-      394  -884994.1 3742106 1442 January 4 3.851673e-07
+                   x       y    i     label        value order
+      110 -1364968.8 4862048  936 January 4 5.290266e-07     1
+      155 -1284973.0 4862048  937 January 4 9.646747e-06     1
+      394  -884994.1 3742106 1442 January 4 3.851673e-07     1
 
diff --git a/tests/testthat/test-animate_distr.R b/tests/testthat/test-animate_distr.R
new file mode 100644
index 00000000..64b86d3c
--- /dev/null
+++ b/tests/testthat/test-animate_distr.R
@@ -0,0 +1,15 @@
+test_that("animate_distr() works with default args", {
+  bf <- BirdFlowModels::amewoo
+
+  expect_no_error(a1 <- animate_distr(get_distr(bf, c(1, 10, 20)), bf) )
+  expect_no_error(print(a1))
+
+})
+
+test_that("animate_distr() works with dynamic masking", {
+  bf <- BirdFlowModels::amewoo
+
+  expect_no_error(a2 <- animate_distr(get_distr(bf, c(1,10, 20)), bf,
+                                      show_dynamic_mask = TRUE) )
+  expect_no_error(print(a2))
+})
diff --git a/tests/testthat/test-plot_distr.R b/tests/testthat/test-plot_distr.R
new file mode 100644
index 00000000..e712abc1
--- /dev/null
+++ b/tests/testthat/test-plot_distr.R
@@ -0,0 +1,70 @@
+test_that("plot_distr() works with default parameters", {
+  bf <- BirdFlowModels::amewoo
+
+  ## Default parameters
+
+  # 1. single plot
+  expect_no_error(p1 <- plot_distr(get_distr(bf, 1), bf) )
+  expect_no_error(print(p1))
+
+  # 2 facet
+  expect_no_error(p2 <- plot_distr(get_distr(bf, 1:3), bf) )
+  expect_no_error(print(p2))
+})
+
+test_that("plot_distr() works with dynamic mask", {
+  bf <- BirdFlowModels::amewoo
+
+  ## With dynamic mask
+
+  # 3. single plot
+  expect_no_error(p3 <- plot_distr(get_distr(bf, 1), bf, show_dynamic_mask = TRUE) )
+  expect_no_error(print(p3))
+
+  # 4 facet
+  expect_no_error(p4 <- plot_distr(get_distr(bf, 1:3), bf, show_dynamic_mask = TRUE) )
+  expect_no_error(print(p4))
+
+})
+
+test_that("plot_distr() works with llimits and dynamic scaling", {
+  bf <- BirdFlowModels::amewoo
+
+  ### Using a prediction spread so range of data varies across distributions
+  point <- data.frame(x = -90, y = 35)
+  d1 <- as_distr(point, bf, crs = "EPSG:4326" )
+  distr <- predict(bf, d1, season = "prebreeding")
+  # Select ~1/3 of the weeks
+  sel_col <- seq_len(ncol(distr))[seq_len(ncol(distr)) %% 3 == 0]
+  distr <- distr[ , sel_col, drop = FALSE]
+
+  # 5. facet with limits
+  expect_no_error(p5 <- plot_distr(distr, bf,
+                                   show_dynamic_mask = TRUE,
+                                   limits = c(0, 0.02)) )
+  expect_no_error(print(p5))
+
+  # 6. facet with dynamic dynamic scaling (of density values)
+
+  expect_no_error(p6 <- plot_distr(distr, bf,
+                                   show_dynamic_mask = TRUE,
+                                   dynamic_scale = TRUE))
+  expect_no_error(print(p6))
+})
+
+test_that("plot_distr() works with custom colors and custom value label and subset", {
+  # This test is somewhat of a catch all for remaining minor options
+  bf <- BirdFlowModels::amewoo
+
+  expect_no_error(p7 <- plot_distr(get_distr(bf, 1:4), bf,
+                                   subset = 2,
+                                   gradient_colors = c("red","yellow", "blue"),
+                                   show_dynamic_mask = TRUE,
+                                   value_label = "Red-yellow-blue"))
+
+  expect_no_error(print(p7))
+
+
+})
+
+
diff --git a/tests/testthat/test-rasterize_distr.R b/tests/testthat/test-rasterize_distr.R
index cfc81d3f..d6c672af 100644
--- a/tests/testthat/test-rasterize_distr.R
+++ b/tests/testthat/test-rasterize_distr.R
@@ -30,18 +30,18 @@ test_that("rasterize_distr() with data.frame output", {
   expect_no_error(df <- rasterize_distr(get_distr(bf, 1), bf, format = "data.frame"))
   expect_false(any(is.na(df$x)))
   expect_false(any(is.na(df$y)))
-  hdf <- head(df[!is.na(df$density) & df$density != 0, ], 3)
+  hdf <- head(df[!is.na(df$value) & df$value != 0, ], 3)
   expect_snapshot(hdf)
   # Convert back to distribution and compare
   vals <- df[!is.na(df$i), ]
-  vals <- vals$density[order(vals$i)]
+  vals <- vals$value[order(vals$i)]
   expect_equal(vals, as.numeric(get_distr(bf,1)))
 
   # Multiple distributions
-  expect_no_error(df <- rasterize_distr(get_distr(bf, 1:3), bf, format = "data.frame"))
+  expect_no_error(df <- rasterize_distr(get_distr(bf, 1:3), bf, format = "dataframe"))
   expect_false(any(is.na(df$x)))
   expect_false(any(is.na(df$y)))
-  hdf <- head(df[!is.na(df$density) & df$density != 0, ], 3)
+  hdf <- head(df[!is.na(df$value) & df$value != 0, ], 3)
   expect_snapshot(hdf)
 })
 
@@ -80,7 +80,7 @@ test_that("rasterize_distr() to dataframe works", {
   d <- get_distr(bf, 1)
   expect_no_error(df <- rasterize_distr(d, bf = bf, format =  "dataframe"))
 
-  d <- head(df[!df$density == 0 & !is.na(df), ], 3)
+  d <- head(df[!df$value == 0 & !is.na(df$value), ], 3)
   expect_snapshot(d)
 })