diff --git a/Boreal_LBMRDataPrep.R b/Boreal_LBMRDataPrep.R
index b802e8f..9765470 100644
--- a/Boreal_LBMRDataPrep.R
+++ b/Boreal_LBMRDataPrep.R
@@ -178,14 +178,19 @@ 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)
   ecoregionFiles <- Cache(ecoregionProducer,
-                          studyAreaRaster = initialCommFiles$initialCommunityMap,
-                          ecoregionMapFull = sim$ecoDistrict,
+                          #studyAreaRaster = initialCommFiles$initialCommunityMap,
+                          #ecoregionMapFull = sim$ecoDistrict,
+                          ecoregionMap = ecoregionMap,
                           ecoregionName = "ECODISTRIC",
                           ecoregionActiveStatus = ecoregionstatus,
-                          studyArea = sim$studyArea,
-                          rstStudyArea = rstStudyRegionBinary,
-                          maskFn = fastMask,
+                          rasterToMatch = initialCommFiles$initialCommunityMap, #sim$rasterToMatch,
+                          #studyArea = sim$studyArea,
+                          #rstStudyArea = rstStudyRegionBinary,
+                          #maskFn = fastMask,
                           userTags = "stable")
 
   message("3: ", Sys.time())
diff --git a/R/ecoregionProducers.R b/R/ecoregionProducers.R
index be1cb67..4289de1 100644
--- a/R/ecoregionProducers.R
+++ b/R/ecoregionProducers.R
@@ -1,48 +1,15 @@
-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, fixErrors(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)]))