From f89edfdba024227981424696e84bac9a433fbc32 Mon Sep 17 00:00:00 2001 From: Eliot McIntire Date: Fri, 9 Nov 2018 17:38:32 -0800 Subject: [PATCH] who removes these? --- Boreal_LBMRDataPrep.R | 133 +++++++++++++++++++++--------------------- 1 file changed, 67 insertions(+), 66 deletions(-) diff --git a/Boreal_LBMRDataPrep.R b/Boreal_LBMRDataPrep.R index a1319b0..dd22310 100644 --- a/Boreal_LBMRDataPrep.R +++ b/Boreal_LBMRDataPrep.R @@ -105,7 +105,7 @@ 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) sim <- scheduleEvent(sim, P(sim)$.plotInitialTime, "Boreal_LBMRDataPrep", "plot") sim <- scheduleEvent(sim, P(sim)$.saveInitialTime, "Boreal_LBMRDataPrep", "save") @@ -133,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, @@ -148,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, @@ -164,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 @@ -181,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, @@ -190,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, @@ -199,7 +199,7 @@ estimateParameters <- function(sim) { pctCoverMinThresh = 50, userTags = "stable") if (ncell(sim$rasterToMatch) > 3e6) .gc() - + message("5: Derive Species Establishment Probability (SEP) from sim$specieslayers", Sys.time()) septable <- Cache(obtainSEP, ecoregionMap = simulationMaps$ecoregionMap, speciesLayers = sim$specieslayers, @@ -207,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, @@ -242,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, @@ -261,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) @@ -274,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", @@ -298,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) { @@ -311,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) @@ -351,19 +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)) } @@ -393,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") @@ -411,30 +411,31 @@ 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)) - + + browser() 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)) { @@ -453,7 +454,7 @@ Save <- function(sim) { 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$shpStudyArea, rasterToMatch = sim$rasterToMatch, @@ -463,7 +464,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 @@ -478,7 +479,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)) { @@ -488,24 +489,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, @@ -519,7 +520,7 @@ Save <- function(sim) { datatype = "INT2U", filename2 = TRUE, overwrite = TRUE, userTags = currentModule(sim)) - + projection(sim$LCC2005) <- projection(sim$rasterToMatch) } if (!suppliedElsewhere("ecoDistrict", sim)) { @@ -536,7 +537,7 @@ Save <- function(sim) { filename2 = TRUE, userTags = cacheTags) } - + if (!suppliedElsewhere("ecoRegion", sim)) { sim$ecoRegion <- Cache(prepInputs, targetFile = asPath(ecoregionFilename), @@ -551,7 +552,7 @@ Save <- function(sim) { filename2 = TRUE, userTags = cacheTags) } - + if (!suppliedElsewhere("ecoZone", sim)) { sim$ecoZone <- Cache(prepInputs, #notOlderThan = Sys.time(), targetFile = asPath(ecozoneFilename), @@ -566,7 +567,7 @@ Save <- function(sim) { filename2 = TRUE, userTags = cacheTags) } - + # stand age map if (!suppliedElsewhere("standAgeMap", sim)) { sim$standAgeMap <- Cache(prepInputs, #notOlderThan = Sys.time(), @@ -574,7 +575,7 @@ Save <- function(sim) { archive = asPath(c("kNN-StructureStandVolume.tar", "NFI_MODIS250m_kNN_Structure_Stand_Age_v0.zip")), destinationPath = dPath, - #url = extractURL("standAgeMap"), + url = extractURL("standAgeMap"), fun = "raster::raster", studyArea = sim$shpStudyAreaLarge, rasterToMatch = sim$rasterToMatch, @@ -583,7 +584,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( @@ -591,7 +592,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, @@ -603,13 +604,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, @@ -618,28 +619,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)) }