Skip to content

Commit

Permalink
Merge pull request #315 from USEPA/issue_303
Browse files Browse the repository at this point in the history
Update result object names
  • Loading branch information
bl-young authored Dec 5, 2024
2 parents d7958f8 + 07641f0 commit 5a6a62a
Show file tree
Hide file tree
Showing 5 changed files with 65 additions and 64 deletions.
83 changes: 42 additions & 41 deletions R/CalculationFunctions.R
Original file line number Diff line number Diff line change
Expand Up @@ -147,50 +147,51 @@ calculateResultsWithExternalFactors <- function(model, perspective = "FINAL", de
}

## Description of result components apply to both FINAL and DIRECT perspectives
# r1 - Domestic emissions from domestic production
# r2 - Emissions from imported goods consumed as intermediate products
# r3 - Emissions from imported goods consumed as final products
# see equations 4-7 and 9-12 in EPA 600/R-24/116
# G_d - Domestic emissions from domestic production
# G_mi - Emissions from imported goods consumed as intermediate products by domestic industries
# G_mf - Emissions from imported goods consumed as final products

if(perspective == "FINAL") {
# Calculate Final Perspective LCI (a matrix with total impacts in form of sector x flows)
logging::loginfo("Calculating Final Perspective LCI and LCIA with external import factors...")
subscript <- "f"
r1 <- model$B %*% model$L_d %*% diag(as.vector(y_d))
r2 <- model$M_m %*% model$A_m %*% model$L_d %*% diag(as.vector(y_d))
subscript <- "l"
G_d <- model$B %*% model$L_d %*% diag(as.vector(y_d))
G_mi <- model$M_m %*% model$A_m %*% model$L_d %*% diag(as.vector(y_d))

} else { # Calculate direct perspective results.
# Calculate Direct Perspective LCI (a matrix with total impacts in form of sector x flows)
logging::loginfo("Calculating Direct + Imported Perspective LCI and LCIA with external import factors...")
subscript <- "d"
subscript <- "r"
s <- getScalingVector(model$L_d, y_d)
r1 <- t(calculateDirectPerspectiveLCI(model$B, s))
r2 <- t(calculateDirectPerspectiveLCI(model$M_m, (model$A_m %*% model$L_d %*% y_d)))
G_d <- t(calculateDirectPerspectiveLCI(model$B, s))
G_mi <- t(calculateDirectPerspectiveLCI(model$M_m, (model$A_m %*% model$L_d %*% y_d)))
}
r3 <- model$M_m %*% diag(as.vector(y_m))
G_mf <- model$M_m %*% diag(as.vector(y_m))

if (use_domestic_requirements) {
# zero out the import results
r2[] <- 0
r3[] <- 0
G_mi[] <- 0
G_mf[] <- 0
}

if(show_RoW) {
if(model$specs$IODataSource=="stateior") {
# collapse third term for SoI and RoUS
r3 <- r3[, 1:sector_count] + r3[, (sector_count+1):(sector_count*2)]
G_mf <- G_mf[, 1:sector_count] + G_mf[, (sector_count+1):(sector_count*2)]

if(perspective == "DIRECT") {
# collapse second and third term for SoI and RoUS
r2 <- r2[, 1:sector_count] + r2[, (sector_count+1):(sector_count*2)]
G_mi <- G_mi[, 1:sector_count] + G_mi[, (sector_count+1):(sector_count*2)]
}
}
if(perspective == "DIRECT") {
LCI <- cbind(r1, r2 + r3) # Term 2 and Term 3 are assigned to RoW
LCI <- cbind(G_d, G_mi + G_mf) # Term 2 and Term 3 are assigned to RoW
} else {
LCI <- cbind(r1 + r2, r3) # Term 3 is assigned to RoW
LCI <- cbind(G_d + G_mi, G_mf) # Term 3 is assigned to RoW
}
} else {
LCI <- r1 + r2 + r3 # All three terms combined and regions do not change
LCI <- G_d + G_mi + G_mf # All three terms combined and regions do not change
}

