diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index a3ac6182..8c9e8da1 100644
--- a/.github/workflows/R-CMD-check.yaml
+++ b/.github/workflows/R-CMD-check.yaml
@@ -22,7 +22,7 @@ jobs:
           - {os: windows-latest, r: 'release'}
           - {os: ubuntu-latest,   r: 'devel', http-user-agent: 'release'}
           - {os: ubuntu-latest,   r: 'release'}
-          - {os: ubuntu-latest,   r: 'oldrel-1'}
+#          - {os: ubuntu-latest,   r: 'oldrel-1'}
 
     env:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
diff --git a/NEWS.md b/NEWS.md
index bb4f792a..6c5d1b45 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,30 @@
+
+# Birdflow 0.1.0.9057
+2024-04-17
+
+* `import_birdflow()` now works with preprocessed hdf5 files. See #177
+* `preprocess_birdflow()` and `import_birdflow()` are now more consistent 
+  with some dimension names:
+    * `bf$geom$ext` is now an unnamed vector in objects from either function.
+    Previously one was a named vector. The `ext(bf)` return object is unchanged.
+    * `distr`, `dynamic_mask`, `uci`, and `lci` now always have `NULL` 
+    rownames and 
+    `"t<timestep>"` column names, and the names
+    of those two dimensions are `"i"` and `"time"`.
+    * `import_birdflow()` no longer drops two metadata items that it was losing
+    before: `ebirdst_version`  and `birdflowr_preprocess_version` these 
+    will exist in old `.hdf5` files but not in old `.rds` files.  
+* `compare_lists()` (internal helper) return format changed slightly.
+* `import_birdflow()` and `export_birdflow()` now work with sparse models.
+   Marginals from these models are converted to standard matrices before 
+   writing so some of the benefits of sparsification are lost when writing to
+   hdf5 files - but compression will probably mitigate this somewhat. Sparse
+   Matrices are re-created when sparse models are re-imported.
+   
+
 # BirdFlowR 0.1.0.9056
 2024-04-04
+
 ## Flux III 
 * `is_between()` now uses a `SparseArray::SparseArray()` for the logical array.
   to reduce memory usage
diff --git a/R/clean_hdf5_dataframe.R b/R/clean_hdf5_dataframe.R
new file mode 100644
index 00000000..fc85be06
--- /dev/null
+++ b/R/clean_hdf5_dataframe.R
@@ -0,0 +1,12 @@
+
+# This strips weird extra within column attributes from dataframes returned
+# by h5read
+# It's used by import_birdflow_v3.
+clean_hdf5_dataframe <- function(df) {
+  if (!inherits(df, "data.frame"))
+    return(df)
+
+  for (i in seq_len(ncol(df)))
+    df[[i]] <- as.vector(df[[i]])
+ return(df)
+}
diff --git a/R/compare_list_item_names.R b/R/compare_list_item_names.R
index 9486285c..8062ac65 100644
--- a/R/compare_list_item_names.R
+++ b/R/compare_list_item_names.R
@@ -34,14 +34,13 @@ compare_list_item_names <- function(x, y, map = "x", differences) {
       differences <- rbind(
         differences,
         data.frame(where = map,
-                   differences = paste0("extra:",
-                                        paste0(lost, collapse = ", "))))
+                   differences = paste0("extra:", lost)))
+
     if (length(gained) != 0)
       differences <- rbind(
         differences,
         data.frame(where = map,
-                   differences = paste0("missing:",
-                                        paste0(gained, collapse = ", "))))
+                   differences = paste0("missing:", gained)))
     return(differences)
   }
 
diff --git a/R/export_birdflow.R b/R/export_birdflow.R
index 623103ae..68201e7c 100644
--- a/R/export_birdflow.R
+++ b/R/export_birdflow.R
@@ -58,6 +58,14 @@ export_birdflow <- function(bf, file = NULL,
     return()
   }
 
