Merge remote-tracking branch 'upstream/development' into development

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

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,

R/helpers.R

@@ -29,4 +29,4 @@ fasterizeFromSp <- function(sp, raster, fieldName) {
     fasterize::fasterize(tempSf, raster)
   } else 
     fasterize::fasterize(tempSf, raster, field = fieldName)
-}
+}

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