# Calculate LCIA (matrix with direct impacts in form of sector x impacts)
Expand All @@ -210,8 +211,8 @@ calculateResultsWithExternalFactors <- function(model, perspective = "FINAL", de
LCI <- rbind(LCI, hh)
LCIA <- rbind(LCIA, hh_lcia)
}
result[[paste0("LCI_", subscript)]] <- LCI
result[[paste0("LCIA_", subscript)]] <- LCIA
result[[paste0("G_", subscript)]] <- LCI
result[[paste0("H_", subscript)]] <- LCIA
return(result)

}
Expand Down Expand Up @@ -255,24 +256,24 @@ calculateStandardResults <- function(model, perspective, demand, use_domestic_re
# Calculate Direct Perspective LCI (a matrix with direct impacts in form of sector x flows)
logging::loginfo("Calculating Direct Perspective LCI...")
s <- getScalingVector(L, f)
result$LCI_d <- calculateDirectPerspectiveLCI(model$B, s)
result$G_r <- calculateDirectPerspectiveLCI(model$B, s)
# Calculate Direct Perspective LCIA (matrix with direct impacts in form of sector x impacts)
logging::loginfo("Calculating Direct Perspective LCIA...")
result$LCIA_d <- calculateDirectPerspectiveLCIA(model$D, s)
result$H_r <- calculateDirectPerspectiveLCIA(model$D, s)
if(household_emissions) {
result$LCI_d <- rbind(result$LCI_d, hh)
result$LCIA_d <- rbind(result$LCIA_d, hh_lcia)
result$G_r <- rbind(result$G_r, hh)
result$H_r <- rbind(result$H_r, hh_lcia)
}
} else if (perspective=="FINAL") {
# Calculate Final Perspective LCI (a matrix with total impacts in form of sector x flows)
logging::loginfo("Calculating Final Perspective LCI...")
result$LCI_f <- calculateFinalPerspectiveLCI(M, f)
result$G_l <- calculateFinalPerspectiveLCI(M, f)
# Calculate Final Perspective LCIA (matrix with total impacts in form of sector x impacts)
logging::loginfo("Calculating Final Perspective LCIA...")
result$LCIA_f <- calculateFinalPerspectiveLCIA(N, f)
result$H_l <- calculateFinalPerspectiveLCIA(N, f)
if(household_emissions) {
result$LCI_f <- rbind(result$LCI_f, hh)
result$LCIA_f <- rbind(result$LCIA_f, hh_lcia)
result$G_l <- rbind(result$G_l, hh)
result$H_l <- rbind(result$H_l, hh_lcia)
}
}

