diff --git a/Boreal_LBMRDataPrep.R b/Boreal_LBMRDataPrep.R index 25b1293..df0129b 100644 --- a/Boreal_LBMRDataPrep.R +++ b/Boreal_LBMRDataPrep.R @@ -72,7 +72,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 = ""), @@ -92,6 +92,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") ) @@ -102,9 +103,9 @@ defineModule(sim, list( doEvent.Boreal_LBMRDataPrep <- function(sim, eventTime, eventType, debug = FALSE) { if (eventType == "init") { - sim$specieslayers <- checkSpeciesLayerNames(sim$specieslayers) + names(sim$specieslayers) <- equivalentName(names(sim$specieslayers), sim$speciesEquivalency, "latinNames") sim <- estimateParameters(sim) - + # schedule future event(s) sim <- scheduleEvent(sim, P(sim)$.plotInitialTime, "Boreal_LBMRDataPrep", "plot") sim <- scheduleEvent(sim, P(sim)$.saveInitialTime, "Boreal_LBMRDataPrep", "save") @@ -124,35 +125,6 @@ doEvent.Boreal_LBMRDataPrep <- function(sim, eventTime, eventType, debug = FALSE ### template initialization -checkSpeciesLayerNames <- function(specieslayers) { - if (!all(names(specieslayers) %in% c("Pice_mar", "Pice_gla", "Abie_sp", "Pinu_sp", "Popu_tre", "Mixed"))) { - speciesEquivalencyTable <- data.table(Common = c("Black.Spruce", "White.Spruce", "Fir", "Pine", "Deciduous", "Mixed"), - Latin = c("Pice_mar", "Pice_gla", "Abie_sp", "Pinu_sp", "Popu_tre", "Mixed")) - matches <- pmatch(names(specieslayers), speciesEquivalencyTable$Common) - if (any(is.na(matches))) { - stop("specieslayers is expecting the species names to be Pice_mar, Pice_gla, Abie_sp, Pinu_sp, Popu_tre, Mixed.", - "Please rename them to these names") - } - names(specieslayers) <- speciesEquivalencyTable$Latin[matches] - } - - - if (FALSE) { - # not needed, but this is for calculating vegetation type maps, from a species abundances stack - VTM <- Cache(pemisc::makeVegTypeMap, sim$specieslayers, sim$vegLeadingProportion) - tableCC <- table(factorValues(CCvtm, CCvtm[], att = "Species")) - tablePP <- table(factorValues(VTM, VTM[], att = "Species")) - - propCC <- round(tableCC/sum(tableCC),2) - propPP <- round(tablePP/sum(tablePP),2) - names(propPP) <- speciesEquivalencyTable$Common[pmatch(names(propPP), speciesEquivalencyTable$Latin)] - - propCC[match(names(propPP), names(propCC))] - propPP - } - - return(specieslayers) -} estimateParameters <- function(sim) { # # ! ----- EDIT BELOW ----- ! # @@ -161,12 +133,12 @@ estimateParameters <- function(sim) { sim$ecoDistrict <- spTransform(sim$ecoDistrict, crs(sim$specieslayers)) sim$ecoRegion <- spTransform(sim$ecoRegion, crs(sim$specieslayers)) sim$ecoZone <- spTransform(sim$ecoZone, crs(sim$specieslayers)) - + message("1: ", Sys.time()) rstStudyRegionBinary <- raster(sim$rasterToMatch) rstStudyRegionBinary[] <- NA rstStudyRegionBinary[!is.na(sim$rasterToMatch[])] <- 1 - + message("2: ", Sys.time()) initialCommFiles <- Cache(initialCommunityProducer, speciesLayers = sim$specieslayers, @@ -176,7 +148,7 @@ estimateParameters <- function(sim) { userTags = "stable") ecoregionstatus <- data.table(active = "yes", 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, @@ -192,7 +164,7 @@ estimateParameters <- function(sim) { #rstStudyArea = rstStudyRegionBinary, #maskFn = fastMask, userTags = "stable") - + message("3: ", Sys.time()) # LCC05 -- land covers 1 to 15 are forested with tree dominated... 34 and 35 are recent burns # this is based on description in LCC05 @@ -209,7 +181,7 @@ estimateParameters <- function(sim) { "then try something like:\n", "reproducible::clearCache(userTags = c('LandWeb_dataPrep', 'init'), x = 'cache/SMALL_All')") } - + simulationMaps <- Cache(nonActiveEcoregionProducer, nonactiveRaster = sim$LCC2005, activeStatus = activeStatusTable, ecoregionMap = ecoregionFiles$ecoregionMap, @@ -218,7 +190,7 @@ estimateParameters <- function(sim) { initialCommunity = initialCommFiles$initialCommunity, userTags = "stable") if (ncell(sim$rasterToMatch) > 3e6) .gc() - + message("4: ", Sys.time()) speciesEcoregionTable <- Cache(obtainMaxBandANPP, speciesLayers = sim$specieslayers, biomassLayer = sim$biomassMap, @@ -227,8 +199,7 @@ estimateParameters <- function(sim) { pctCoverMinThresh = 50, userTags = "stable") if (ncell(sim$rasterToMatch) > 3e6) .gc() - - browser() + message("5: Derive Species Establishment Probability (SEP) from sim$specieslayers", Sys.time()) septable <- Cache(obtainSEP, ecoregionMap = simulationMaps$ecoregionMap, speciesLayers = sim$specieslayers, @@ -236,23 +207,23 @@ estimateParameters <- function(sim) { userTags = "stable") septable[, SEP := round(SEP, 4)] if (ncell(sim$rasterToMatch) > 3e6) .gc() - + message("6: ", Sys.time()) speciesEcoregionTable[, species := as.character(species)] septable[, species := as.character(species)] speciesEcoregionTable <- septable[speciesEcoregionTable, on = c("ecoregion", "species")] # speciesEcoregionTable <- left_join(speciesEcoregionTable, septable, by = c("ecoregion", "species")) %>% # data.table() - + # Fill in 0 for maxBiomass and maxANPP when SEP was estimated to be 0 speciesEcoregionTable[SEP == 0, ':='(maxBiomass = 0, maxANPP = 0)] NON_NAdata <- speciesEcoregionTable[!is.na(maxBiomass),] NAdata <- speciesEcoregionTable[is.na(maxBiomass),] - + if (nrow(NAdata) > 1) { # # replace NA values with ecoregion value #biomassFrombiggerMap <- sim$obtainMaxBandANPPFromBiggerEcoArea(speciesLayers = sim$specieslayers, - + message(" 6a obtainMaxBandANPPFromBiggerEcoArea: ", Sys.time()) biomassFrombiggerMap <- Cache(obtainMaxBandANPPFromBiggerEcoArea, speciesLayers = sim$specieslayers, @@ -271,7 +242,7 @@ estimateParameters <- function(sim) { NAdata <- biomassFrombiggerMap$addData[is.na(maxBiomass), .(ecoregion, species, maxBiomass, maxANPP, SEP)] } if (ncell(sim$rasterToMatch) > 3e6) .gc() - + message("7: ", Sys.time()) if (nrow(NAdata) > 1) { #biomassFrombiggerMap <- sim$obtainMaxBandANPPFromBiggerEcoArea(speciesLayers = sim$specieslayers, @@ -290,7 +261,7 @@ estimateParameters <- function(sim) { .(ecoregion, species, maxBiomass, maxANPP, SEP)] } if (ncell(sim$rasterToMatch) > 3e6) .gc() - + message("8: ", Sys.time()) NAdata[, ':='(maxBiomass = 0, maxANPP = 0, SEP = 0)] speciesEcoregion <- rbind(NON_NAdata, NAdata) @@ -303,15 +274,15 @@ estimateParameters <- function(sim) { sim$speciesEcoregion <- speciesEcoregion sim$ecoregion <- simulationMaps$ecoregion sim$ecoregionMap <- simulationMaps$ecoregionMap - + sim$initialCommunitiesMap <- Cache(createInitCommMap, simulationMaps$initialCommunityMap, as.integer(simulationMaps$initialCommunityMap[]), file.path(outputPath(sim), "initialCommunitiesMap.tif"), userTags = "stable") if (ncell(sim$rasterToMatch) > 3e6) .gc() - + message("9: ", Sys.time()) - + # species traits inputs speciesTable <- sim$speciesTable names(speciesTable) <- c("species", "Area", "longevity", "sexualmature", "shadetolerance", "firetolerance", @@ -327,12 +298,12 @@ estimateParameters <- function(sim) { "_", as.character(substring(species2, 1, 1)), tolower(as.character(substring(species2, 2, nchar(species2)))), sep = ""))] - + speciesTable$species <- toSentenceCase(speciesTable$species) speciesTable[species == "Pinu_con.con", species := "Pinu_con"] speciesTable[species == "Pinu_con.lat", species := "Pinu_con"] speciesTable[species == "Betu_all", species := "Betu_sp"] - + ## convert species names to match user-input list speciesList <- sim$speciesList rownames(speciesList) <- sapply(strsplit(speciesList[,1], "_"), function(x) { @@ -340,36 +311,36 @@ estimateParameters <- function(sim) { x[2] <- substring(x[2], 1, 3) paste(x, collapse = "_") }) - + ## replace eventual "spp" and "all" by sp (currently used instead of spp) rownames(speciesList) <- sub("_spp*", "_sp", rownames(speciesList)) rownames(speciesList) <- sub("_all", "_sp", rownames(speciesList)) - + ## match rownames to speciesTable$species rownames(speciesList) <- toSentenceCase(rownames(speciesList)) - + ## find matching names to replace in speciesTable matchNames <- speciesTable[species %in% rownames(speciesList), species] speciesTable[species %in% rownames(speciesList), species := speciesList[matchNames, 2]] - + ## filter table to existing species layers speciesTable <- speciesTable[species %in% names(sim$specieslayers)] - + 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"] - + initialCommunities <- simulationMaps$initialCommunity[, .(mapcode, description = NA, species)] set(initialCommunities, NULL, paste("age", 1:15, sep = ""), NA) initialCommunities <- data.frame(initialCommunities) message("11: ", Sys.time()) - + ## filter communities to species that have traits initialCommunities <- initialCommunities[initialCommunities$species %in% speciesTable$species,] - + initialCommunitiesFn <- function(initialCommunities, speciesTable) { for (i in 1:nrow(initialCommunities)) { agelength <- sample(1:15, 1) @@ -380,16 +351,19 @@ estimateParameters <- function(sim) { data.table::data.table(initialCommunities) } message("12: ", Sys.time()) - + sim$initialCommunities <- Cache(initialCommunitiesFn, initialCommunities, speciesTable, 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 ----- ! # return(invisible(sim)) } @@ -419,17 +393,17 @@ Save <- function(sim) { # ! ----- EDIT BELOW ----- ! # cPath <- cachePath(sim) dPath <- asPath(dataPath(sim), 1) - + # 1. test if all input objects are already present (e.g., from inputs, objects or another module) a <- depends(sim) whThisMod <- which(unlist(lapply(a@dependencies, function(x) x@name)) == "Boreal_LBMRDataPrep") objNames <- a@dependencies[[whThisMod]]@inputObjects$objectName objExists <- !unlist(lapply(objNames, function(x) is.null(sim[[x]]))) names(objExists) <- objNames - - + + crsUsed <- P(sim)[[".crsUsed"]] - + # Filenames ecoregionFilename <- file.path(dPath, "ecoregions.shp") ecodistrictFilename <- file.path(dPath, "ecodistricts.shp") @@ -437,30 +411,30 @@ Save <- function(sim) { biomassMapFilename <- file.path(dPath, "NFI_MODIS250m_kNN_Structure_Biomass_TotalLiveAboveGround_v0.tif") lcc2005Filename <- file.path(dPath, "LCC2005_V1_4a.tif") standAgeMapFilename <- file.path(dPath, "NFI_MODIS250m_kNN_Structure_Stand_Age_v0.tif") - + # Also extract fexts <- c("dbf", "prj", "sbn", "sbx", "shx") ecoregionAE <- basename(paste0(tools::file_path_sans_ext(ecoregionFilename), ".", fexts)) ecodistrictAE <- basename(paste0(tools::file_path_sans_ext(ecodistrictFilename), ".", fexts)) ecozoneAE <- basename(paste0(tools::file_path_sans_ext(ecozoneFilename), ".", fexts)) - + 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) - + seedToKeep <- .GlobalEnv$.Random.seed set.seed(1234) sim$shpStudyAreaLarge <- SpaDES.tools::randomPolygon(x = polyCenter, hectares = 10000) .GlobalEnv$.Random.seed <- seedToKeep } - + if (!suppliedElsewhere("shpStudyArea", sim)) { message("'shpStudyArea' was not provided by user. Using the same as 'shpStudyAreaLarge'") sim$shpStudyArea <- sim$shpStudyAreaLarge } - + needRTM <- FALSE if (is.null(sim$rasterToMatch)) { if (!suppliedElsewhere("rasterToMatch", sim)) { @@ -473,7 +447,7 @@ Save <- function(sim) { " or in a module that gets loaded prior to ", currentModule(sim)) } } - + if (!suppliedElsewhere("biomassMap", sim) || needRTM) { sim$biomassMap <- Cache(prepInputs, targetFile = asPath(basename(biomassMapFilename)), @@ -489,7 +463,7 @@ Save <- function(sim) { 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 @@ -504,7 +478,7 @@ Save <- function(sim) { } 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)) { @@ -514,24 +488,24 @@ Save <- function(sim) { 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, @@ -545,7 +519,7 @@ Save <- function(sim) { datatype = "INT2U", filename2 = TRUE, overwrite = TRUE, userTags = currentModule(sim)) - + projection(sim$LCC2005) <- projection(sim$rasterToMatch) } if (!suppliedElsewhere("ecoDistrict", sim)) { @@ -562,7 +536,7 @@ Save <- function(sim) { filename2 = TRUE, userTags = cacheTags) } - + if (!suppliedElsewhere("ecoRegion", sim)) { sim$ecoRegion <- Cache(prepInputs, targetFile = asPath(ecoregionFilename), @@ -577,7 +551,7 @@ Save <- function(sim) { filename2 = TRUE, userTags = cacheTags) } - + if (!suppliedElsewhere("ecoZone", sim)) { sim$ecoZone <- Cache(prepInputs, #notOlderThan = Sys.time(), targetFile = asPath(ecozoneFilename), @@ -592,7 +566,7 @@ Save <- function(sim) { filename2 = TRUE, userTags = cacheTags) } - + # stand age map if (!suppliedElsewhere("standAgeMap", sim)) { sim$standAgeMap <- Cache(prepInputs, #notOlderThan = Sys.time(), @@ -609,7 +583,7 @@ Save <- function(sim) { filename2 = TRUE, overwrite = TRUE, userTags = c("stable", currentModule(sim))) } - + if (!suppliedElsewhere("speciesList", sim)) { ## default to 6 species, one changing name, and two merged into one sim$speciesList <- as.matrix(data.frame( @@ -617,7 +591,7 @@ Save <- function(sim) { speciesNamesEnd = c("Abie_sp", "Pice_gla", "Pice_mar", "Pinu_sp", "Pinu_sp", "Popu_tre") )) } - + if (!suppliedElsewhere("specieslayers", sim)) { #opts <- options(reproducible.useCache = "overwrite") specieslayersList <- Cache(loadkNNSpeciesLayers, @@ -629,13 +603,13 @@ Save <- function(sim) { 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 } - + # 3. species maps sim$speciesTable <- Cache(prepInputs, "speciesTraits.csv", destinationPath = dPath, @@ -644,28 +618,28 @@ Save <- function(sim) { header = TRUE, stringsAsFactors = FALSE, userTags = c(cacheTags, "speciesTable")) %>% data.table() - + sim$sufficientLight <- data.frame(speciesshadetolerance = 1:5, X0 = 1, X1 = c(0.5, rep(1, 4)), X2 = c(0, 0.5, rep(1, 3)), X3 = c(rep(0, 2), 0.5, rep(1, 2)), X4 = c(rep(0, 3), 0.5, 1), X5 = c(rep(0, 4), 1)) - + if (!suppliedElsewhere("seedingAlgorithm", sim)) sim$seedingAlgorithm <- "wardDispersal" - + if (!suppliedElsewhere("successionTimestep", sim)) sim$successionTimestep <- 10 - + if (!suppliedElsewhere(sim$studyArea)) { sim$studyArea <- sim$shpStudyAreaLarge } - - + + if (!suppliedElsewhere("speciesThreshold", sim = sim)) { sim$speciesThreshold <- 50 } - + return(invisible(sim)) }