+  # Convert sparse matrices to standard (for hdf5 only)
+  if (has_marginals(bf) && bf$metadata$is_sparse) {
+    mn <- setdiff(names(bf$marginals), "index")
+    for (m in mn) {
+      bf$marginals[[m]] <- as.matrix(bf$marginals[[m]])
+    }
+  }
+
   # Write HDF5
   ns <- names(bf)
   for (i in seq_along(ns)) {
diff --git a/R/import_birdflow_v3.R b/R/import_birdflow_v3.R
index 18c60511..80edcc11 100644
--- a/R/import_birdflow_v3.R
+++ b/R/import_birdflow_v3.R
@@ -49,7 +49,6 @@ import_birdflow_v3 <- function(hdf5) {
     "geom/crs",
     "geom/mask",
     "transitions",
-    "marginals",
     "dates",
     "distr",
     "species",
@@ -72,8 +71,6 @@ import_birdflow_v3 <- function(hdf5) {
     "metadata/has_marginals",
     "metadata/has_transitions",
     "metadata/has_distr",
-    "metadata/hyperparameters",
-    "metadata/loss_values",
     "metadata/n_transitions",
     "metadata/n_active",
     "metadata/n_timesteps",
@@ -82,13 +79,26 @@ import_birdflow_v3 <- function(hdf5) {
     "metadata/ebird_access_end_date",
     "metadata/birdflow_preprocess_date",
     "metadata/birdflow_model_date",
-    "metadata/is_sparse",
-    "marginals")
+    "metadata/is_sparse"
+  )
+
+  fit_model_items <- c(
+    "metadata/hyperparameters",
+    "metadata/loss_values",
+    "marginals"
+  )
+
 
   # Check HDF5 for version consistency and missing contents
   contents <- h5ls(hdf5)
   contents <- paste0(contents$group, "/", contents$name)
   contents <- gsub("^/*", "", contents)
+
+  is_fitted_model <- "marginals" %in% contents
+  if (is_fitted_model) {
+    expected_contents <- c(expected_contents, fit_model_items)
+  }
+
   absent <- setdiff(expected_contents, contents)
   extra <- setdiff(contents, expected_contents) # nolint: object_usage_linter
 
@@ -97,6 +107,7 @@ import_birdflow_v3 <- function(hdf5) {
          paste(absent, collapse = "', '"), "'")
   }
 
+
   expected_version <- 3  # of HDF5 BirdFlow export
   version <- as.vector(h5read(hdf5, "metadata/birdflow_version"))
   if (version != expected_version) {
@@ -127,31 +138,34 @@ import_birdflow_v3 <- function(hdf5) {
   }
 
   # hyperparameters
-  hp <- h5read(hdf5, "metadata/hyperparameters")
-  # hdf5 seems to store logical as a factor or at least R reads them as such.
-  # The code below looks for factors that store logical values and
-  # explicitly converts them to logical
-  for (i in seq_along(hp)) {
-    a <- hp[[i]] # this hyper parameter
-    if (is.factor(a) && all(tolower(levels(a)) %in% c("true", "false"))) {
-      a <- as.logical(a)
+  if (is_fitted_model) {
+
+    hp <- h5read(hdf5, "metadata/hyperparameters")
+    # hdf5 seems to store logical as a factor or at least R reads them as such.
+    # The code below looks for factors that store logical values and
+    # explicitly converts them to logical
+    for (i in seq_along(hp)) {
+      a <- hp[[i]] # this hyper parameter
+      if (is.factor(a) && all(tolower(levels(a)) %in% c("true", "false"))) {
+        a <- as.logical(a)
+      }
+      if (inherits(a, "array")) {
+        a <- as.vector(a)
+      }
+      hp[[i]] <- a
     }
-    if (inherits(a, "array")) {
-      a <- as.vector(a)
+    bf$metadata$hyperparameters <- hp
+
+    # loss values
+    lv <- as.data.frame(h5read(hdf5, "metadata/loss_values"))
+    for (i in seq_len(ncol(lv))) {
+      # IF R re-exports an imported hdf5 the loss values columns are each
+      # arrays.  This returns them to standard data.frame columns
+      lv[[i]] <- as.vector(lv[[i]])
     }
-    hp[[i]] <- a
-  }
-  bf$metadata$hyperparameters <- hp
-
-  # loss values
-  lv <- as.data.frame(h5read(hdf5, "metadata/loss_values"))
-  for (i in seq_len(ncol(lv))) {
-    # IF R re-exports an imported hdf5 the loss values columns are each
-    # arrays.  This returns them to standard data.frame columns
-    lv[[i]] <- as.vector(lv[[i]])
-  }
-  bf$metadata$loss_values <- lv
+    bf$metadata$loss_values <- lv
 
+  } # end only for fitted models
 
 
   # dates
@@ -162,77 +176,134 @@ import_birdflow_v3 <- function(hdf5) {
   colnames(dates) <- gsub("^week_", "", colnames(dates))
   bf$dates <- dates
 
-  # Save marginals into list
-  marg <- h5read(hdf5, "marginals", native = TRUE)
-  nt <- length(marg[!names(marg) == "index"])
-  bf$metadata$n_transitions <- nt
-  if (is.null(bf$metadata$timestep_padding))
-    bf$metadata$timestep_padding <- nchar(nt)
-  circular <- nt == length(unique(dates$date))
-  bf$marginals <- vector(mode = "list", length = nt)
-
-  # If the hdf5 has been re-exported from R than we just copy the marginals over
-  if ("index" %in% names(marg)) {
-    bf$marginals <- marg
-  } else {
-    # If this hdf5 was written by python then we need to copy and rename
-    # marginals
-    for (i in seq_len(nt)) {
-      python_label <- paste0("Week", i, "_to_", i + 1)
-      if (circular && i == nt) {
-        label <- paste0("M_", pad_timestep(i, bf), "-", pad_timestep(1, bf))
-      } else {
-        label <- paste0("M_", pad_timestep(i, bf), "-", pad_timestep(i + 1, bf))
+  if (is_fitted_model) {
+
+    # Save marginals into list
+    marg <- h5read(hdf5, "marginals", native = TRUE)
+    nt <- length(marg[!names(marg) == "index"])
+    bf$metadata$n_transitions <- nt
+    if (is.null(bf$metadata$timestep_padding))
+      bf$metadata$timestep_padding <- nchar(nt)
+    circular <- nt == length(unique(dates$date))
+    bf$marginals <- vector(mode = "list", length = nt)
+
+    # If the hdf5 has been re-exported from R, just copy the marginals over
+    if ("index" %in% names(marg)) {
+      bf$marginals <- marg
+    } else {
+      # If this hdf5 was written by python then we need to copy and rename
+      # marginals
+      for (i in seq_len(nt)) {
+        python_label <- paste0("Week", i, "_to_", i + 1)
+        if (circular && i == nt) {
+          label <- paste0("M_", pad_timestep(i, bf), "-", pad_timestep(1, bf))
+        } else {
+          label <- paste0("M_", pad_timestep(i, bf), "-",
+                          pad_timestep(i + 1, bf))
+        }
+        bf$marginals[[i]] <- marg[[python_label]]
+        names(bf$marginals)[i] <- label
       }
-      bf$marginals[[i]] <- marg[[python_label]]
-      names(bf$marginals)[i] <- label
+      bf$metadata$has_marginals <- TRUE
     }
-    bf$metadata$has_marginals <- TRUE
   }
 
   # Save distributions
   bf$distr <- h5read(hdf5, "distr", native = TRUE)
 
-  # Cleanup duplicated distribution, dynamic_mask row, and date added to
-  # input to force circular model fitting.
+  if (is_fitted_model) {
+    # Cleanup duplicated distribution, dynamic_mask row, and date added to
+    # input to force circular model fitting.
+
+    # Cleanup duplicated date row
+    sv <- duplicated(bf$dates$date)
+    if (any(sv)) {
+      bf$dates <- bf$dates[!sv, ]
+    }
+    bf$metadata$n_timesteps <- nrow(bf$dates)
+
+    # Delete duplicated distribution
+    d <- bf$distr
+    if (ncol(d) == n_timesteps(bf) + 1) {
+      if (!all(d[, 1] == d[, ncol(d)]))
+        stop("Expected extra distribution to match first distribution")
+      d <- d[, 1:(ncol(d) - 1)]
+    }
+
+    bf$distr <- d
+
+    # Delete duplicated dynamic mask row
+    dm <- bf$geom$dynamic_mask
+    if (ncol(dm) == n_timesteps(bf) + 1) {
+      if (!all(dm[, 1] == dm[, ncol(dm)]))
+        stop("Expected first and last dynamic mask columns to matrch in ",
+             "circular BirdFlow model")
+      dm <- dm[, 1:(ncol(dm) - 1)]
+    }
+
+    bf$geom$dynamic_mask <- dm
+
+    # Make and save marginal index
+    bf$marginals$index <- make_marginal_index(bf)
+
+  } # end fitted model only
 
-  # Cleanup duplicated date row
-  sv <- duplicated(bf$dates$date)
-  if (any(sv)) {
-    bf$dates <- bf$dates[!sv, ]
-  }
-  bf$metadata$n_timesteps <- nrow(bf$dates)
-
-  # Delete duplicated distribution
-  d <- bf$distr
-  if (ncol(d) == n_timesteps(bf) + 1) {
-    if (!all(d[, 1] == d[, ncol(d)]))
-      stop("Expected extra distribution to match first distribution")
-    d <- d[, 1:(ncol(d) - 1)]
-  }
 
-  ### back compatibility code
+  # Restore distr and dynamic mask dimnames (lost in hdf5 write+read)
   ts_col <- ifelse(bf$metadata$ebird_version_year < 2022,
                    "interval",
                    "timestep"
-  )
-  dimnames(d) <- list(i = NULL, time = paste0("t", bf$dates[[ts_col]]))
-  bf$distr <- d
-
-  # Delete duplicated dynamic mask row
-  dm <- bf$geom$dynamic_mask
-  if (ncol(dm) == n_timesteps(bf) + 1) {
-    if (!all(dm[, 1] == dm[, ncol(dm)]))
-      stop("Expected first and last dynamic mask columns to matrch in circular",
-           "BirdFlow model")
-    dm <- dm[, 1:(ncol(dm) - 1)]
+  )      ### back compatibility code
+
+  # Set dimnames for distr and dynamic mask
+  dn <- list(i = NULL, time = paste0("t", bf$dates[[ts_col]]))
+  dimnames(bf$distr) <- dn
+  dimnames(bf$geom$dynamic_mask) <- dn
+
+  if (!is_fitted_model) {
+    # Need to import some stuff here
+    bf$marginals <- NULL
+    bf$distances <- as.numeric(h5read(hdf5, "distances"))
+    bf$uci <- h5read(hdf5, "uci", native = TRUE)
+    bf$lci <- h5read(hdf5, "lci", native = TRUE)
+    dimnames(bf$uci) <- dn
+    dimnames(bf$lci) <- dn
   }
 
-  dimnames(dm) <- list(i = NULL, time = paste0("t", bf$dates[[ts_col]]))
-  bf$geom$dynamic_mask <- dm
+  # Convert sparse matricies (if present) back into sparse matrices
+  # (only relevant if reimporting a previously imported and exported model)
+  if (has_marginals(bf) && bf$metadata$is_sparse) {
+    mn <- setdiff(names(bf$marginals), "index")
+    for (m in mn) {
+      bf$marginals[[m]] <- Matrix::Matrix(bf$marginals[[m]], sparse = TRUE)
+    }
+
+    # Clean up metadata$sparse attriubtes and order problems
+
+    sparse <- bf$metadata$sparse
+
+    # Restore standard order to sparse if names are as expected
+    sparse_order <-
+      c("fix_stats", "method", "arguments", "stats",
+        "pct_zero", "pct_density_lost")
+    if (setequal(sparse_order, names(sparse))) {
+      sparse <- sparse[sparse_order]
+    }
+
+    # Remove extra attributes hidden in data frame columns and vectors
+    for (i in seq_along(bf$metadata$sparse)) {
+      if (inherits(sparse[[i]], "data.frame")) {
+        sparse[[i]] <- clean_hdf5_dataframe(sparse[[i]])
+      } else {
+        sparse[[i]] <- as.vector(sparse[[i]])
+       }
+    }
+    # Cleanup extra attributes in argument list
+    if ("arguments" %in% names(sparse))
+      sparse$arguments <- lapply(sparse$arguments, as.vector)
 
-  # Make and save marginal index
-  bf$marginals$index <- make_marginal_index(bf)
+    bf$metadata$sparse <- sparse
+  } # end if sparse
 
   return(bf)
 }
diff --git a/R/new_birdflow.R b/R/new_birdflow.R
index ed6458e1..3e99d706 100644
--- a/R/new_birdflow.R
+++ b/R/new_birdflow.R
@@ -45,9 +45,11 @@ new_BirdFlow <- function() {
                    ebird_version_year = NA_integer_,
                    ebird_release_year = NA_integer_,
                    ebird_access_end_date = NA,
+                   ebirdst_version = NA,
                    birdflow_preprocess_date = NA,
                    birdflow_model_date = NA,
                    birdflow_version = 3,
+                   birdflowr_preprocess_version = NA,
                    birdflowr_version =
                      as.character(utils::packageVersion("BirdFlowR")),
                    is_sparse = FALSE,
diff --git a/R/preprocess_species.R b/R/preprocess_species.R
index b4767ed2..5344ad45 100644
--- a/R/preprocess_species.R
+++ b/R/preprocess_species.R
@@ -399,7 +399,7 @@ preprocess_species <- function(species = NULL,
   geom <- list(nrow = nrow(mask),
                ncol = ncol(mask),
                res = res(mask),
-               ext = as.vector(ext(mask)),
+               ext = as.numeric(as.vector(ext(mask))),
                crs = crs(mask),
                mask = NA)
   m <- terra::values(mask)
@@ -497,9 +497,9 @@ preprocess_species <- function(species = NULL,
     uci <- cbind(uci, uci[, 1, drop = FALSE])
     lci <- cbind(lci, lci[, 1, drop = FALSE])
     dynamic_mask <- cbind(dynamic_mask, dynamic_mask[, 1, drop = FALSE])
-    colnames(uci) <- colnames(lci) <-
-      colnames(distr) <- colnames(dynamic_mask) <-
-      paste0("t", seq_len(ncol(distr)))
+    dimnames(uci) <- dimnames(lci) <-
+      dimnames(distr) <- dimnames(dynamic_mask) <-
+      list(i = NULL, time = paste0("t", seq_len(ncol(distr))))
 
     export$distr <- distr
     export$uci <- uci
diff --git a/R/validate_BirdFlow.R b/R/validate_BirdFlow.R
index 50a1dcb1..c80e5c22 100644
--- a/R/validate_BirdFlow.R
+++ b/R/validate_BirdFlow.R
@@ -55,26 +55,34 @@ validate_BirdFlow <- function(x, error = TRUE, allow_incomplete = FALSE) {
   }
 
   report_problems <- function() {
-    # This function uses p, error,  and allow_incomplete
+    # Use p, error,  and allow_incomplete
     # from the parent function to either throw an error
     # OR return from the parent function.
     message <- NA
-    if (error) {
-      if (allow_incomplete) {
-        if (any(p$type == "error"))
-          message <- paste0("Problems found by validate_BirdFlow:\n\t",
-                            paste(p$problem[p$type == "error"],
-                                  collapse = ";\n\t"))
-      } else { # Don't allow incomplete:
-        if (nrow(p) > 0)
-          message <- paste0("Problems found by validate_BirdFlow:\n\t",
+
+    # If not throwing error return (visibly) the table
+    if (!error) {
+      do.call("return", args = list(p), envir = parent.frame(n = 1))
+    }
+
+
+    # If throwing an error compose message
+
+    if (allow_incomplete) {
+      if (any(p$type == "error"))
+        message <- paste0("Problems found by validate_BirdFlow:\n\t",
+                          paste(p$problem[p$type == "error"],
+                                collapse = ";\n\t"))
+    } else { # Don't allow incomplete:
+      if (nrow(p) > 0)
+        message <- paste0("Problems found by validate_BirdFlow:\n\t",
                           paste(p$problem, collapse = "; \n\t"))
-      }
     }
     if (!is.na(message))
       stop(message, call. =  FALSE)
 
-    do.call("return", args = list(p), envir = parent.frame(n = 1))
+    # error is TRUE but there are no errors: invisibly return empty df:
+    do.call("invisible", args = list(p), envir = parent.frame(n = 1))
 
   }
 
@@ -101,16 +109,17 @@ validate_BirdFlow <- function(x, error = TRUE, allow_incomplete = FALSE) {
     c("x$marginals should not be a list", # ok to have marginals
       "x$dates should not be a list", # ok to have dates
       "x$transitions should not be a list", # ok to have transitions
-      "x extra:uci, lci",  # Ok to have these (included in preprocessing output)
-      "x extra:lci, uci",  # ok in this order too
+      "x extra:uci",  # (included in preprocessing output)
+      "x extra:lci",  # (included in preprocessing output)
       "x missing:marginals", # ok to be missing marginals
-      "x missing:marginals, distances",
       "x missing:distances",
       "x$metadata$sparse_stats should not be a list",  # Having them is fine
       "x$metadata$sparse should not be a list",  #
-      "x$metadata missing:birdflow_version",
-      "x$metadata extra:hyperparameters",
-      "x$metadata extra:hyperparameters, loss_values"
+      "x$metadata missing:birdflow_version", ### back compatibility
+      "x$metadata extra:hyperparameters", # added by python
+      "x$metadata extra:loss_values", # added by python
+      "x$metadata missing:ebirdst_version", ### back compatibility
+      "x$metadata missing:birdflowr_preprocess_version" ### back compatibility
     ))
 
 
@@ -392,7 +401,7 @@ validate_BirdFlow <- function(x, error = TRUE, allow_incomplete = FALSE) {
       p <- add_prob("Not all marginals have a sum of one.", "error", p)
   }
 
-  return(report_problems())
+  report_problems()
 
 } # end validation function
 # nolint end
diff --git a/tests/testthat/_snaps/lookup_date.md b/tests/testthat/_snaps/lookup_date.md
new file mode 100644
index 00000000..9e390e95
--- /dev/null
+++ b/tests/testthat/_snaps/lookup_date.md
@@ -0,0 +1,7 @@
+# lookup dates works with transitions and marginals
+
+    Code
+      d
+    Output
+      [1] "2021-12-24" "2021-12-31" "2021-01-07" "2021-01-14"
+
diff --git a/tests/testthat/test-calc_flux.R b/tests/testthat/test-calc_flux.R
index d8d9f0c8..8ea0257c 100644
--- a/tests/testthat/test-calc_flux.R
+++ b/tests/testthat/test-calc_flux.R
@@ -1,5 +1,5 @@
 test_that("calc_flux() works without directionality", {
-
+  local_quiet()
   # Sparsify and truncate to speed things up
   bf <- BirdFlowModels::amewoo
   bf <- truncate_birdflow(bf, start = 1, end = 5)
@@ -20,7 +20,7 @@ test_that("calc_flux() works without directionality", {
 
 
 test_that("Test sensativity of flux to radius", {
-
+  local_quiet()
   testthat::skip("In depth flux radius analysis - always skipped")
 
   # This is an exploration of how flux varies across a very wide range
diff --git a/tests/testthat/test-extend_birdflow.R b/tests/testthat/test-extend_birdflow.R
index ee92fd2c..4e0eff09 100644
--- a/tests/testthat/test-extend_birdflow.R
+++ b/tests/testthat/test-extend_birdflow.R
@@ -17,6 +17,8 @@ test_that("extend_birdflow() works with hdf5", {
   skip_on_cran()
   skip_on_ci()
 
+  local_quiet()
+
   dir <- withr::local_tempdir("extend_hdf5")
 
   bf <- BirdFlowModels::amewoo
@@ -36,9 +38,15 @@ test_that("extend_birdflow() works with hdf5", {
   expect_true(extend_birdflow(hdf, e))
   bfe2 <- import_birdflow(hdf)
 
-  # For weird historical reasons:
+  ### Back compatibility code
+  ### (delete metadata items that only exist in some model versions)
   names(bfe1$metadata)[names(bfe1$metadata) == "birdFlowr_version"] <-
     "birdflowr_version"
+  bfe1$metadata$ebirdst_version <- NULL
+  bfe2$metadata$ebirdst_version <- NULL
+  bfe1$metadata$birdflowr_preprocess_version <- NULL
+  bfe2$metadata$birdflowr_preprocess_version <- NULL
+
 
    # Reimporting changes the import version so nuke both
   bfe1$metadata$birdflowr_version <- ""
diff --git a/tests/testthat/test-get_naturalearth.R b/tests/testthat/test-get_naturalearth.R
index bd0927cb..a22a6564 100644
--- a/tests/testthat/test-get_naturalearth.R
+++ b/tests/testthat/test-get_naturalearth.R
@@ -1,4 +1,5 @@
 test_that("get_coastline returns expected objects", {
+  local_quiet()
   bf <- new_BirdFlow()
   bf$geom$crs <- terra::crs(paste0("+proj=moll +lon_0=-90 +x_0=0 +y_0=0",
                             " +ellps=WGS84 +units=m +no_defs"))
@@ -34,7 +35,7 @@ test_that("get_coastline returns expected objects", {
 })
 
 test_that("get_countries returns expected objects", {
-
+  local_quiet()
   bf <- BirdFlowModels::amewoo
   expect_s3_class(countries <- get_countries(bf, scale = "small"),
                   class = c("sf", "data.frame"))
@@ -44,6 +45,7 @@ test_that("get_countries returns expected objects", {
 })
 
 test_that("get_states returns expected objects", {
+  local_quiet()
   # States requires rnaturalearthhires which is a large download
   # Use  devtools::install_github("ropensci/rnaturalearthhires")
   # to install.
@@ -66,6 +68,7 @@ test_that("get_states returns expected objects", {
 })
 
 test_that("get_naturalearth downloads and returns expected objects", {
+  local_quiet()
   # This downloads data and writes it to disk so skip everywhere but local
   # machine
   skip_on_ci()
@@ -86,7 +89,7 @@ test_that("get_naturalearth downloads and returns expected objects", {
 })
 
 test_that("get_naturalearth() works at edge of WGS84", {
-
+  local_quiet()
   # Construct a psuedo BirdFlow object that has a crs centered on the edge
   # of the wgs84 projection (used by rnaturalearth)
   # mollweide centered on 180 deg lon:
@@ -117,7 +120,7 @@ test_that("get_naturalearth() works at edge of WGS84", {
 
 
 test_that("get_naturalearth() works at edge of WGS84 with old method", {
-
+  local_quiet()
   # Construct a psuedo BirdFlow object that has a crs centered on the edge
   # of the wgs84 projection (used by rnaturalearth)
   # mollweide centered on 180 deg lon.
@@ -156,6 +159,7 @@ test_that("get_naturalearth() works at edge of WGS84 with old method", {
 })
 
 test_that("get_naturalearth() works at edge of WGS84 with double wrapping", {
+  local_quiet()
   # mollweide centered on 180 deg lon:
   seam_crs <- crs(paste0("+proj=moll +lon_0=180 +x_0=0 +y_0=0 +datum=WGS84 ",
                          "+units=m +no_defs +type=crs"))
@@ -172,6 +176,7 @@ test_that("get_naturalearth() works at edge of WGS84 with double wrapping", {
 })
 
 test_that("get_naturalearth() works with mollweide and broken bounding box", {
+  local_quiet()
   # Construct a psuedo BirdFlow object with the extent and
   #  projection that a user submitted.  This is a mollweide where
   # one corner of the bounding box is not on the map.
@@ -199,6 +204,7 @@ test_that("get_naturalearth() works with mollweide and broken bounding box", {
 
 
 test_that("get_naturalearth() works with lambert equal area (laea)", {
+  local_quiet()
   # Construct a psuedo BirdFlow object
   bf <- new_BirdFlow()
   bf$geom$crs <- crs(paste0("+proj=laea +lat_0=39.161 +lon_0=-85.094 +x_0=0 ",
@@ -218,7 +224,7 @@ test_that("get_naturalearth() works with lambert equal area (laea)", {
 })
 
 test_that("get_naturalearth() issues appropriate warning with empty extent", {
-
+  local_quiet()
   crs <- crs(paste0("+proj=moll +lon_0=180 +x_0=0 +y_0=0 +datum=WGS84 +units=m",
                     " +no_defs +type=crs")) # mollweide centered on 180 deg lon.
   bf <- new_BirdFlow()
@@ -250,6 +256,7 @@ test_that("get_naturalearth() issues appropriate warning with empty extent", {
 
 test_that(paste0("get_naturalearth with default keep_buffer = FALSE, crops ",
                  "to input extent"), {
+  local_quiet()
   bf <- new_BirdFlow()
   bf$geom$crs <- terra::crs(paste0("+proj=moll +lon_0=-90 +x_0=0 +y_0=0",
                                    " +ellps=WGS84 +units=m +no_defs"))
@@ -268,7 +275,7 @@ test_that(paste0("get_naturalearth with default keep_buffer = FALSE, crops ",
 
 
 test_that("get_naturalearth() works with projections that have +units=us-ft", {
-
+  local_quiet()
   bf <- new_BirdFlow()
   bf$geom$crs <- terra::crs("EPSG:2249")
   if (interactive()) {
@@ -298,6 +305,7 @@ test_that("get_naturalearth() works with projections that have +units=us-ft", {
 
 
 test_that("Double wrapped buffer works.", {
+  local_quiet()
   # Note this only sort of works because there are artifacts but
   # if you don't force the old method it works perfectly
   bf <- new_BirdFlow()
@@ -319,6 +327,7 @@ test_that("Double wrapped buffer works.", {
 
 
 test_that("Left wrapped buffer works.", {
+  local_quiet()
   bf <- new_BirdFlow()
   bf$geom$crs <- terra::crs(paste0("+proj=moll +lon_0=-90 +x_0=0 +y_0=0",
                                    " +ellps=WGS84 +units=m +no_defs"))
@@ -333,6 +342,7 @@ test_that("Left wrapped buffer works.", {
 
 
 test_that("Right wrapped buffer works.", {
+  local_quiet()
   bf <- new_BirdFlow()
   bf$geom$crs <- terra::crs(paste0("+proj=moll +lon_0=+90 +x_0=0 +y_0=0",
                                    " +ellps=WGS84 +units=m +no_defs"))
diff --git a/tests/testthat/test-import_birdflow.R b/tests/testthat/test-import_birdflow.R
new file mode 100644
index 00000000..d6116e10
--- /dev/null
+++ b/tests/testthat/test-import_birdflow.R
@@ -0,0 +1,48 @@
+test_that("import_birdflow() works with preprocessed species", {
+  local_quiet()
+  skip_on_cran()
+  dir <- withr::local_tempdir(pattern = "import_bf")
+  bf1 <- preprocess_species("example_data", hdf5 = TRUE, out_dir = dir,
+                            res = 200)
+  f <- list.files(dir, "\\.hdf5$", full.names = TRUE)[1]
+
+  bf2 <- import_birdflow(f)
+  expect_equal(bf1, bf2)
+
+  rexport <-  file.path(dir, "rexport.hdf5")
+  export_birdflow(bf2, rexport)
+
+  bf3 <- import_birdflow(rexport)
+  expect_equal(bf1, bf3)
+
+  expect_no_condition(validate_BirdFlow(bf3, allow_incomplete = TRUE))
+
+})
+
+test_that("export_birdflow() and import_birdflow() work with sparse models", {
+  local_quiet()
+  skip_on_cran()
+
+  bf <- BirdFlowModels::amewoo |> truncate_birdflow(start = 10, end = 15)
+  sbf <- sparsify(bf, method = "conditional", p = 0.99)
+
+  ### Back compatability, Needed for BirdFlowModels 0.0.2.9002
+  names(sbf$metadata)[names(sbf$metadata) == "birdFlowr_version"] <-
+    "birdflowr_version"
+
+
+  file <- withr::local_tempfile(fileext = ".hdf5")
+  expect_no_error(export_birdflow(sbf, file = file))
+
+  expect_no_error(sbf2 <- import_birdflow(file))
+
+  ### back compatibility these will have to be deleted when we update the
+  ### BirdFlowModels::amewoo model, Needed for BirdFlowModels 0.0.2.9002
+  sbf2$metadata$ebirdst_version <- NULL
+  sbf2$metadata$birdflowr_preprocess_version <- NULL
+
+
+  expect_equal(sbf, sbf2)
+
+
+})
diff --git a/tests/testthat/test-is_between.R b/tests/testthat/test-is_between.R
index f6b70feb..fe60cfd2 100644
--- a/tests/testthat/test-is_between.R
+++ b/tests/testthat/test-is_between.R
@@ -1,5 +1,7 @@
 test_that("is_between() works", {
 
+  local_quiet()
+
   # Sparsifying and truncating to speed things up
   bf <- BirdFlowModels::amewoo
   bf <- truncate_birdflow(bf, start = 1, end = 5)
diff --git a/tests/testthat/test-snap_to_birdflow.R b/tests/testthat/test-snap_to_birdflow.R
index 2adc07b9..87a1125f 100644
--- a/tests/testthat/test-snap_to_birdflow.R
+++ b/tests/testthat/test-snap_to_birdflow.R
@@ -55,15 +55,22 @@ test_that("snap_to_birdflow() works with preprocessed models", {
 
   local_quiet()  # to suppress preprocess chatter
 
-  # The required preprocesing takes a while
+  # The required preprocesing takes a long time:
   skip_on_cran()
   skip_on_covr()
   skip_on_ci()
 
-  bf <- preprocess_species("amewoo", res = 150, hdf5 = FALSE)
+  # And it requires and ebirdst access key
+  err <- tryCatch(bf <- preprocess_species("amewoo", res = 150, hdf5 = FALSE),
+                  error = identity)
+  if (inherits(err, "error") &&
+     grepl("Cannot access Status and Trends", err$message)) {
+    skip("Skipping - missing valid ebirdst key")
+  }
+
   d <- make_fake_move_data(bf)
   expect_no_error(
-    s <- snap_to_birdflow(d, x = "x", y = "y", crs = crs(bf),  bf = bf,
+    s <- snap_to_birdflow(d, x = "x", y = "y", crs = crs(bf), bf = bf,
                           id_cols = c("bird_id", "track_id"))
   )
 })