Expand Down Expand Up @@ -304,9 +305,9 @@ getScalingVector <- function(L, demand) {
#' Journal of Cleaner Production 158 (August): 308–18. https://doi.org/10.1016/j.jclepro.2017.04.150.
#' SI1, Equation 9.
calculateDirectPerspectiveLCI <- function(B, s) {
lci_d <- t(B %*% diag(as.vector(s), nrow(s)))
rownames(lci_d) <- rownames(s)
return(lci_d)
G_r <- t(B %*% diag(as.vector(s), nrow(s)))
rownames(G_r) <- rownames(s)
return(G_r)
}

#' The final perspective LCI aligns flows with sectors consumed by final users.
Expand All @@ -319,10 +320,10 @@ calculateDirectPerspectiveLCI <- function(B, s) {
#' Journal of Cleaner Production 158 (August): 308–18. https://doi.org/10.1016/j.jclepro.2017.04.150.
#' SI1, Equation 10.
calculateFinalPerspectiveLCI <- function(M, y) {
lci_f <- t(M %*% diag(as.vector(y)))
colnames(lci_f) <- rownames(M)
rownames(lci_f) <- colnames(M)
return(lci_f)
G_l <- t(M %*% diag(as.vector(y)))
colnames(G_l) <- rownames(M)
rownames(G_l) <- colnames(M)
return(G_l)
}

#' The direct perspective LCIA aligns impacts with sectors consumed by direct use.
Expand All @@ -335,9 +336,9 @@ calculateFinalPerspectiveLCI <- function(M, y) {
#' Journal of Cleaner Production 158 (August): 308–18. https://doi.org/10.1016/j.jclepro.2017.04.150.
#' SI1, Equation 9.
calculateDirectPerspectiveLCIA <- function(D, s) {
lcia_d <- t(D %*% diag(as.vector(s), nrow(s)))
rownames(lcia_d) <- rownames(s)
return(lcia_d)
H_r <- t(D %*% diag(as.vector(s), nrow(s)))
rownames(H_r) <- rownames(s)
return(H_r)
}

#' The final perspective LCIA aligns impacts with sectors consumed by final users.
Expand All @@ -350,10 +351,10 @@ calculateDirectPerspectiveLCIA <- function(D, s) {
#' Journal of Cleaner Production 158 (August): 308–18. https://doi.org/10.1016/j.jclepro.2017.04.150.
#' SI1, Equation 10.
calculateFinalPerspectiveLCIA <- function(N, y) {
lcia_f <- t(N %*% diag(as.vector(y)))
colnames(lcia_f) <- rownames(N)
rownames(lcia_f) <- colnames(N)
return(lcia_f)
H_l <- t(N %*% diag(as.vector(y)))
colnames(H_l) <- rownames(N)
rownames(H_l) <- colnames(N)
return(H_l)
}

#' Divide/Normalize a sector x flows matrix by the total of respective flow (column sum)
Expand Down
4 changes: 2 additions & 2 deletions R/ValidateModel.R
Original file line number Diff line number Diff line change
Expand Up @@ -456,7 +456,7 @@ validateImportFactorsApproach <- function(model, demand = "Consumption"){
print(all.equal(LCI, LCI_dm))

# Calculate LCIA using standard approach
LCIA <- t(model$C %*% M %*% diag(as.vector(y))) #same as result_std_consumption$LCIA_f, above
LCIA <- t(model$C %*% M %*% diag(as.vector(y)))
colnames(LCIA) <- rownames(model$N_m)
rownames(LCIA) <- colnames(model$N_m)

Expand Down Expand Up @@ -491,7 +491,7 @@ validateHouseholdEmissions <- function(model) {
flows <- setNames(flows$FlowAmount, flows$Flow)

cat("\nTesting that LCI emissions from households are equivalent to calculated result from Total Consumption.\n")
result <- r$LCI_f[codes, names(flows)]
result <- r$G_l[codes, names(flows)]
all.equal(flows, result)
}

Expand Down
18 changes: 9 additions & 9 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -146,34 +146,34 @@ Note: `S00402/US - Used and secondhand goods` and `S00300/US - Noncomparable imp

#### Examples

Rank sectors based a composite score of selected total impacts (LCIA_d or LCIA_f) associated with total US demand (US production or consumption vector).
Rank sectors based a composite score of selected total impacts (H_r or H_l) associated with total US demand (US production or consumption vector).
Comparing rankings may also be used as another form of model validation that incorporates the demand vectors and the indicators as well as the model result matrices.

```r
# Calculate model LCIA_d and LCIA_f
# Calculate model LCIA_d (H_r) and LCIA_f (H_l)
result <- c(calculateEEIOModel(model, perspective = 'DIRECT', demand = "Production"),
calculateEEIOModel(model, perspective = 'FINAL', demand = "Consumption"))
colnames(result$LCIA_d) <- model$Indicators$meta[match(colnames(result$LCIA_d),
model$Indicators$meta$Name),
colnames(result$H_r) <- model$Indicators$meta[match(colnames(result$H_r),
model$Indicators$meta$Name),
"Code"]
colnames(result$LCIA_f) <- colnames(result$LCIA_d)
colnames(result$H_l) <- colnames(result$H_r)
# Define indicators
indicators <- c("ACID", "CCDD", "CMSW", "CRHW", "ENRG", "ETOX", "EUTR", "GHG",
"HRSP", "HTOX", "LAND", "MNRL", "OZON", "SMOG", "WATR")
# Create figure on the left
heatmapSectorRanking(model,
matrix = result$LCIA_d,
matrix = result$H_r,
indicators,
sector_to_remove = "",
N_sector = 20,
x_title = "LCIA_d (DIRECT perspective) & US production demand")
x_title = "H_r (DIRECT perspective) & US production demand")
# Create figure on the right
heatmapSectorRanking(model,
matrix = result$LCIA_f,
matrix = result$H_l,
indicators,
sector_to_remove = "",
N_sector = 20,
x_title = "LCIA_f (FINAL perspective) & US consumption demand")
x_title = "H_l (FINAL perspective) & US consumption demand")
```

![](inst/img/ranking_direct_prod_final_cons_v2.0.1.png)
Expand Down
10 changes: 5 additions & 5 deletions format_specs/Calculation.md
Original file line number Diff line number Diff line change
@@ -1,13 +1,13 @@
## Calculation Result
A result object from a Model calculation (`calculateEEIOModel()`) contains an LCI and an LCIA matrix. The matrices have a suffix appended to them in the form `_x`, where `x` indicates the calculation perspective that was used to prepare them.
A result object from a Model calculation (`calculateEEIOModel()`) contains an LCI (G) and an LCIA (H) matrix. The matrices have a suffix appended to them in the form `_x`, where `x` indicates the calculation perspective that was used to prepare them.

| Matrix | Shape | Description |
| --- | ---- | ---|
| LCI | sector x flow | Direct + indirect resource use and emissions by sector |
| LCIA | sector x indicator | Direct + indirect impacts by sector |
| G | sector x flow | Direct + indirect resource use and emissions by sector (LCI) |
| H | sector x indicator | Direct + indirect impacts by sector (LCIA) |

## Calculation Perspectives
|Suffix |Perspective|Definition|
|---|---|---|
|_d|DIRECT|Associates results with sector where emissions occur.|
|_f|FINAL|Associates results with sector is a final consumer. Final consumption is use by a final demand component.|
|_r|DIRECT|Associates results with sector where emissions occur.|
|_l|FINAL|Associates results with sector is a final consumer. Final consumption is use by a final demand component.|
14 changes: 7 additions & 7 deletions inst/doc/Example.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -48,18 +48,18 @@ print(paste("Sectors failing:", paste(econval$Failure$rownames, collapse = ", ")
```{r include=FALSE}
result <- c(useeior::calculateEEIOModel(model, perspective = 'DIRECT', demand = "Production"),
useeior::calculateEEIOModel(model, perspective = 'FINAL', demand = "Consumption"))
colnames(result$LCIA_d) <- model$Indicators$meta[match(colnames(result$LCIA_d), model$Indicators$meta$Name), "Code"]
colnames(result$LCIA_f) <- colnames(result$LCIA_d)
colnames(result$H_r) <- model$Indicators$meta[match(colnames(result$H_r), model$Indicators$meta$Name), "Code"]
colnames(result$H_l) <- colnames(result$H_r)
indicators <- c("ACID", "CCDD", "CMSW", "CRHW", "ENRG", "ETOX", "EUTR", "GHG",
"HRSP", "HTOX", "LAND", "MNRL", "OZON", "SMOG", "WATR")
model_list <- list("USEEIOv2.0.1-411" = model)
```

```{r "ranking_direct_prod_final_cons_v2.0.1", fig.width = 20, fig.height = 12}
p1 <- heatmapSectorRanking(model, matrix = result$LCIA_d, indicators,
sector_to_remove = "", N_sector = 20, x_title = "LCIA_d (DIRECT perspective) & US production demand")
p2 <- heatmapSectorRanking(model, matrix = result$LCIA_f, indicators,
sector_to_remove = "", N_sector = 20, x_title = "LCIA_f (FINAL perspective) & US consumption demand")
p1 <- heatmapSectorRanking(model, matrix = result$H_r, indicators,
sector_to_remove = "", N_sector = 20, x_title = "H_r (DIRECT perspective) & US production demand")
p2 <- heatmapSectorRanking(model, matrix = result$H_l, indicators,
sector_to_remove = "", N_sector = 20, x_title = "H_l (FINAL perspective) & US consumption demand")
gridExtra::grid.arrange(p1, p2, nrow = 1)
```

Expand All @@ -71,7 +71,7 @@ plotMatrixCoefficient(model_list, matrix_name = "N", coefficient_name = coeffs,
```{r "domestic_proportion_impact_USconsumption_v2.0.1", include=FALSE, fig.height = 12, fig.width = 12}
fullcons <- calculateEEIOModel(model, perspective = "DIRECT", demand = "Consumption", use_domestic_requirements = FALSE)
domcons <- calculateEEIOModel(model, perspective = "DIRECT", demand = "Consumption", use_domestic_requirements = TRUE)
barplotFloworImpactFractionbyRegion(R1_calc_result = domcons$LCIA_d, Total_calc_result = fullcons$LCIA_d, x_title = "")
barplotFloworImpactFractionbyRegion(R1_calc_result = domcons$H_r, Total_calc_result = fullcons$H_r, x_title = "")
```

```{r "indicator_score_v2.0.1", fig.height = 12, fig.width = 12}
Expand Down

0 comments on commit 5a6a62a

Please sign in to comment.