+ Show code- Hide code
library(OpasnetUtils)
library(xtable)
library(ggplot2)
library(reshape2)
objects.latest("Op_en6007", code_name = "answer") # [[OpasnetUtils/Drafts]]. We need ovaMatrixProduct.
limit <- 0.001
show_bins <- 10 # How many different direct inputs to show?
# Create the damage factor table based on data from [[Damage vector for life-cycle analysis]]
damagesPerImpact <- opbase.data("Op_en5902.damagefactors")
damagesPerImpact$Obs <- NULL
damagesPerImpact <- Ovariable("damagesPerImpact", damagesPerImpact)
# Take the impact factor table from the database. Do the same procedures as with damagesPerImpact.
impactsPerDollar <- Ovariable("impactsPerDollar", data = opbase.data("Op_en5902"))
# Take the direct requirements of an activity (in this case, producing a cup of coffee).
activity <- opbase.data("Op_en5902.coffee_cup_inputs")
activity$Obs <- NULL # Remove the redundant Obs column.
cat("Primary prosesses related to a cup of coffee (in Euro)\n")
head(activity)
# Combine the direct requirements of a coffee cup with a full vector of requirements and fill empty cells with 0.
activity <- merge(
unique(EvalOutput(impactsPerDollar)@output["Purchasing_sector"]),
activity,
all.x = TRUE
)
activity$Result[is.na(activity$Result)] <- 0
activity <- Ovariable("activity", data = activity)
# Get the normalisation data for different damages and make an ovariable out of it.
normalisation <- opbase.data("Op_en5904") # [[Normalisation data for life cycle assessments]]
normalisation$Obs <- NULL
head(normalisation)
normalisation <- Ovariable("normalisation", data = normalisation)
##########################333
damagesPerImpact <- EvalOutput(damagesPerImpact, N = 1)
impactsPerDollar <- EvalOutput(impactsPerDollar, N = 1)
activity <- EvalOutput(activity, N = 1)
normalisation <- EvalOutput(normalisation, N = 1)
#damagesPerImpact2 <- damagesPerImpact@output
#impactsPerDollar2 <- impactsPerDollar@output
#activity2 <- activity@output
#normalisation2 <- normalisation@output
#############################33
#damagesPerImpact2 <- reshape( # Reshape it into the wide format.
# damagesPerImpact2[ , colnames(damagesPerImpact2) != "Obs"], # Data to be reshaped
# times = "Result", # Column(s) to be reshaped
# timevar = "Unique_categories", # Column identifiers
# idvar = "Damage_categories", # Row identifiers
# direction = "wide" # Reshape from long to wide format
#)
#colnames(damagesPerImpact2) <- gsub("Result.", "", colnames(damagesPerImpact2)) # Remove extra "Result." from colnames.
#rownames(damagesPerImpact2) <- damagesPerImpact2[[1]] # Make the first column the rownames.
#damagesPerImpact2[1] <- NULL # Remove the first column.
#damagesPerImpact2[1:5, 1:5]
#damagesPerImpact2 <- t(as.matrix(damagesPerImpact2)) # Turn the data.frame into a matrix and transpose it.
#impactsPerDollar2 <- reshape(
# impactsPerDollar2,
# times = "Result",
# timevar = "Unique_categories",
# idvar = "Purchasing_sector",
# direction = "wide"
#)
#colnames(impactsPerDollar2) <- gsub("Result.", "", colnames(impactsPerDollar2))
#rownames(impactsPerDollar2) <- impactsPerDollar2[[1]]
#impactsPerDollar2[1] <- NULL
#impactsPerDollar2[1:5, 1:5]
#impactsPerDollar2 <- as.matrix(impactsPerDollar2)
#activity2 <- merge(data.frame(directRequirements2 = rownames(impactsPerDollar2)), activity2, all.x = TRUE)
#activity2$Result[is.na(activity2$Result)] <- 0
impactsPerDollar2 <- impactsPerDollar * activity # Multiply data matrix with activities.
head(impactsPerDollar@output)
head(damagesPerImpact@output)
out <- ovaMatrixProduct(
impactsPerDollar2,
damagesPerImpact,
c("Purchasing_sector", "Unique_categories"),
c("Unique_categories", "Damage_categories")
) # Do a matrix multiplication
# After matrix operations, turn out into a data.frame for graphics.
#out <- as.data.frame(out)
#out$Purchasing_sector <- rownames(out)
#out <- melt(out, idvars = "Purchasing_sector", variable_name = "damages") # Reshape into long format
# Merge all but show_bins largest bins to 'Other'.
sums <- as.data.frame(as.table(tapply(out@output$Result, out@output["Purchasing_sector"], sum)))
limit <- sort(sums$Freq, decreasing = TRUE)[show_bins]
keeps <- sums[sums$Freq > limit , "Purchasing_sector"]
keeps <- levels(keeps)[keeps]
levels(out@output$Purchasing_sector) <- ifelse(
levels(out@output$Purchasing_sector) %in% keeps,
levels(out@output$Purchasing_sector),
"Other"
)
#out <- dropall(out)
# Rename the columns to reflect actual things.
#colnames(out)[colnames(out) == "value"] <- "Result"
#damages <- EvalOutput(new("ovariable", name = "damages", data = out)) # Make an ovariable
damageFractions <- out / normalisation * 365 # Normalise and scale to daily values.
# Plot results on a bar graph.
cat("Effects smaller than or equal to ", limit / sum(sums$Freq) * 100, " % of the total effect are not shown.\n")
#oprint(damageFractions)
ggplot(damageFractions@output, aes(x = Damage, weight = Result, fill = Purchasing_sector)) + geom_bar() +
theme_grey(base_size = 18) +
theme(axis.text.x = element_text(angle = 45)) +
labs(
title = "Life cycle impacts of a cup of coffee",
x = "Damage",
y = "Amount"
)
| |