Skip to content

Commit

Permalink
Merge branch 'development' into LandRAssertionFix
Browse files Browse the repository at this point in the history
  • Loading branch information
@eliotmcintire
eliotmcintire authored Jun 6, 2024
2 parents 6f61feb + ec1f66e commit 29b1222
Showing 1 changed file with 107 additions and 70 deletions.
177 changes: 107 additions & 70 deletions Biomass_core.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -572,6 +572,7 @@ Init <- function(sim, verbose = getOption("LandR.verbose", TRUE)) {
## prepare species ------------------------------------------------
if (is.null(sim$species))
stop("'species' object must be provided")

species <- as.data.table(sim$species) # The former setDT actually changed the vector
LandR::assertSpeciesTable(species)
set(species, NULL, "speciesCode", factor(species$species, levels = unique(species$species))) # supply levels for speed
Expand Down Expand Up @@ -624,7 +625,7 @@ Init <- function(sim, verbose = getOption("LandR.verbose", TRUE)) {
rstLCC <- makeDummyRstLCC(sim$rasterToMatch)

## make sure speciesLayers match RTM (they may not if they come from another module's init.)
if (!.compareRas(sim$speciesLayers, sim$rasterToMatch)) {
if (!.compareRas(sim$speciesLayers, sim$rasterToMatch, stopOnError = FALSE)) {
message(blue("'speciesLayers' and 'rasterToMatch' do not match. "),
red("'speciesLayers' will be cropped/masked/reprojected to 'rasterToMatch'. "),
blue("If this is wrong, provide matching 'speciesLayers' and 'rasterToMatch'"))
Expand Down Expand Up @@ -780,6 +781,35 @@ Init <- function(sim, verbose = getOption("LandR.verbose", TRUE)) {
age = "integer",
B = "integer"))

## hamornize to simulated species -- we assume species is the correct set
## (it has been filtered by B_borealDP and B_speciesParams)
## leave sim$sppEquiv untouched for future reference - sim$sppColorVect can be changed
mod$sppEquiv <- sim$sppEquiv[get(P(sim)$sppEquivCol) %in% unique(sim$species$speciesCode)]
sppColorVect <- sim$sppColorVect[c(unique(as.character(sim$species$speciesCode)), "Mixed")]
sppColorVect <- sppColorVect[complete.cases(sppColorVect)]

sppOuts <- sppHarmonize(mod$sppEquiv, unique(sim$species$speciesCode), sppEquivCol = P(sim)$sppEquivCol,
sppColorVect = sppColorVect, vegLeadingProportion = P(sim)$vegLeadingProportion)

## TODO: it'd be great to functionize this:
if (length(setdiff(sim$sppColorVect, sppOuts$sppColorVect))) {
message(blue(
"sim$sppColorVect will be filtered to simulated species only (sim$species$speciesCode)"
))
}
sim$sppColorVect <- sppOuts$sppColorVect

if (length(setdiff(sim$sppNameVector, sppOuts$sppNameVector))) {
message(blue(
"sim$sppNameVector will be filtered to simulated species only (sim$species$speciesCode)"
))
}
sim$sppNameVector <- sppOuts$sppNameVector

assertSppVectors(sppEquiv = sim$species, sppEquivCol = "speciesCode",
sppColorVect = sim$sppColorVect)
assertSppVectors(sppEquiv = mod$sppEquiv, sppEquivCol = P(sim)$sppEquivCol, sppColorVect = sim$sppColorVect)

rasterNamesToCompare <- c("ecoregionMap", "pixelGroupMap")
if (!identical(P(sim)$initialBiomassSource, "cohortData")) {
rasterNamesToCompare <- c(rasterNamesToCompare, "biomassMap")
Expand Down Expand Up @@ -892,36 +922,38 @@ Init <- function(sim, verbose = getOption("LandR.verbose", TRUE)) {
sim$regenerationOutput <- data.table(seedingAlgorithm = character(), species = character(),
Year = numeric(), numberOfReg = numeric())
}
} else if (grepl("biomassMap", tolower(P(sim)$initialBiomassSource))) {
stop("'biomassMap as a value for P(sim)$initialBiomassSource is not working currently; ",
"please use 'cohortData'")
if (verbose > 0)
message("Skipping spinup and using the sim$biomassMap * SpeciesLayers pct as initial biomass values")
biomassTable <- data.table(biomass = as.vector(values(sim$biomassMap)),
pixelGroup = as.vector(values(pixelGroupMap)))
biomassTable <- na.omit(biomassTable)
maxBiomass <- maxValue(sim$biomassMap)
if (maxBiomass < 1e3) {
if (verbose > 0) {
message(crayon::green(" Because biomassMap values are all below 1000, assuming that these are on tonnes/ha.\n",
" Converting to $g/m^2$ by multiplying by 100"))
} else {
if (grepl("biomassMap", tolower(P(sim)$initialBiomassSource))) {
stop("'biomassMap as a value for P(sim)$initialBiomassSource is not working currently; ",
"please use 'cohortData'")
if (verbose > 0)
message("Skipping spinup and using the sim$biomassMap * SpeciesLayers pct as initial biomass values")
biomassTable <- data.table(biomass = as.vector(values(sim$biomassMap)),
pixelGroup = as.vector(values(pixelGroupMap)))
biomassTable <- na.omit(biomassTable)
maxBiomass <- maxValue(sim$biomassMap)
if (maxBiomass < 1e3) {
if (verbose > 0) {
message(crayon::green(" Because biomassMap values are all below 1000, assuming that these are on tonnes/ha.\n",
" Converting to $g/m^2$ by multiplying by 100"))
}
biomassTable[, `:=`(biomass = biomass * 100)]
}
biomassTable[, `:=`(biomass = biomass * 100)]
}

## In case there are non-identical biomasses in each pixelGroup -- this should be irrelevant with
## improved biomass_borealDataPrep.R (Jan 6, 2019 -- Eliot)
biomassTable <- biomassTable[, list(Bsum = mean(biomass, na.rm = TRUE)), by = pixelGroup]
if (!is.integer(biomassTable[["Bsum"]]))
set(biomassTable, NULL, "Bsum", asInteger(biomassTable[["Bsum"]]))

## Delete the B from cohortData -- it will be joined from biomassTable
set(cohortData, NULL, "B", NULL)
cohortData[, totalSpeciesPresence := sum(speciesPresence), by = "pixelGroup"]
cohortData <- cohortData[biomassTable, on = "pixelGroup"]
cohortData[, B := Bsum * speciesPresence / totalSpeciesPresence, by = c("pixelGroup", "speciesCode")]
if (!is.integer(cohortData[["B"]]))
set(cohortData, NULL, "B", asInteger(cohortData[["B"]]))
## In case there are non-identical biomasses in each pixelGroup -- this should be irrelevant with
## improved biomass_borealDataPrep.R (Jan 6, 2019 -- Eliot)
biomassTable <- biomassTable[, list(Bsum = mean(biomass, na.rm = TRUE)), by = pixelGroup]
if (!is.integer(biomassTable[["Bsum"]]))
set(biomassTable, NULL, "Bsum", asInteger(biomassTable[["Bsum"]]))

## Delete the B from cohortData -- it will be joined from biomassTable
set(cohortData, NULL, "B", NULL)
cohortData[, totalSpeciesPresence := sum(speciesPresence), by = "pixelGroup"]
cohortData <- cohortData[biomassTable, on = "pixelGroup"]
cohortData[, B := Bsum * speciesPresence / totalSpeciesPresence, by = c("pixelGroup", "speciesCode")]
if (!is.integer(cohortData[["B"]]))
set(cohortData, NULL, "B", asInteger(cohortData[["B"]]))
}
}

pixelAll <- cohortData[, .(uniqueSumB = sum(B, na.rm = TRUE)), by = pixelGroup]
Expand Down Expand Up @@ -965,7 +997,7 @@ Init <- function(sim, verbose = getOption("LandR.verbose", TRUE)) {
if (!is.null(P(sim)$calcSummaryBGM))
sim$vegTypeMap <- vegTypeMapGenerator(sim$cohortData, sim$pixelGroupMap,
P(sim)$vegLeadingProportion, mixedType = P(sim)$mixedType,
sppEquiv = sim$sppEquiv, sppEquivCol = P(sim)$sppEquivCol,
sppEquiv = mod$sppEquiv, sppEquivCol = P(sim)$sppEquivCol,
colors = sim$sppColorVect,
doAssertion = getOption("LandR.assertions", TRUE))

Expand Down Expand Up @@ -1064,7 +1096,7 @@ SummaryBGM <- compiler::cmpfun(function(sim) {
if (!is.null(P(sim)$calcSummaryBGM))
sim$vegTypeMap <- vegTypeMapGenerator(sim$cohortData, sim$pixelGroupMap,
P(sim)$vegLeadingProportion, mixedType = P(sim)$mixedType,
sppEquiv = sim$sppEquiv, sppEquivCol = P(sim)$sppEquivCol,
sppEquiv = mod$sppEquiv, sppEquivCol = P(sim)$sppEquivCol,
colors = sim$sppColorVect,
doAssertion = getOption("LandR.assertions", TRUE))

Expand Down Expand Up @@ -1098,8 +1130,8 @@ MortalityAndGrowth <- compiler::cmpfun(function(sim) {
## This tests for available memory and tries to scale the groupSize accordingly.
## It is, however, a very expensive operation. It now only does it once per simulation
groupSize <- maxRowsDT(maxLen = 1e7, maxMem = P(sim)$.maxMemory,
startClockTime = sim$._startClockTime, groupSize = groupSize,
modEnv = mod)
startClockTime = sim$._startClockTime, groupSize = groupSize,
modEnv = mod)

numGroups <- ceiling(length(pgs) / groupSize)
groupNames <- paste0("Group", seq(numGroups))
Expand Down Expand Up @@ -1631,8 +1663,8 @@ summaryRegen <- compiler::cmpfun(function(sim) {
})

plotSummaryBySpecies <- compiler::cmpfun(function(sim) {
LandR::assertSpeciesPlotLabels(sim$species$species, sim$sppEquiv)
assertSppVectors(sppEquiv = sim$sppEquiv, sppEquivCol = P(sim)$sppEquivCol, sppColorVect = sim$sppColorVect)
LandR::assertSpeciesPlotLabels(sim$species$species, mod$sppEquiv)
assertSppVectors(sppEquiv = mod$sppEquiv, sppEquivCol = P(sim)$sppEquivCol, sppColorVect = sim$sppColorVect)

## BIOMASS, WEIGHTED AVERAGE AGE, AVERAGE ANPP
## AND AGE OF OLDEST COHORT PER SPECIES
Expand Down Expand Up @@ -1677,11 +1709,11 @@ plotSummaryBySpecies <- compiler::cmpfun(function(sim) {
stringsAsFactors = FALSE)

whMixedLeading <- which(summaryBySpecies1$leadingType == "Mixed")
summaryBySpecies1$leadingType <- equivalentName(summaryBySpecies1$leadingType, sim$sppEquiv,
summaryBySpecies1$leadingType <- equivalentName(summaryBySpecies1$leadingType, mod$sppEquiv,
"EN_generic_short")
summaryBySpecies1$leadingType[whMixedLeading] <- "Mixed"

colours <- equivalentName(names(sim$sppColorVect), sim$sppEquiv, "EN_generic_short")
colours <- equivalentName(names(sim$sppColorVect), mod$sppEquiv, "EN_generic_short")
whMixedSppColors <- which(names(sim$sppColorVect) == "Mixed")
colours[whMixedSppColors] <- "Mixed"

Expand All @@ -1694,7 +1726,7 @@ plotSummaryBySpecies <- compiler::cmpfun(function(sim) {

if (length(unique(summaryBySpecies1$year)) > 1) {
df <- sim$species[, list(speciesCode, species)][summaryBySpecies, on = "speciesCode"]
df$species <- equivalentName(df$species, sim$sppEquiv, "EN_generic_short")
df$species <- equivalentName(df$species, mod$sppEquiv, "EN_generic_short")

colorIDs <- match(df$species, colours)
df$cols <- sim$sppColorVect[colorIDs]
Expand All @@ -1706,7 +1738,7 @@ plotSummaryBySpecies <- compiler::cmpfun(function(sim) {
unqCols2 <- unqdf$cols
names(unqCols2) <- unqdf$species

assertSppVectors(sppEquiv = sim$sppEquiv, sppEquivCol = "EN_generic_short", sppColorVect = unqCols2)
assertSppVectors(sppEquiv = mod$sppEquiv, sppEquivCol = "EN_generic_short", sppColorVect = unqCols2)

## although Plots can deal with .plotInitialTime == NA by not plotting, we need to
## make sure the plotting windows are not changed/opened if .plotInitialTime == NA
Expand Down Expand Up @@ -1784,8 +1816,8 @@ plotSummaryBySpecies <- compiler::cmpfun(function(sim) {
})

plotVegAttributesMaps <- compiler::cmpfun(function(sim) {
LandR::assertSpeciesPlotLabels(sim$species$species, sim$sppEquiv)
assertSppVectors(sppEquiv = sim$sppEquiv, sppEquivCol = P(sim)$sppEquivCol,
LandR::assertSpeciesPlotLabels(sim$species$species, mod$sppEquiv)
assertSppVectors(sppEquiv = mod$sppEquiv, sppEquivCol = P(sim)$sppEquivCol,
sppColorVect = sim$sppColorVect)

## these plots are not saved.
Expand Down Expand Up @@ -1813,7 +1845,7 @@ plotVegAttributesMaps <- compiler::cmpfun(function(sim) {
## Other species in levs[[levelsName]] are already "Leading",
## but it needs to be here in case it is not Leading in the future.
# The ones we want
sppEquiv <- sim$sppEquiv[!is.na(sim$sppEquiv[[P(sim)$sppEquivCol]]),]
sppEquiv <- mod$sppEquiv[!is.na(mod$sppEquiv[[P(sim)$sppEquivCol]]),]

levsLeading <- equivalentName(levs[[levelsName]], sppEquiv, "Leading")

Expand Down Expand Up @@ -1890,7 +1922,7 @@ plotVegAttributesMaps <- compiler::cmpfun(function(sim) {
})

plotAvgVegAttributes <- compiler::cmpfun(function(sim) {
LandR::assertSpeciesPlotLabels(sim$species$species, sim$sppEquiv)
LandR::assertSpeciesPlotLabels(sim$species$species, mod$sppEquiv)

## AVERAGE STAND BIOMASS/AGE/ANPP
## calculate acrosS pixels
Expand Down Expand Up @@ -2001,40 +2033,40 @@ CohortAgeReclassification <- function(sim) {
}

if (needRTM) {
if (is.null(sim$rawBiomassMap)) {
rawBiomassMapURL <- paste0("http://ftp.maps.canada.ca/pub/nrcan_rncan/Forests_Foret/",
if (is.null(sim$rawBiomassMap)) {
rawBiomassMapURL <- paste0("http://ftp.maps.canada.ca/pub/nrcan_rncan/Forests_Foret/",
"canada-forests-attributes_attributs-forests-canada/",
"2001-attributes_attributs-2001/",
"NFI_MODIS250m_2001_kNN_Structure_Biomass_TotalLiveAboveGround_v1.tif")

httr::with_config(config = httr::config(ssl_verifypeer = P(sim)$.sslVerify), {
rawBiomassMap <- prepRawBiomassMap(url = rawBiomassMapURL,
studyAreaName = P(sim)$.studyAreaName,
cacheTags = cacheTags,
to = sim$studyArea,
projectTo = NA, ## don't project to SA
destinationPath = dPath)
})
} else {
rawBiomassMap <- sim$rawBiomassMap
if (!.compareCRS(sim$rawBiomassMap, sim$studyArea)) {
## note that extents may never align if the resolution and projection do not allow for it
rawBiomassMap <- Cache(postProcess,
rawBiomassMap,
method = "bilinear",
to = sim$studyAreaLarge,
projectTo = NA, ## don't project to SA
overwrite = TRUE)
cacheTags = cacheTags,
to = sim$studyArea,
projectTo = NA, ## don't project to SA
destinationPath = dPath)
})
} else {
rawBiomassMap <- sim$rawBiomassMap
if (!.compareCRS(sim$rawBiomassMap, sim$studyArea)) {
## note that extents may never align if the resolution and projection do not allow for it
rawBiomassMap <- Cache(postProcess,
rawBiomassMap,
method = "bilinear",
to = sim$studyAreaLarge,
projectTo = NA, ## don't project to SA
overwrite = TRUE)
}
}
}

RTMs <- prepRasterToMatch(studyArea = sim$studyArea,
studyAreaLarge = sim$studyArea,
rasterToMatch = NULL,
rasterToMatchLarge = NULL,
destinationPath = dPath,
templateRas = rawBiomassMap,
studyAreaName = P(sim)$.studyAreaName,
RTMs <- prepRasterToMatch(studyArea = sim$studyArea,
studyAreaLarge = sim$studyArea,
rasterToMatch = NULL,
rasterToMatchLarge = NULL,
destinationPath = dPath,
templateRas = rawBiomassMap,
studyAreaName = P(sim)$.studyAreaName,
cacheTags = cacheTags)
sim$rasterToMatch <- RTMs$rasterToMatch
rm(RTMs)
Expand Down Expand Up @@ -2120,11 +2152,16 @@ CohortAgeReclassification <- function(sim) {
if (!suppliedElsewhere("species", sim)) {
speciesTable <- getSpeciesTable(dPath = dPath, url = extractURL("species"),
cacheTags = c(cacheTags, "speciesTable"))
tempSppEquiv <- if (!is.null(sim$speciesLayers)) { ## may be supplied *after* and not exist yet
sim$sppEquiv[get(P(sim)$sppEquivCol) %in% names(sim$speciesLayers)]
} else {
copy(sim$sppEquiv)
}
sim$species <- prepSpeciesTable(speciesTable = speciesTable,
# speciesLayers = sim$speciesLayers,
sppEquiv = sim$sppEquiv[get(P(sim)$sppEquivCol) %in%
names(sim$speciesLayers)],
sppEquiv = tempSppEquiv,
sppEquivCol = P(sim)$sppEquivCol)
rm(tempSppEquiv)
}

## if not using LandR growth/mortality drivers... (assumes LandR.CS)
Expand Down

0 comments on commit 29b1222

Please sign in to comment.